C Count cubes in cubic lattice cut by sphere around origin
C A cube with an inscribed sphere is divided into n*n*n
C smaller cubes. The program counts how many of the small
C cubes are cut or touched at their vertices.
C Cubes only touched by the cutting sphere, but all points
C either inside or outside the sphere are counted only half,
C because they occur in inside/outside pairs.
C Author: Hugo Pfoertner http://www.pfoertner.org/
C Version history:
C 2003 Jul 16 Initial version
C 2016 Apr  9 Fixed overflow in comparison of radii
C
      IMPLICIT INTEGER (A-Z)
      PARAMETER (M=1000)
      INTEGER*8 MI, MA
      DO 100 R = 2,M
      RR = R * R
C L counts cubes fully cut
      L = 0
C N counts cubes touched at one of its vertices
      N = 0
C Loops over lattice points
      DO 110 I = -R, R-2, 2
      II = I*I
      IPP = (I+2)*(I+2)
      DO 120 J = -R, R-2, 2
      JJ = J*J
      JPP = (J+2)*(J+2)
      DO 130 K = -R, R-2, 2
      KK = K*K
      KPP = (K+2)*(K+2)
C Differences of squared radii of 8 vertices from cutting sphere
      D1 = II  + JJ  + KK  - RR
      D2 = II  + JPP + KK  - RR
      D3 = IPP + JJ  + KK  - RR
      D4 = IPP + JPP + KK  - RR
      D5 = II  + JJ  + KPP - RR
      D6 = II  + JPP + KPP - RR
      D7 = IPP + JJ  + KPP - RR
      D8 = IPP + JPP + KPP - RR
C      WRITE(*,*)R,I,J,D1,D2,D3,D4
      MI = MIN (D1,D2,D3,D4,D5,D6,D7,D8)
      MA = MAX (D1,D2,D3,D4,D5,D6,D7,D8)
C Check if sign change of squared differences occurs
      IF ( MI * MA .LE. 0 ) THEN
        IF ( MI * MA .EQ. 0 ) THEN
C Vertex hit
          N = N + 1
        ELSE
C Cube cut
          L = L + 1
        ENDIF
C      ELSE
      ENDIF
130   CONTINUE
120   CONTINUE
110   CONTINUE
C Only half of the vertex hits are counted
      WRITE (*,*) R, ',', L + N/2
100   CONTINUE
      END