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