C Construction of A071383, A071384 and A071385
C derived from John Conway's analysis of my
C brute force results in the seqfan mailing list
C
C Author: Hugo Pfoertner, http://www.pfoertner.org/
C
C Change history:
C 11.06.2002 Check with formula for number of points on circle
C 09.06.2002 11 primes, avoid exponentiation, extended to
C            cover all combinations up to 2**63
C            (range of signed 64 bit integers)
C 06.06.2002 Overflow protection for P=89
C 05.06.2002 Extended to 10 primes of the form 4*N+1
C 03.06.2002 Initial version
C
      IMPLICIT INTEGER*8 (A-Z)
      PARAMETER (NPRIM=11)
      DIMENSION P(NPRIM), Q(NPRIM), QQ(NPRIM)
      DOUBLEPRECISION LMAX, L(NPRIM), LP(NPRIM)
      EXTERNAL ICIRLT
C
C Primes of the form 4*N+1
      DATA P /5, 13, 17, 29, 37, 41, 53, 61, 73, 89, 97/
C
C Logarithms for overflow protection
C (this limits the range to be scanned)
      LMAX = LOG ( 2.0D0**63 - 1.0D0 )
C      LMAX = LOG ( 2.0D0**31 - 1.0D0 )
C
      DO 5 I = 1, NPRIM
      LP(I) = LOG ( DBLE ( P(I) ) )
      L(I) = 0.0D0
      Q(I) = 1
5     CONTINUE
C
C Loops over all products P_k**E_k with nonincreasing E_k
C Range of first loop given by 64 bit limitation,
C because 5**28~=2**65.01, whereas 5**27~=2**62.69
      DO 10 I1 = 1, 27
      L(1) = L(1) + LP(1)
C This condition truncates the processing when the required range
C is exceeded (similar for subsequent inner loops)
      IF ( L(1) .GT. LMAX ) GOTO 10
      Q(1) = Q(1) * P(1)
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
      DO 20 I2 = 0, I1
C Loop index 0 contributes the factor 1 to the product
C (similar for subsequent inner loops)
      IF ( I2 .GT. 0 ) THEN
        Q(2) = Q(2) * P(2)
        L(2) = L(2) + LP(2)
        IF ( L(1) + L(2) .GT. LMAX ) GOTO 20
      ENDIF
      QQ(2) = Q(1) * Q(2) 
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
      DO 30 I3 = 0, I2
      IF ( I3 .GT. 0 ) THEN
        Q(3) = Q(3)*P(3)
        L(3) = L(3) + LP(3)
        IF ( L(1) + L(2) + L(3) .GT. LMAX ) GOTO 30
      ENDIF
      QQ(3) = QQ(2) * Q(3) 
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
      DO 40 I4 = 0, I3
      IF ( I4 .GT. 0 ) THEN
        Q(4) = Q(4)*P(4)
        L(4) = L(4) + LP(4)
        IF ( L(1)+L(2)+L(3)+L(4) .GT. LMAX ) GOTO 40
      ENDIF
      QQ(4) = QQ(3) * Q(4)
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
      DO 50 I5 = 0, I4
      IF ( I5 .GT. 0 ) THEN
        Q(5) = Q(5) * P(5)
        L(5) = L(5) + LP(5)
        IF ( L(1)+L(2)+L(3)+L(4)+L(5) .GT. LMAX ) GOTO 50
      ENDIF
      QQ(5) = QQ(4) * Q(5)
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
      DO 60 I6 = 0, I5
      IF ( I6 .GT. 0 ) THEN
        Q(6) = Q(6) * P(6)
        L(6) = LP(6) * DBLE(I6)
        IF ( L(1)+L(2)+L(3)+L(4)+L(5)+L(6) .GT. LMAX ) GOTO 60
      ENDIF
      QQ(6) = QQ(5) * Q(6)
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
      DO 70 I7 = 0, I6
      IF ( I7 .GT. 0 ) THEN
        Q(7) = Q(7) * P(7)
        L(7) = L(7) + LP(7)
        IF ( L(1)+L(2)+L(3)+L(4)+L(5)+L(6)+L(7) .GT. LMAX ) GOTO 70
      ENDIF
      QQ(7) = QQ(6) * Q(7)
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
      DO 80 I8 = 0, I7
      IF ( I8 .GT. 0 ) THEN
        Q(8) = Q(8) * P(8)
        L(8) = L(8) + LP(8)
        IF ( L(1)+L(2)+L(3)+L(4)+L(5)+L(6)+L(7)+L(8)
     &     .GT. LMAX ) GOTO 80
      ENDIF
      QQ(8) = QQ(7) * Q(8)
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
      DO 90 I9 = 0, I8
      IF ( I9 .GT. 0 ) THEN
        Q(9) = Q(9) * P(9)
        L(9) = L(9) + LP(9)
        IF ( L(1)+L(2)+L(3)+L(4)+L(5)+L(6)+L(7)+L(8)+L(9)
     &     .GT. LMAX ) GOTO 90
      ENDIF
      QQ(9) = QQ(8) * Q(9)
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
C I10 has to be limited to 9, because 89^10~=2^64.76 > 2^63
      DO 100 I10 = 0, MIN(9,I9)
      IF ( I10 .GT. 0 ) THEN
        Q(10) = Q(10) * P(10)
        L(10) = L(10) + LP(10)
        IF ( L(1)+L(2)+L(3)+L(4)+L(5)+L(6)+L(7)+L(8)+L(9)+L(10)
     &     .GT. LMAX ) GOTO 100
      ENDIF
      QQ(10) = QQ(9) * Q(10)
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
C I11 has to be limited to 9, because 97^10~=2^66.00 > 2^63
C This limitation is automatically achieved by the previous limitation
C of I10 (I11<=I10)
      DO 110 I11 = 0, I10
      IF ( I11 .GT. 0 ) THEN
        Q(11) = Q(11) * P(11)
        L(11) = L(11) + LP(11)
      ENDIF
