C Test program for primality test in Fortran 77 C C written by Hugo Pfoertner http://www.pfoertner.org/ C C Expected result: Number of primes < 10^n C See C http://www.research.att.com/projects/OEIS?Anum=A006880 C 10^1 2 3 4 5 6 7 8 9 C [4,25,168,1229,9592,78498,664579,5761455,50847534] C C Test ISPRIME LOGICAL ISPRIM EXTERNAL ISPRIM L = 4 N10 = 10 DO 10 I = 11, 1000000001, 2 IF ( I .GT. N10 ) THEN WRITE (*,*) N10, L N10 = N10 * 10 ENDIF IF ( ISPRIM(I) ) L = L + 1 10 CONTINUE END C C...+....1....+....2....+....3....+....4....+....5....+....6....+....7.. C Adaptation of prime.cpp by Olivier Langlois C http://www3.sympatico.ca/michel.langlois2/archive/prime_cpp.htm C Change history: C Jul 3 2003 Initial version C...+....1....+....2....+....3....+....4....+....5....+....6....+....7.. C LOGICAL FUNCTION ISPRIM ( N ) INTEGER N INTEGER SP (1:10) LOGICAL SPRP EXTERNAL SPRP C C LIST OF SMALL PRIMES DATA SP / 2, 3, 5, 7, 11, 13, 17, 19, 23, 29 / ISPRIM = .FALSE. C C LOOK UP FIRST 10 PRIMES C Test may be deactivated of only primes > 29 shall be tested C Observed CPU time savings: ~=10% for primes in the range 2*10^9 IF ( N. LE. 29 ) THEN DO 5 I = 1, 10 IF ( N .EQ. SP(I) ) THEN ISPRIM = .TRUE. RETURN ENDIF 5 CONTINUE ENDIF C C EXCLUDE MULTIPLES OF FIRST 10 PRIMES (inline) C Multiples of 2 C Test may be deactivated if only odd numbers are to be tested. C Observed savings in CPU time are negligible (-0.5%) when C bit test function BTEST is used. IF ( .NOT. BTEST (N,0) ) RETURN C if BTEST is not available, replace by: C IF ( MOD ( N, 2 ) .EQ. 0 ) RETURN C The observed saving in CPU time by using BTEST instead of C MOD(...,2) was -12% for primes in the range 2*10^9. This was C the cumulative effect of using BTEST also in the subsequent C auxiliary functions. C C Next small primes IF ( MOD ( N, 3 ) .EQ. 0 ) RETURN IF ( MOD ( N, 5 ) .EQ. 0 ) RETURN IF ( MOD ( N, 7 ) .EQ. 0 ) RETURN IF ( MOD ( N, 11 ) .EQ. 0 ) RETURN IF ( MOD ( N, 13 ) .EQ. 0 ) RETURN IF ( MOD ( N, 17 ) .EQ. 0 ) RETURN IF ( MOD ( N, 19 ) .EQ. 0 ) RETURN IF ( MOD ( N, 23 ) .EQ. 0 ) RETURN IF ( MOD ( N, 29 ) .EQ. 0 ) RETURN C ISPRIM = .TRUE. C Range already covered by sieving IF ( N .LT. 853 ) RETURN C C The following part is a modified version of the original C code C using the optimized parameters found by Jaeschke C see: Finding primes & proving primality at C http://www.utm.edu/research/primes/prove/prove2_3.html C ISPRIM = .FALSE. IF ( N .LT. 9080191 ) THEN IF ( .NOT. SPRP(N,31) ) RETURN IF ( .NOT. SPRP(N,73) ) RETURN ISPRIM = .TRUE. ELSE C Valid for N < 4_759_123_141 (sufficient to cover unsigned 32bit) IF ( .NOT. SPRP(N,2) ) RETURN IF ( .NOT. SPRP(N,7) ) RETURN IF ( .NOT. SPRP(N,61) ) RETURN ISPRIM = .TRUE. ENDIF END C...+....1....+....2....+....3....+....4....+....5....+....6....+....7.. LOGICAL FUNCTION SPRP ( N, B ) C Translation of b_SPRP, Strong probable Prime test from prime.cpp C by Olivier Langlois INTEGER N, B INTEGER T, NM1, A, TEST T = N - 1 NM1 = T A = 0 10 CONTINUE C Using bit test function IF ( .NOT. BTEST (T,0) ) THEN C If BTEST (or similar function) is not available, replace by C (see remark in ISPRIM): C IF ( MOD ( N, 2 ) .EQ. 0 ) THEN T = T / 2 A = A + 1 GOTO 10 ENDIF C TEST = MODEXP ( B, T, N ) IF ( TEST .EQ. 1 .OR. TEST .EQ. NM1 ) THEN SPRP = .TRUE. RETURN ENDIF 20 CONTINUE TEST = MULMOD ( TEST, TEST, N ) IF ( TEST .EQ. NM1 ) THEN SPRP = .TRUE. RETURN ENDIF A = A - 1 IF ( A .GT. 0 ) GOTO 20 SPRP = .FALSE. END C...+....1....+....2....+....3....+....4....+....5....+....6....+....7.. INTEGER FUNCTION MODEXP ( A, B, M ) C A^B MOD M INTEGER A, B, M INTEGER N, BB, AA N = 1 BB = B AA = A 10 CONTINUE IF ( BB .LE. 0 ) GOTO 20 C Using bit test function IF ( BTEST ( BB, 0 ) ) N = MULMOD ( N, AA, M ) C If not available replace by (see remark in ISPRIM): C IF ( MOD ( BB, 2 ) .NE. 0 ) N = MULMOD ( N, AA, M ) BB = BB / 2 AA = MULMOD ( AA, AA, M ) GOTO 10 20 CONTINUE MODEXP = N END C...+....1....+....2....+....3....+....4....+....5....+....6....+....7.. INTEGER FUNCTION MULMOD ( A, B, M ) C A*B MOD M INTEGER A, B, M C C 64 bit Integer required (e.g. Digital Alpha, IRIX64, etc.) C To achieve a similar functionality on 32bit architectures see C e.g. Olivier Langlois' prime_cpp.html, function mulmod C C Local 64 bit variables INTEGER*8 AA, BB, P, MM AA = A BB = B C Product stays within 64 bit range for all (signed) 32bit integers P = AA * BB MM = M MULMOD = MOD ( P, MM ) END