; This is Scheme code. See scan for finding primes. ; The following primality test is from second edition of ; volume 2 of Knuth's "The Art of Computer Programming", ; page 379. n is the number to test and x is some random positive integer < n. (define (pt x n) (or (= n 2) (let* ((nm1 (- n 1))(pr (let z ((q nm1)(k 0))(if (odd? q)(cons q k) (z (quotient q 2)(+ k 1))))) (q (car pr))(k (cdr pr))) (let lp ((j 0)(y (mod-exp x q n))) (or (and (= j 0)(= y 1)) (= y nm1) (and (< (+ j 1) k)(lp (+ j 1)(modulo (* y y) n)))))))) ; 10^100-797, 10^200-189, 10^299-171, 10^300-69, 10^500+961 are prime. ; 2^1024 + 643 and 2^1024 - 105 too, (define (ex x)(write x)x) (define (rand31 seed) (lambda()(let* ((hi (quotient seed 127773)) (lo (- seed (* 127773 hi)))(test (- (* 16807 lo) (* 2836 hi)))) (set! seed (if (> test 0) test (+ test 2147483647))) seed))) (define (gbrand max seed)(let* ((rq (rand31 seed)) (nq (let v ((c max)(n 0))(if (< c (expt 2 31)) (cons (let w ((cx c)(nx 1))(if (< cx 2) nx (w (quotient cx 2)(+ nx 1))))n) (v (quotient c (expt 2 31))(+ n 1))))) (z (car nq))(n (cdr nq))(k (+ (* 31 n) z))) (write(list max k z n)) ; max, k, n, z, L and m are non-negative integers. ; 2^(k-1)<=max<2^k. k=31*n+z. 0= a P) (j (modulo a P) P) (if (even? a) (let((q (j (/ a 2) P)) (m (modulo P 8))) (if (or (= m 3)(= m 5))(- q) q)) (let((q (j (modulo P a) a))) (if(or (= (modulo P 4) 1) (= (modulo a 4) 1)) q (- q))))))) (define (mod-exp b p m)(cond ((= p 0) 1) ((even? p)(let ((x (mod-exp b (/ p 2) m)))(modulo (* x x)m))) (#t (modulo (* b (mod-exp b (- p 1) m)) m)))) ; The following is by Solovay & Strassen as presented in Knuth page 396. (define (p-test a P)(if(zero? a)(cdr 2))(and (odd? P) (= (gcd a P) 1) (zero? (modulo (- (j a P)(mod-exp a (/(- P 1) 2) P)) P)))) (quote "The function scan below returns the first prime in the") (quote "arithmetic sequence a + n*b") (define (scan a b)(let* ((g (gcd a b))(n (modulo b 30030)) (random (gbrand (+ a (if (positive? b) 0 (* 2000 b))) 228765)) (probe (+ 1 (random)))) (if (> g 1) (list "Always divisible by" g) (let more ((a1 a)(m (modulo a 30030))(cn 0)) (if (and (let all ((l (list 2 3 5 7 11 13)))(or (null? l) (and (positive? (remainder m (car l))) (all (cdr l))))) (let all ((l 20)(p probe)) (or (zero? l) (and (ww cn (pt p a1)) (all (- l 1)(+ 1 (random))))))) a1 (begin (display ",")(more (+ a1 b)(modulo (+ m n) 30030)(+ cn 1)))))))) (quote "The function scand below like scan but requires that both n & 2n+1 be prime") (define (scand a b)(let* ((ps (list 2 3 5 7 11 13))(pp (apply * ps)) (g (gcd a b))(g1 (gcd (+ a a 1)(+ b b)))(n (modulo b pp)) (random (gbrand (+ a (if (positive? b) 0 (* 2000 b))) 228765)) (probe (+ 1 (random)))) (cond ((> g 1) (list "Always divisible by" g)) ((> g1 1) (list "2(a+nb) + 1 always divisible by" g1)) (#t (let more ((a1 a)(m (modulo a pp))(cn 0)) (if (and (let all ((l ps))(or (null? l) (and (positive? (remainder m (car l))) (positive? (remainder (+ 1 m m) (car l))) (all (cdr l))))) (let all ((l 20)(p probe)) (or (zero? l) (and (ww cn (pt p a1)) (pt (+ p 15) (+ 1 a1 a1)) (all (- l 1)(+ 1 (random))))))) a1 (begin (display ",")(more (+ a1 b)(modulo (+ m n) pp)(+ cn 1))))))))) (define (pc wx) (let zz ((q (- wx 1))(s 0))(if (zero? q) s (zz (- q 1)(+ (if (p-test q wx) 1 0) s)))))