C Diagnostic output
C**      WRITE (*,*) (Q(K),K=1,NPRIM)
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
      IF ( L(1)+L(2)+L(3)+L(4)+L(5)+L(6)+L(7)+L(8)+L(9)+L(10)+L(11)
     & .LE. LMAX ) THEN
C final result:
C        Z = Q(1)*Q(2)*Q(3)*Q(4)*Q(5)*Q(6)*Q(7)*Q(8)*Q(9)*Q(10)*Q(11)
         Z = QQ(10) * Q(11)
C KIDNNT takes a REAL*8 argument and returns an INTEGER*8 result
C on DIGITAL Fortran, the generic Fortran didn't work
C        SQZ = NINT (SQRT ( DBLE ( Z ) ) )
        SQZ = KIDNNT (SQRT ( DBLE ( Z ) ) )
        IF ( SQZ*SQZ .LT. Z ) SQZ = SQZ + 1
C Function ICIRLT counts lattice points hit by circle with given
C squared radius (the hard way)
         M = ICIRLT ( Z )
C John Conway found
C        M = 4 * (I1+1)*(I2+1)*(I3+1)*(I4+1)*(I5+1)*(I6+1)
C     &         *(I7+1)*(I8+1)*(I9+1)*(I10+1)*(I11+1)
C Output to terminal
        WRITE (*,1000) I1, I2, I3, I4, I5, I6, I7, I8, I9, I10, I11,
     &                 Z, SQZ, M
C Output to file on channel 11
        WRITE (11,1000) I1, I2, I3, I4, I5, I6, I7, I8, I9, I10, I11,
     &                 Z, SQZ, M
      ENDIF
1000  FORMAT ( 11I3, I20, I11, I7 )
C
110   CONTINUE
C Reset products and log sums
      Q(11) = 1
      L(11) = 0.0D0
100   CONTINUE
      Q(10) = 1
      L(10) = 0.0D0
90    CONTINUE
      Q(9) = 1
      L(9) = 0.0D0
80    CONTINUE
      Q(8) = 1
      L(8) = 0.0D0
70    CONTINUE
      Q(7) = 1
      L(7) = 0.0D0
60    CONTINUE
      Q(6) = 1
      L(6) = 0.0D0
50    CONTINUE
      Q(5) = 1
      L(5) = 0.0D0
40    CONTINUE
      Q(4) = 1
      L(4) = 0.0D0
30    CONTINUE
      Q(3) = 1
      L(3) = 0.0D0
20    CONTINUE
      Q(2) = 1
      L(2) = 0.0D0
10    CONTINUE
      END
