; 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<z<=31. 
; I claim this code is efficient if 
; (quotient .. (expt 2 ..)) and (modulo .. (expt 2 ..)) are efficient.
  (if (= n 0)(lambda () (let r ()(let ((q (modulo (rq) (expt 2 z))))
     (if (< q max) q (r)))))
   (let ((L (quotient max (expt 2 (- k 31)))))(write L)
; max = L*2^(k-31) + m; 0<=m<2^(k-31). L is the high 31 bits of max.
    (lambda()(+ (* (expt 2 z)(let g ((nx n))(if (= nx 1)
      (let r ()(let ((q (rq)))(if (< q L) q (r))))
      (+ (* (expt 2 31)(g (- nx 1)))(rq)))))(modulo (rq)(expt 2 z))))))))


(define ww (lambda(n x)(display (list n x)) x))
; The following is the "Jacobi Symbol". (a|P)
(define (j a P) (if (< a 2) (if (= a 1) 1
      (if (= a 0) 0 (j (modulo a P) P)))
   (if (>= 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)))))