C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      PROGRAM PRFFT
      PARAMETER( M =16 )
      PARAMETER( N = 2**M )
      PARAMETER( NHALF = 2**(M-1) )
      PARAMETER( NQRTR = 2**(M-2) )
      PARAMETER( ISGBIT = 16 )
      PARAMETER( IBASE = 2**ISGBIT )
      REAL*8 DATAR(N), DATAI(N)
      INTEGER FLAG
      INTEGER CKFLAG
      INTEGER PRIME
      INTEGER SYSOUT
      INTEGER SYSIN
      INTEGER RNDCHK
      REAL*8 CWR(NHALF,M-2), CWI(NHALF,M-2)
      COMMON /OMEGA/ CWR, CWI
      REAL*8 CHKSUM, CHKOUT, T1, T2
      COMMON /CKS/ CHKSUM, CHKOUT
C
      SYSIN = 5
      SYSOUT = 6
C
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C           INITIALIZE OMEGA TABLES
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      CALL KNUTHW
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C           CHECK VALIDITY OF ARRAY SIZES AND CHOICE OF BASE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      READ( SYSIN, 99992 ) RNDCHK
99992 FORMAT ( I6 )
      CALL CHECK( DATAR, DATAI, RNDCHK )
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C           GET MERSENNE NUMBER TO CHECK
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2000  CONTINUE
      READ( SYSIN, 99994 ) PRIME
99994 FORMAT ( I6 )
      IF( PRIME .LT. 3 ) STOP
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C           CHECK MERSENNE NUMBER AGAINST BOUNDS
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      IF ( PRIME .LT. NHALF*ISGBIT ) GO TO 3000
      WRITE( SYSOUT, 99995 ) PRIME, NHALF*ISGBIT-1
99995 FORMAT( 1X, I7, ' IS TOO LARGE, MAXIMUM IS ', I8 )
      GO TO 2000
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C           INITIALIZE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3000  WRITE( SYSOUT, 99996 ) PRIME
99996 FORMAT( ///,' PERFORMING LUCAS-LEHMER TEST ON 2 ** ', I6 )
      DO 4000 I = 1,N
            DATAR(I) = 0.0D0
            DATAI(I) = 0.0D0
4000  CONTINUE
      DATAR(N) = 4.0D0
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C           COMPUTE INITIAL CHECKSUM
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      CHKSUM = DATAR(N) * DATAR(N)
5000  CHKSUM = DINT( DMOD( CHKSUM + 0.5D0, DBLE(IBASE)) )
     |       + DINT( (CHKSUM + 0.5D0) / DBLE(IBASE) )
      IF ( CHKSUM .GE. DBLE(IBASE) ) GO TO 5000
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C           CRUNCH CRUNCH CRUNCH CRUNCH CRUNCH CRUNCH CRUNCH CRUNCH
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      CALL CLOCK(T1)
      DO 6000 LOOP = 1, PRIME-2
            CALL SQUARE( DATAR, DATAI )
            CALL CLEAN( DATAR, DATAI, CKFLAG )
            IF ( CKFLAG .NE. 0 ) GO TO 8000
            CALL LESS2( DATAR )
            CALL FOLD( DATAR, DATAI, PRIME, CKFLAG )
            IF ( CKFLAG .NE. 0 ) GO TO 9000
6000  CONTINUE
      CALL CLOCK(T2)
      CALL WINNER( DATAR, PRIME, FLAG )
      IF ( FLAG .EQ. 1 ) GOTO 7000
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C           YET ANOTHER LOSER
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      NCELL = 0
      RESIDUE = 0
140   NCELL = NCELL + 1
      RESIDUE = RESIDUE + DATAR(N-NCELL+1) * (2**(ISGBIT*(NCELL-1)))
      IF ( NCELL*ISGBIT .LT. 15 ) GOTO 140
141   FORMAT(' TESTING TIME WAS: ',G11.5,' SECONDS')
      WRITE( SYSOUT, 150 ) PRIME,IFIX(SNGL((RESIDUE)))
      WRITE( SYSOUT, 141) T2-T1
150   FORMAT( ' 2 ** ', I6, ' is NOT a MERSENNE prime & the RESIDUE
     $= ',O5)
      GOTO 2000

C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C           IMPOSSIBLE  WE'VE FOUND ONE.  BETTER CHECK THE HARDWARE.
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
7000  WRITE( SYSOUT, 99999 ) PRIME
      WRITE( SYSOUT, 141) T2-T1
99999 FORMAT( ' 2 ** ', I6, ' IS A MERSENNE PRIME' )
      GO TO 2000
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C           CHECKSUM FAILURE.  HARDWARE PROBLEM?
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
8000  WRITE( SYSOUT, 99900 ) CHKOUT, CHKSUM
99900 FORMAT(' CHECKSUM FAILURE, WAS ',D23.16,', SHOULD BE ',D23.16,/)
      STOP
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C           THIS IS SERIOUS.
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
9000  WRITE( SYSOUT, 99901 )
99901 FORMAT( ' FOLD ROUTINE IS SICK. CANNOT CONTINUE', / )
      STOP
      END
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      SUBROUTINE CHECK( DATAR, DATAI, RNDCHK )
      PARAMETER( M =16 )
      PARAMETER( N = 2**M )
      PARAMETER( NHALF = 2**(M-1) )
      PARAMETER( NQRTR = 2**(M-2) )
      PARAMETER( ISGBIT = 16 )
      PARAMETER( IBASE = 2**ISGBIT )
      REAL*8 DATAR(N), DATAI(N)
      INTEGER RNDCHK
      INTEGER CKFLAG
      INTEGER PRIME
      REAL*8 CARRY
      INTEGER SYSOUT
      REAL*8 CWR(NHALF,M-2), CWI(NHALF,M-2)
      COMMON /OMEGA/ CWR, CWI
      REAL*8 CHKSUM, CHKOUT
      COMMON /CKS/ CHKSUM, CHKOUT
C
      SYSOUT = 6
C
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     CHECK OF 'BASE-1' SQUARED
C           NOTE:          9 :=       81           F :=       E1
C           NOTE:         99 :=     9801          FF :=     FE01
C           NOTE:        999 :=   998001         FFF :=   FFE001
C           NOTE:       9999 := 99980001        FFFF := FFFE0001
C     OUR 'DIGIT' IS 'IBASE'-1
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      DO 1000 I = 1, NHALF
            DATAR(I) = 0.0D0
            DATAI(I) = 0.0D0
1000  CONTINUE
      DO 1100 I = NHALF+1, N
            DATAR(I) = DBLE(IBASE-1)
            DATAI(I) = 0.0D0
1100  CONTINUE
C
      CALL SQUARE( DATAR, DATAI )
C
      DATAR(N) = DBLE(DINT(DATAR(N)+0.5D0))
      IF ( DATAR(N) .NE. 0.0D0 ) GO TO 3000
      CARRY = 0.0D0
      DO 1200 I = N-1, 1, -1
            DATAR(I) = DBLE(DINT(DATAR(I)+0.5D0)) + CARRY
            CARRY = DINT( DATAR(I) / DBLE(IBASE) )
            DATAR(I+1) = DMOD( DATAR(I), DBLE(IBASE) )
1200  CONTINUE
      DATAR(1) = CARRY
C
      DO 1300 I = 1, NHALF-1
            IF ( DATAR(I) .NE. DBLE(IBASE) - 1.0D0 ) GO TO 3000
1300  CONTINUE
      IF ( DATAR(NHALF) .NE. DBLE(IBASE) - 2.0D0 ) GO TO 3000
      DO 1400 I = NHALF+1, N-1
            IF ( DATAR(I) .NE. 0.0D0 ) GO TO 3000
1400  CONTINUE
      IF ( DATAR(N) .NE. 1.0D0 ) GO TO 3000
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     SET UP FOR RANDOM CHECKS
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      DO 1500 I = 1, N
            DATAR(I) = 0.0D0
            DATAI(I) = 0.0D0
1500  CONTINUE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     RANDOMIZE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      CHKSUM = DINT( DBLE(RAN0(-1)))
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     WALT: THE PARAMETER TO RAN0, "132049", SHOULD BE DIFFERENT ON
C           EACH RUN.  SUGGEST YOU USE THE TIME OF DAY, IN SECONDS.
C           ADD 1 IN CASE SECONDS IS ZERO.
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      CHKSUM = DINT( DBLE(RAN0(132049)) * DBLE(IBASE) )
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     POPULATE WITH RANDOM NUMBERS
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      DATAR(NHALF+1) = DINT( DBLE(RAN0(1)) * DBLE(IBASE) )
      DATAR(NHALF+1) = DINT( DATAR(NHALF+1) / 2.0D0 )
      CHKSUM = DATAR(NHALF+1)
      DO 1600 I = NHALF+2, N
            DATAR(I) = DINT( DBLE(RAN0(1)) * DBLE(IBASE) )
            CHKSUM = CHKSUM + DATAR(I)
1600  CONTINUE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     CALCULATE CHECKSUM
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1700  CHKSUM = DINT( DMOD( CHKSUM + 0.5D0, DBLE(IBASE)) )
     |       + DINT( (CHKSUM + 0.5D0) / DBLE(IBASE) )
      IF ( CHKSUM .GE. DBLE(IBASE) ) GO TO 1700
C
      CHKSUM = CHKSUM * CHKSUM
C
1800  CHKSUM = DINT( DMOD( CHKSUM + 0.5D0, DBLE(IBASE)) )
     |       + DINT( (CHKSUM + 0.5D0) / DBLE(IBASE) )
      IF ( CHKSUM .GE. DBLE(IBASE) ) GO TO 1800
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     CALCULATE 'PRIME' TO AS HIGH AS POSSIBLE
C     IF TIMING, START CLOCK HERE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      PRIME = NHALF * ISGBIT - 1
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     DO RANDOM CHECKS
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      DO 1900 LOOP = 1, RNDCHK
            CALL SQUARE( DATAR, DATAI )
            CALL CLEAN( DATAR, DATAI, CKFLAG )
            IF ( CKFLAG .NE. 0 ) GO TO 2000
            CALL LESS2( DATAR )
            CALL FOLD( DATAR, DATAI, PRIME, CKFLAG )
            IF ( CKFLAG .NE. 0 ) GO TO 4000
1900  CONTINUE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     IF TIMING, STOP CLOCK HERE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      RETURN
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     CHECK ROUTINE FAILURE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2000  WRITE( SYSOUT, 99991 ) CHKSUM, CHKOUT
99991 FORMAT ( ' CHKSUM = ', D23.16, ',  CHKOUT = ', D23.16 )
C
3000  WRITE( SYSOUT, 99992 )
99992 FORMAT( ' CHECK ROUTINE FAILURE: DECREMENT ISGBIT PARAMETER', / )
      STOP
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     THIS IS SERIOUS.
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4000  WRITE( SYSOUT, 99993 )
99993 FORMAT( ' FOLD ROUTINE IS SICK. CANNOT CONTINUE', / )
      STOP
      END
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      SUBROUTINE SQUARE(DATAR,DATAI)
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     IF THIS ROUTINE DOES NOT VECTORIZE 100%, I THINK THAT I WILL
C     BE FORCED TO TAKE IT VERY SERIOUSLY AND KILL SOMEBODY.
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      PARAMETER( M =16 )
      PARAMETER( N = 2**M )
      PARAMETER( NHALF = 2**(M-1) )
      PARAMETER( NQRTR = 2**(M-2) )
      PARAMETER( ISGBIT = 16 )
      PARAMETER( IBASE = 2**ISGBIT )
      REAL*8 DATAR(N),DATAI(N)
      REAL*8 TDATAR(N),TDATAI(N)
      INTEGER PASS
      INTEGER EVNODD, EVEN
      REAL*8 A,B,C,D,E,F,G,H
      REAL*8 WR,WI
      REAL*8 TEMPR,TEMPI
      REAL*8 CWR(NHALF,M-2), CWI(NHALF,M-2)
      COMMON /OMEGA/ CWR, CWI
C
      EVEN = 0
      I = M / 2
      I = I + I
      EVNODD = M - I
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     FFT PASS 1 OF M, SPECIAL CASE
C           IMAGINARY PARTS KNOWN TO BE ZERO
C           DATAR(1:NHALF) KNOWN TO BE ZERO
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      PASS = 1
      K = 1
      DO 1000 J = 1, NHALF
            TEMPR = -DATAR(J+NHALF)
            TDATAR(K) = DATAR(J+NHALF)
            TDATAI(K) = 0.0D0
            TDATAR(K+1) = TEMPR*CWR(J,PASS)
            TDATAI(K+1) = TEMPR*CWI(J,PASS)
            K = K + 2
1000  CONTINUE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     FFT PASSES 2 THRU M-2 IF M IS ODD
C         ELSE PASSES 2 THRU M-3
C            THESE LOOPS, ALONG WITH THEIR IFFT COUNTERPARTS, WILL BURN
C            THE MAJORITY OF THE CPU TIME.
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      L = M-2
      IF ( EVNODD .EQ. EVEN ) L = M-3
      DO 1300 PASS = 2, L, 2
            K = 1
            DO 1100 J = 1, NHALF
                  WR = CWR(J,PASS)
                  WI = CWI(J,PASS)
                  TEMPR = TDATAR(J) - TDATAR(J+NHALF)
                  TEMPI = TDATAI(J) - TDATAI(J+NHALF)
                  DATAR(K) = TDATAR(J) + TDATAR(J+NHALF)
                  DATAI(K) = TDATAI(J) + TDATAI(J+NHALF)
                  DATAR(K+1) = TEMPR*WR - TEMPI*WI
                  DATAI(K+1) = TEMPR*WI + TEMPI*WR
                  K = K + 2
1100        CONTINUE
            K = 1
            DO 1200 J = 1, NHALF
                  WR = CWR(J,PASS+1)
                  WI = CWI(J,PASS+1)
                  TEMPR = DATAR(J) - DATAR(J+NHALF)
                  TEMPI = DATAI(J) - DATAI(J+NHALF)
                  TDATAR(K) = DATAR(J) + DATAR(J+NHALF)
                  TDATAI(K) = DATAI(J) + DATAI(J+NHALF)
                  TDATAR(K+1) = TEMPR*WR - TEMPI*WI
                  TDATAI(K+1) = TEMPR*WI + TEMPI*WR
                  K = K + 2
1200        CONTINUE
1300  CONTINUE
      IF ( EVNODD .NE. EVEN ) GO TO 1399
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     FFT PASS M-2, IF M IS EVEN
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      PASS = M-2
      K = 1
      DO 1305 J = 1, NHALF
            WR = CWR(J,PASS)
            WI = CWI(J,PASS)
            TEMPR = TDATAR(J) - TDATAR(J+NHALF)
            TEMPI = TDATAI(J) - TDATAI(J+NHALF)
            DATAR(K) = TDATAR(J) + TDATAR(J+NHALF)
            DATAI(K) = TDATAI(J) + TDATAI(J+NHALF)
            DATAR(K+1) = TEMPR*WR - TEMPI*WI
            DATAI(K+1) = TEMPR*WI + TEMPI*WR
            K = K + 2
1305  CONTINUE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     FFT PASS M-1 OF M, SPECIAL CASE, M IS EVEN
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      K = 1
      DO 1310 J = 1, NQRTR
            TDATAR(K)         = DATAR(J) + DATAR(J+NHALF)
            TDATAI(K)         = DATAI(J) + DATAI(J+NHALF)
            TDATAR(K+1)       = DATAR(J) - DATAR(J+NHALF)
            TDATAI(K+1)       = DATAI(J) - DATAI(J+NHALF)
            TDATAR(K+NHALF)   = DATAR(J+NQRTR) + DATAR(J+NQRTR+NHALF)
            TDATAI(K+NHALF)   = DATAI(J+NQRTR) + DATAI(J+NQRTR+NHALF)
            TDATAR(K+1+NHALF) = DATAI(J+NQRTR) - DATAI(J+NQRTR+NHALF)
            TDATAI(K+1+NHALF) = DATAR(J+NQRTR+NHALF) - DATAR(J+NQRTR)
            K = K + 2
1310  CONTINUE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     FFT PASS M OF M, M IS EVEN
C     COMBINED WITH SQUARING LOOP
C     COMBINED WITH IFFT PASS 1 OF M
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      DO 1315 J = 1, NHALF
            A = (TDATAR(J) + TDATAR(J+NHALF))
            B = (TDATAI(J) + TDATAI(J+NHALF))
            C = (TDATAR(J) - TDATAR(J+NHALF))
            D = (TDATAI(J) - TDATAI(J+NHALF))
            E = A*A - B*B
            F = A*B
            G = C*C - D*D
            H = C*D
            TDATAR(J) = (E+G) / N
            TDATAI(J) = (F+H) / NHALF
            TDATAR(J+NHALF) = (E-G) / N
            TDATAI(J+NHALF) = (F-H) / NHALF
1315  CONTINUE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     IFFT PASS 2 OF M, SPECIAL CASE, M IS EVEN
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      K = 1
      DO 1320 J = 1, NQRTR
            DATAR(J)             = TDATAR(K)    + TDATAR(K+1)
            DATAI(J)             = TDATAI(K)    + TDATAI(K+1)
            DATAR(J+NHALF)       = TDATAR(K)    - TDATAR(K+1)
            DATAI(J+NHALF)       = TDATAI(K)    - TDATAI(K+1)
            DATAR(J+NQRTR)       = TDATAR(K+NHALF) - TDATAI(K+NHALF+1)
            DATAI(J+NQRTR)       = TDATAI(K+NHALF) + TDATAR(K+NHALF+1)
            DATAR(J+NQRTR+NHALF) = TDATAR(K+NHALF) + TDATAI(K+NHALF+1)
            DATAI(J+NQRTR+NHALF) = TDATAI(K+NHALF) - TDATAR(K+NHALF+1)
            K = K + 2
1320  CONTINUE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     IFFT PASS 3, IF M IS EVEN
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      PASS = M-2
      K = 1
      DO 1325 J = 1, NHALF
            WR = CWR(J,PASS)
            WI = CWI(J,PASS)
            TEMPR=WR*DATAR(K+1)+WI*DATAI(K+1)
            TEMPI=WR*DATAI(K+1)-WI*DATAR(K+1)
            TDATAR(J)=DATAR(K)+TEMPR
            TDATAI(J)=DATAI(K)+TEMPI
            TDATAR(J+NHALF)=DATAR(K)-TEMPR
            TDATAI(J+NHALF)=DATAI(K)-TEMPI
            K = K + 2
1325  CONTINUE
      GO TO 1699
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     FFT PASS M-1 OF M, SPECIAL CASE, M IS ODD
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1399  K = 1
      DO 1400 J = 1, NQRTR
            DATAR(K)         = TDATAR(J) + TDATAR(J+NHALF)
            DATAI(K)         = TDATAI(J) + TDATAI(J+NHALF)
            DATAR(K+1)       = TDATAR(J) - TDATAR(J+NHALF)
            DATAI(K+1)       = TDATAI(J) - TDATAI(J+NHALF)
            DATAR(K+NHALF)   = TDATAR(J+NQRTR) + TDATAR(J+NQRTR+NHALF)
            DATAI(K+NHALF)   = TDATAI(J+NQRTR) + TDATAI(J+NQRTR+NHALF)
            DATAR(K+1+NHALF) = TDATAI(J+NQRTR) - TDATAI(J+NQRTR+NHALF)
            DATAI(K+1+NHALF) = TDATAR(J+NQRTR+NHALF) - TDATAR(J+NQRTR)
            K = K + 2
1400  CONTINUE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     FFT PASS M OF M, M IS ODD
C     COMBINED WITH SQUARING LOOP
C     COMBINED WITH IFFT PASS 1 OF M
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      DO 1500 J = 1, NHALF
            A = (DATAR(J) + DATAR(J+NHALF))
            B = (DATAI(J) + DATAI(J+NHALF))
            C = (DATAR(J) - DATAR(J+NHALF))
            D = (DATAI(J) - DATAI(J+NHALF))
            E = A*A - B*B
            F = A*B
            G = C*C - D*D
            H = C*D
            DATAR(J) = (E+G) / N
            DATAI(J) = (F+H) / NHALF
            DATAR(J+NHALF) = (E-G) / N
            DATAI(J+NHALF) = (F-H) / NHALF
1500  CONTINUE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     IFFT PASS 2 OF M, SPECIAL CASE, M IS EVEN
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      K = 1
      DO 1600 J = 1, NQRTR
            TDATAR(J)             = DATAR(K)    + DATAR(K+1)
            TDATAI(J)             = DATAI(K)    + DATAI(K+1)
            TDATAR(J+NHALF)       = DATAR(K)    - DATAR(K+1)
            TDATAI(J+NHALF)       = DATAI(K)    - DATAI(K+1)
            TDATAR(J+NQRTR)       = DATAR(K+NHALF) - DATAI(K+NHALF+1)
            TDATAI(J+NQRTR)       = DATAI(K+NHALF) + DATAR(K+NHALF+1)
            TDATAR(J+NQRTR+NHALF) = DATAR(K+NHALF) + DATAI(K+NHALF+1)
            TDATAI(J+NQRTR+NHALF) = DATAI(K+NHALF) - DATAR(K+NHALF+1)
            K = K + 2
1600  CONTINUE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     IFFT PASSES 3 THRU M-1, IF M IS ODD
C          ELSE PASSES 4 THRU M-1 IF M IS EVEN
C            THESE LOOPS, ALONG WITH THEIR FFT COUNTERPARTS, WILL BURN
C            THE MAJORITY OF THE CPU TIME.
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1699  L = M-2
      IF ( EVNODD .EQ. EVEN ) L = M-3
      DO 1900 PASS = L, 3, -2
            K = 1
            DO 1700 J = 1, NHALF
                  WR = CWR(J,PASS)
                  WI = CWI(J,PASS)
                  TEMPR=WR*TDATAR(K+1)+WI*TDATAI(K+1)
                  TEMPI=WR*TDATAI(K+1)-WI*TDATAR(K+1)
                  DATAR(J)=TDATAR(K)+TEMPR
                  DATAI(J)=TDATAI(K)+TEMPI
                  DATAR(J+NHALF)=TDATAR(K)-TEMPR
                  DATAI(J+NHALF)=TDATAI(K)-TEMPI
                  K = K + 2
1700        CONTINUE
            K = 1
            DO 1800 J = 1, NHALF
                  WR = CWR(J,PASS-1)
                  WI = CWI(J,PASS-1)
                  TEMPR=WR*DATAR(K+1)+WI*DATAI(K+1)
                  TEMPI=WR*DATAI(K+1)-WI*DATAR(K+1)
                  TDATAR(J)=DATAR(K)+TEMPR
                  TDATAI(J)=DATAI(K)+TEMPI
                  TDATAR(J+NHALF)=DATAR(K)-TEMPR
                  TDATAI(J+NHALF)=DATAI(K)-TEMPI
                  K = K + 2
1800        CONTINUE
1900  CONTINUE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     IFFT PASS M OF M, WE NO LONGER CARE ABOUT IMAGINARY PART
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      PASS = 1
      K = 1
      DO 2000 J = 1, NHALF
            WR = CWR(J,PASS)
            WI = CWI(J,PASS)
            TEMPR=WR*TDATAR(K+1)+WI*TDATAI(K+1)
            DATAR(J)=TDATAR(K)+TEMPR
            DATAR(J+NHALF)=TDATAR(K)-TEMPR
            K = K + 2
2000  CONTINUE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     DONE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      RETURN
      END
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     THE GOSPEL ACCORDING TO KNUTH
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      SUBROUTINE KNUTHW
      PARAMETER( M =16 )
      PARAMETER( N = 2**M )
      PARAMETER( NHALF = 2**(M-1) )
      PARAMETER( NQRTR = 2**(M-2) )
      PARAMETER( ISGBIT = 16 )
      PARAMETER( IBASE = 2**ISGBIT )
      REAL*16 QUADXR(M), QUADXI(M)
      REAL*16 QUADYR(NHALF), QUADYI(NHALF)
      INTEGER PASS
      INTEGER I
      INTEGER TWOS
      INTEGER Q
      REAL*8 CWR(NHALF,M-2), CWI(NHALF,M-2)
      COMMON /OMEGA/ CWR, CWI
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     THE SPEED OF THIS ENTIRE ROUTINE IS A "DON'T CARE"
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      QUADXR(1) = -1.0Q0
      QUADXI(1) = 0.0Q0
      DO 1000 I = 2, M
            QUADXR(I) = QSQRT( (1.0Q0 + QUADXR(I-1)) / 2 )
            QUADXI(I) = QSQRT( (1.0Q0 - QUADXR(I-1)) / 2 )
1000  CONTINUE
C
      QUADYR(1) = 1.0Q0
      QUADYI(1) = 0.0Q0
      DO 3000 I = 2, NHALF
            J = (I-1) * N
            TWOS = 0
2000        TWOS = TWOS + 1
            J = J / 2
            KK = J / 2
            LL = KK + KK
            IF ( J .EQ. LL ) GO TO 2000
            KK = (M+M) - TWOS
            Q = ((I-1)*N) / 2**TWOS 
            Q = (2**TWOS)*(Q-1) / N + 1
            QUADYR(I) = QUADYR(Q) * QUADXR(KK) - QUADYI(Q) * QUADXI(KK)
            QUADYI(I) = QUADYR(Q) * QUADXI(KK) + QUADYI(Q) * QUADXR(KK)
3000  CONTINUE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     WE MAY WANT TO TRY AN EXPERIMENT HERE AT A LATER DATE.
C     THIS LOOP MOVES THE QUAD PRECISION OMEGAS TO DOUBLE PRECISION.
C     KNUTH TALKS ABOUT TRUNCATING AS OPPOSED TO ROUNDING.  IT MAY
C     ADVENTAGEOUS TO TRUNCATE HERE.  DOES FORTRAN 77/SX HAVE SOME
C     TYPE OF INTRINSIC FUNCTION TO ACCOMPLISH TRUNCATION HERE?
C     IF NOT, WE SHOULD BE ABLE TO 'AND' OUT THE LOW ORDER BITS OF
C     THE MANTISSAS OF THE QUAD PRECISION VARIABLES.  THIS ENTIRE
C     EXERCISE MAY INCREASE OUR 'ISGBIT' PARAMETER SIGNIFICANTLY.
C     ON THE OTHER HAND, IT MIGHT BE A WASTE OF TIME.
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      DO 4000 I = 1, NHALF
            CWR(I,1) = QUADYR(I)
            CWI(I,1) = -QUADYI(I)
4000  CONTINUE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      L = 1
      JJ = NHALF
      DO 7000 PASS = 1, M-3
            L = L+L
            JJ = JJ / 2
            DO 6000 I = 1, JJ
                  DO 5000 J = 1, L
                        K  = I*L - L + J
                        KK = I*L - L + 1
                        CWR(K,PASS+1) = CWR(KK,PASS)
                        CWI(K,PASS+1) = CWI(KK,PASS)
5000              CONTINUE
6000        CONTINUE
7000  CONTINUE
C
      RETURN
      END
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      SUBROUTINE LESS2( DATAR )
      PARAMETER( M =16 )
      PARAMETER( N = 2**M )
      PARAMETER( NHALF = 2**(M-1) )
      PARAMETER( NQRTR = 2**(M-2) )
      PARAMETER( ISGBIT = 16 )
      PARAMETER( IBASE = 2**ISGBIT )
      REAL*8 DATAR( N )
C
      DO 3000 J = 1, 2
            I = N
1000        IF ( DATAR(I) .NE. 0.0D0 ) GO TO 2000
            DATAR(I) = DBLE(IBASE) - 1.0D0
            I = I - 1
            GO TO 1000
2000        DATAR(I) = DATAR(I) - 1.0D0
3000  CONTINUE
      RETURN
      END
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      SUBROUTINE WINNER( DATAR, PRIME, FLAG )
      PARAMETER( M =16 )
      PARAMETER( N = 2**M )
      PARAMETER( NHALF = 2**(M-1) )
      PARAMETER( NQRTR = 2**(M-2) )
      PARAMETER( ISGBIT = 16 )
      PARAMETER( IBASE = 2**ISGBIT )
      REAL*8 DATAR(N)
      INTEGER PRIME
      INTEGER TRULEN
      INTEGER FLAG
      REAL*8 MASK
C
      FLAG = 0
      TRULEN = PRIME / ISGBIT + 1
      DO 1000 I = N, N-TRULEN+2, -1
            IF ( DATAR(I) .NE. DBLE( IBASE - 1.0D0 ) ) RETURN
1000  CONTINUE
C
      J = PRIME / ISGBIT
      J = PRIME - J * ISGBIT
      MASK = 2**J
      IF ( DATAR(N-TRULEN+1 ) .NE. (MASK - 1.0D0 ) ) RETURN
      FLAG = 1
      RETURN
      END
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     FUNCTION RAN FROM 'NUMERICAL RECIPES'
C           EVIDENTLY, THE RANDOM NUMBER GENERATOR ON MY
C           68000 IF FUBAR.
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      FUNCTION RAN(ISEED)
      PARAMETER(IA=7141,IC=54773,IM=259200)
      ISEED=MOD(ISEED*IA+IC,IM)
      RAN=FLOAT(ISEED)/FLOAT(IM)
      RETURN
      END
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     FUNCTION RAN0 FROM 'NUMERICAL RECIPES'
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      FUNCTION RAN0(IDUM)
      DIMENSION V(97)
      DATA IFF /0/
      IF(IDUM.LT.0.OR.IFF.EQ.0)THEN
            IFF=1
            ISEED=ABS(IDUM)
            IDUM=1
            DO 11 J=1,97
                  DUM=RAN(ISEED)
11          CONTINUE
            DO 12 J=1,97
                  V(J)=RAN(ISEED)
12          CONTINUE
            Y=RAN(ISEED)
      ENDIF
      J=1+INT(97.*Y)
      IF(J.GT.97.OR.J.LT.1)PAUSE
      Y=V(J)
      RAN0=Y
      V(J)=RAN(ISEED)
      RETURN
      END
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      SUBROUTINE CLEAN( DATAR, DATAI, CKFLAG )
      PARAMETER( M =16 )
      PARAMETER( N = 2**M )
      PARAMETER( NHALF = 2**(M-1) )
      PARAMETER( NQRTR = 2**(M-2) )
      PARAMETER( ISGBIT = 16 )
      PARAMETER( IBASE = 2**ISGBIT )
      REAL*8 DATAR( N ), DATAI( N )
      REAL*8 CARRY( N )
      REAL*8 TEMP
      INTEGER CKFLAG
      REAL*8 CHKSUM, CHKOUT
      COMMON /CKS/ CHKSUM, CHKOUT
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     NOTE THAT A FUNNY THING HAS HAPPENED NOW THAT WE ARE SQUARING
C           USING FFT'S, INSTEAD OF SCHOOLBOY OR KARATSUBA.
C           OUR ENTIRE NUMBER IS 'SHIFTED' LEFT BY ONE WORD INTO THE
C           DATAR ARRAY.  IN OTHER WORDS, BEFORE CARRYING, OUR VALUE
C           NOW LIVES IN DATAR(1:N-1), WHEREAS BEFORE IT LIVED IN
C           DATAR(1:N).  NOTE THAT DATAR(N) SHOULD = ZERO.  PERHAPS
C           WE SHOULD PUT IN A CHECK AFTER THE FOLLOWING INSTRUCTION
C           THAT WILL CAUSE AN ABORT "IF ( DATAR(N) .NE. 0.0D0 )".
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      DATAR(N) = DINT( DATAR(N) + 0.5D0 )
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     COMPLEX CARRY PROPOGATION
C     STRIP OUT CARRIES, PASS 1
C          CARRIES ARE PUT IN DATAR
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      DATAI(1) = 0.0D0
      TEMP = DATAR(1) + 0.5D0 
      DATAI(2) = DINT(DMOD( TEMP, DBLE(IBASE) ))
      DATAR(1) = DINT( TEMP / DBLE(IBASE) )
      CARRY(N) = 0.0D0
      CARRY(N-1) = 0.0D0
      DO 1000 I = 2, N-1
            TEMP = DATAR(I) + 0.5D0 
            DATAI(I+1) = DINT(DMOD( TEMP, DBLE(IBASE) ))
            TEMP = DINT( TEMP / DBLE(IBASE) )
            DATAR(I) = DINT(DMOD( TEMP, DBLE(IBASE) ))
            CARRY(I-1) = DINT( TEMP / DBLE(IBASE) )
1000  CONTINUE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     ADD IN CARRIES, PASS 1
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      DO 2000 I = 1, N
            DATAR(I) = DATAI(I) + DATAR(I) + CARRY(I)
2000  CONTINUE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     SIMPLE CARRY PROPOGATION
C         WE NEED TO HIT THESE LOOPS AT LEAST ONCE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     CORRECTION
C     WE DON'T NEED TO HIT THESE LOOPS.
C         AT THIS POINT THE CARRIES ARE NEARLY COMPLETELY PROPOGATED.
C         BY A STROKE OF GOOD LUCK, THIS INCOMPLETENESS DOES NOT AFFECT
C         THE CHECKSUM  AND BY MORE GOOD LUCK, THE REST OF THE
C         PROPOGATION IS PICKED UP BY THE LOGIC IN ROUTINE FOLD
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3000  DO 4000 I = N, 2, -2
            TEMP = DATAR(I)
            DATAR(I) = DMOD( TEMP, DBLE(IBASE) )
            DATAR(I-1) = DATAR(I-1) + DINT( TEMP / DBLE(IBASE) )
4000  CONTINUE
      DO 5000 I = N-1, 3, -2
            TEMP = DATAR(I)
            DATAR(I) = DMOD( TEMP, DBLE(IBASE) )
            DATAR(I-1) = DATAR(I-1) + DINT( TEMP / DBLE(IBASE) )
5000  CONTINUE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     AS IN SIMILAR LOGIC IN ROUTINE FOLD, THE BRANCH IN DO LOOP 6000
C        IS EXCEEDEDINGLY RARELY TAKEN.   SO RARELY IN FACT THAT FOR
C        ALL PRACTICAL PURPOSES, WE CAN CONSIDER IT AS NEVER TAKEN.
C        HOWEVER, THERE IS ALWAYS A VERY SLIM POSSIBILITY THAT IT
C        WILL BE TAKEN, SO WE MUST LEAVE THE LOOP AS IS.
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      DO 6000 I = N, 2, -1
            IF ( DATAR(I) .GE. DBLE(IBASE) ) GO TO 3000
6000  CONTINUE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     CALCULATE CHECKSUM OF SQUARE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      CHKOUT = 0.0D0
      DO 8000 I = 1, N
            CHKOUT = CHKOUT + DATAR(I)
8000  CONTINUE
C
9000  CHKOUT = DINT( DMOD( CHKOUT + 0.5D0, DBLE(IBASE)) )
     |       + DINT( (CHKOUT + 0.5D0) / DBLE(IBASE) )
      IF ( CHKOUT .GE. DBLE(IBASE) ) GO TO 9000
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     RETURN 'OK' IF CHECKSUMS MATCH
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      CKFLAG = 0
      IF ( CHKSUM .EQ. CHKOUT ) RETURN
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     CHECKSUM FAILURE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      CKFLAG = 1
      RETURN
      END
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      SUBROUTINE FOLD( DATAR, DATAI, PRIME, FLAG )
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     BEFORE SQUARING WE HAD A NUMBER 'PRIME' BITS LONG.  WE NOW HAVE
C           A NUMBER '2*PRIME' BITS LONG.  TO ACCOMPLISH MODULAR
C           DIVISION BY 2**PRIME-1,  BREAK THE STRING OF 2*PRIME BITS
C           INTO TWO STRINGS OF 'PRIME' BITS EACH.  THEN ADD THE TWO
C           RESULTANT STRINGS TOGETHER, BEING VERY AWARE THAT A SINGLE
C           BIT 'CARRY' CAN OCCUR.  IF WE DO HAVE A CARRY, ADD IT BACK
C           IN, BUT THIS CAN HAPPEN ONLY ONCE.  IF IT HAPPENS TWICE,
C           SOMETHING IS SERIOUSLY WRONG.
C     DATAI IS USED AS A SCRATCH ARRAY IN THIS ROUTINE.
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      PARAMETER( M =16 )
      PARAMETER( N = 2**M )
      PARAMETER( NHALF = 2**(M-1) )
      PARAMETER( NQRTR = 2**(M-2) )
      PARAMETER( ISGBIT = 16 )
      PARAMETER( IBASE = 2**ISGBIT )
      REAL*8 DATAR(N), DATAI(N)
      INTEGER PRIME
      INTEGER TRULEN
      INTEGER FLAG
      REAL*8 MASK
      REAL*8 SFTLFT
      REAL*8 TEMP
      REAL*8 CHKSUM, CHKOUT
      COMMON /CKS/ CHKSUM, CHKOUT
C
      FLAG = 0
      TRULEN = PRIME / ISGBIT + 1
      J = PRIME / ISGBIT
      J = PRIME - J * ISGBIT
      MASK = 2**J
      I = ISGBIT-J
      SFTLFT = 2**I
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     RE-STRUCTURED CODE SO THAT LOOP DO 1000 COULD BE ELIMINATED
C     MERGED LOOPS DO 1100 AND DO 1200 TOGETHER
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      DO 1100 I = 1, TRULEN
            DATAI(NHALF+1-I) = DINT( DATAR( N-TRULEN+2-I ) / MASK )
     |                       + DMOD( DATAR( N-TRULEN+1-I), MASK )
     |                       * SFTLFT
1100  CONTINUE
      DATAR( N+1-TRULEN) = DMOD(DATAR(N+1-TRULEN), MASK )
      DO 1300 I = NHALF, NHALF-TRULEN +1, -1
            DATAR(NHALF+I) = DATAR( NHALF+I) + DATAI(I)
1300  CONTINUE

1400  L = N-TRULEN+1
      DO 1500 I = N, L, -2
            TEMP = DATAR(I)
            DATAR(I) = DMOD( TEMP, DBLE(IBASE) )
            DATAR(I-1) = DATAR(I-1) + DINT( TEMP / DBLE(IBASE) )
1500  CONTINUE
      DO 1510 I = N-1, L, -2
            TEMP = DATAR(I)
            DATAR(I) = DMOD( TEMP, DBLE(IBASE) )
            DATAR(I-1) = DATAR(I-1) + DINT( TEMP / DBLE(IBASE) )
1510  CONTINUE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     LUCKY FOR US, IT SEEMS LIKE THE BRANCH INSIDE THE FOLLOWING LOOP
C           IS TAKEN VERY RARELY.  BUT WE STILL NEED TO CHECK ANYWAY.
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      DO 1520 I = N, N-TRULEN+1, -1
            IF ( DATAR(I) .GE. DBLE(IBASE) ) GO TO 1400
1520  CONTINUE
      IF ( DATAR( N-TRULEN+ 1 ) .LT. MASK ) GO TO 1800
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     WE HAVE A SINGLE BIT CARRY.
C     THIS LOOP DOES NOT NEED TO VECTORIZE AS IT WILL ALMOST ALWAYS
C          BAIL OUT IMMEDIATELY.
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      DATAR( N-TRULEN+1 ) = DATAR( N-TRULEN+1 ) - MASK
      DO 1700 I = N, N-TRULEN +2, -1
            IF ( DATAR(I) .EQ. DBLE(IBASE-1) ) GO TO 1600
            DATAR(I) = DATAR(I) + 1.0D0
            GO TO 1800
1600        DATAR(I) = 0.0D0
1700  CONTINUE
      DATAR( N-TRULEN+1 ) = DATAR( N-TRULEN+1 ) + 1.0D0
      IF ( DATAR( N-TRULEN+1 ) .LT. MASK ) GO TO 1800
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     WE NEED TO CARRY TWICE.  SERIOUS FLAW IN PROGRAM.  GO DIE.
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      FLAG = 1
      RETURN
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     NO NEED TO ZERO OUT DATAR1:NHALF! AS THE VERY FIRST LOOP IN
C          SQUARE ASSUMES THAT PART TO BE ZERO.
C          LOOP CONTROL WAS --> 1800  DO 1900 I = 1, N - TRULEN
C     ZERO OUT HIGH ORDER PART OF REAL ARRAY
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1800  DO 1900 I = NHALF+1, N - TRULEN
            DATAR(I) = 0.0D0
1900  CONTINUE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     IT IS OK TO ELIMINATE THIS LOOP BECAUSE THE VERY FIRST LOOP
C          IN THE FFT SQUARE ROUTINE ASSUMES THAT THE IMAGINARY PART
C          IS ALREADY ZERO.  IF SAID LOOP CHANGES, THEN IT MAY BE
C          NECCESSARY TO RE-INSTATE THE FOLLOWING DO 2000 LOOP.
C     RE-ZERO IMAGINARY PART
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C      DO 2000 I = 1,N
C            DATAI(I) = 0.0D0
C2000  CONTINUE
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C     CALCULATE NEW CHECKSUM
C     THE LOOP CONTROL USED TO BE -->  DO 2100 I = NHALF+1, N
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      CHKSUM = 0.0D0
      DO 2100 I = N-TRULEN+1, N
            CHKSUM = CHKSUM + DATAR(I)
2100  CONTINUE
C
2200  CHKSUM = DINT( DMOD( CHKSUM + 0.5D0, DBLE(IBASE)) )
     |       + DINT( (CHKSUM + 0.5D0) / DBLE(IBASE) )
      IF ( CHKSUM .GE. DBLE(IBASE) ) GO TO 2200
C
      CHKSUM = CHKSUM * CHKSUM
C
2300  CHKSUM = DINT( DMOD( CHKSUM + 0.5D0, DBLE(IBASE)) )
     |       + DINT( (CHKSUM + 0.5D0) / DBLE(IBASE) )
      IF ( CHKSUM .GE. DBLE(IBASE) ) GO TO 2300
      RETURN
      END

