; combined n +/- 1 primality prover
; primes n -- list of primes not greater than n in ascending order
(define (primes n) ; assumes n is an integer greater than one
(let* ((len (quotient (- n 1) 2)) (bits (make-vector len #t)))
(let loop ((i 0) (p 3) (ps (list 2))) ; sieve of eratosthenes
(cond ((< n (* p p))
(do ((i i (+ i 1)) (p p (+ p 2))
(ps ps (if (vector-ref bits i) (cons p ps) ps)))
((= i len) (reverse ps))))
((vector-ref bits i)
(do ((j (+ (* 2 i i) (* 6 i) 3) (+ j p)))
((<= len j) (loop (+ i 1) (+ p 2) (cons p ps)))
(vector-set! bits j #f)))
(else (loop (+ i 1) (+ p 2) ps))))))
; prime? n -- #f if provably composite, else #t if probably prime
(define prime? ; strong pseudoprime to prime bases less than 100
(let ((ps (primes 100))) ; assumes n an integer greater than one
(lambda (n)
(define (expm b e m)
(let loop ((b b) (e e) (x 1))
(if (zero? e) x
(loop (modulo (* b b) m) (quotient e 2)
(if (odd? e) (modulo (* b x) m) x)))))
(define (spsp? n a) ; #t if n is a strong pseudoprime base a
(do ((d (- n 1) (/ d 2)) (s 0 (+ s 1)))
((odd? d) (if (= (expm a d n) 1) #t
(do ((r 0 (+ r 1)))
((or (= (expm a (* d (expt 2 r)) n) (- n 1))
(= r s))
(< r s)))))))
(if (< n 2) #f (if (member n ps) #t
(do ((ps ps (cdr ps)))
((or (null? ps) (not (spsp? n (car ps))))
(null? ps))))))))
; factors n -- list of prime factors of n in ascending order
(define (factors n) ; assumes n is an integer, may be negative
(if (<= -1 n 1) (list n) ; pollard's rho algorithm
(if (< n 0) (cons -1 (factors (- n)))
(let fact ((n n) (c 1) (fs (list)))
(define (f x) (modulo (+ (* x x) c) n))
(if (even? n) (fact (/ n 2) c (cons 2 fs))
(if (= n 1) fs
(if (prime? n) (sort < (cons n fs))
(let loop ((t 2) (h 2) (d 1))
(cond ((= d 1) (let ((t (f t)) (h (f (f h))))
(loop t h (gcd (- t h) n))))
((= d n) (fact n (+ c 1) fs)) ; cyclic
((prime? d)
(fact (/ n d) (+ c 1) (cons d fs)))
(else (fact n (+ c 1) fs)))))))))))
(define next-prime ; uses 2,3,5,7-wheel
(let ((gap (vector 1 10 9 8 7 6 5 4 3 2 1 2 1 4
3 2 1 2 1 4 3 2 1 6 5 4 3 2 1 2 1 6 5 4 3 2
1 4 3 2 1 2 1 4 3 2 1 6 5 4 3 2 1 6 5 4 3 2
1 2 1 6 5 4 3 2 1 4 3 2 1 2 1 6 5 4 3 2 1 4
3 2 1 6 5 4 3 2 1 8 7 6 5 4 3 2 1 4 3 2 1 2
1 4 3 2 1 2 1 4 3 2 1 8 7 6 5 4 3 2 1 6 5 4
3 2 1 4 3 2 1 6 5 4 3 2 1 2 1 4 3 2 1 6 5 4
3 2 1 2 1 6 5 4 3 2 1 6 5 4 3 2 1 4 3 2 1 2
1 4 3 2 1 6 5 4 3 2 1 2 1 6 5 4 3 2 1 4 3 2
1 2 1 4 3 2 1 2 1 10 9 8 7 6 5 4 3 2 1 2)))
(lambda (n)
(if (< n 2) 2 (if (< n 3) 3 (if (< n 5) 5
(let loop ((n (+ n (if (even? n) 1 2))))
(if (prime? n) n (loop (+ n
(vector-ref gap (modulo n 210))))))))))))
(define (expm b e m) ; modular exponentiation
(let loop ((b b) (e e) (x 1))
(if (zero? e) x
(loop (modulo (* b b) m) (quotient e 2)
(if (odd? e) (modulo (* b x) m) x)))))
(define (fermat-test? n ps)
(let loop ((ps ps) (a 1))
(if (null? ps) #t
(if (and (= (expm a (- n 1) n) 1)
(= (gcd (expt a (/ (- n 1) (car ps))) n) 1))
(loop (cdr ps) 1)
(loop ps (+ a 1))))))
(define (jacobi a m) ; jacobi symbol
(let loop1 ((a (modulo a m)) (m m) (t 1))
(if (zero? a) (if (= m 1) t 0)
(let ((z (if (member (modulo m 8) (list 3 5)) -1 1)))
(let loop2 ((a a) (t t))
(if (even? a) (loop2 (/ a 2) (* t z))
(loop1 (modulo m a) a
(if (and (= (modulo a 4) 3)
(= (modulo m 4) 3))
(- t) t))))))))
(define (selfridge n) ; initialize lucas sequence
(let loop ((d-abs 5) (sign 1))
(let ((d (* d-abs sign)))
(cond ((< 1 (gcd d n)) (values d 0 0))
((= (jacobi d n) -1) (values d 1 (/ (- 1 d) 4)))
(else (loop (+ d-abs 2) (- sign)))))))
(define (lucas p q m n) ; lucas sequences u[n] and v[n] and q^n (mod m)
(define (even e o) (if (even? n) e o))
(define (mod n) (if (zero? m) n (modulo n m)))
(let ((d (- (* p p) (* 4 q))))
(let loop ((un 1) (vn p) (qn q) (n (quotient n 2))
(u (even 0 1)) (v (even 2 p)) (k (even 1 q)))
(if (zero? n) (values u v k)
(let ((u2 (mod (* un vn))) (v2 (mod (- (* vn vn) (* 2 qn))))
(q2 (mod (* qn qn))) (n2 (quotient n 2)))
(if (even? n) (loop u2 v2 q2 n2 u v k)
(let* ((uu (+ (* u v2) (* u2 v)))
(vv (+ (* v v2) (* d u u2)))
(uu (if (and (positive? m) (odd? uu)) (+ uu m) uu))
(vv (if (and (positive? m) (odd? vv)) (+ vv m) vv))
(uu (mod (/ uu 2))) (vv (mod (/ vv 2))))
(loop u2 v2 q2 n2 uu vv (* k q2)))))))))
(define (lucas-test? n ps)
(call-with-values
(lambda () (selfridge n))
(lambda (d pp qq)
(let loop ((ps ps) (p pp) (q qq))
(if (null? ps) #t
(call-with-values
(lambda () (lucas p q n (+ n 1)))
(lambda (un+1 vn+1 qn+1)
(call-with-values
(lambda () (lucas p q n (/ (+ n 1) (car ps))))
(lambda (un+1/p vn+1/p qn+1/p)
(if (and (zero? un+1) (= (gcd un+1/p n) 1))
(loop (cdr ps) pp qq)
(loop ps (+ p 2) (+ p q 1))))))))))))
(define-syntax while
(syntax-rules ()
((while pred? body ...)
(do () ((not pred?)) body ...))))
(define (bls-prime? n . args)
(define (last-pair xs) (if (null? (cdr xs)) xs (last-pair (cdr xs))))
(define (cycle xs) (set-cdr! (last-pair xs) xs) xs)
(define wheel (list 2 4 2 4 6 2 6 4 2 4 6 6 2 6 4 2 6 4 6 8
4 2 4 2 4 8 6 4 6 2 4 6 2 6 6 4 2 4 6 2 6 4 2 4 2 10 2 10))
(define (rho n c b)
(define (f y) (modulo (+ (* y y) c) n))
(define (g p x y) (modulo (* p (abs (- x y))) n))
(let stage1 ((x 2) (y (+ 4 c)) (z (+ 4 c)) (j 1) (q 2) (p 1))
(if (= j b) #f ; timeout
(if (= x y) (rho n (+ c 1) (- b j)) ; cycle
(if (= j q) (let ((t (f y))) (stage1 y (f y) z (+ j 1) (* q 2) (g p y t)))
(if (positive? (modulo j 100)) (stage1 x (f y) z (+ j 1) q (g p x y))
(let ((d (gcd p n)))
(if (= d 1) (stage1 x (f y) y (+ j 1) q (g p x y))
(if (and (< 1 d n) (bls-prime? d)) d ; factor
(let stage2 ((k 1) (z (f z)))
(if (= k 100) (rho n (+ c 1) (- b j))
(let ((d (gcd (- z x) n)))
(if (= d 1) (stage2 (+ k 1) (f z))
(if (= d n) (rho n (+ c 1) (- b j))
(if (bls-prime? d) d ; factor
(rho n (+ c 1) (- b j)))))))))))))))))
(define (remove-twos n)
(let loop ((f 1) (r n))
(if (odd? r) (values f (list 2) r)
(loop (* f 2) (/ r 2)))))
(define (join x xs) (if (member x xs) xs (cons x xs)))
(define (enough? n p f1 f2)
(< n (* (max (+ (* p f1) 1) (- (* p f2) 1)) (+ (* p p f1 f2 1/2) 1))))
(if (not (prime? n)) #f ; sanity check
(let-values (((b1) (if (pair? args) (car args) 20000000))
((b2) (if (and (pair? args) (pair? (cdr args)))
(cadr args) 10000000))
((f1 f1s r1) (remove-twos (- n 1)))
((f2 f2s r2) (remove-twos (+ n 1))))
(let loop ((p 3) (ws (cons 2 (cons 2 (cons 4 (cycle wheel))))))
(cond ((and (< p b1) (not (enough? n p f1 f2)))
(case (modulo (+ n 1) p)
((0) (while (zero? (modulo r2 p))
(set! f2 (* f2 p)) (set! f2s (join p f2s)) (set! r2 (/ r2 p))))
((2) (while (zero? (modulo r1 p))
(set! f1 (* f1 p)) (set! f1s (join p f1s)) (set! r1 (/ r1 p)))))
(when (< r1 (* p p))
(set! f1 (- n 1)) (set! f1s (join r1 f1s)) (set! r1 1))
(when (< r2 (* p p))
(set! f2 (+ n 1)) (set! f2s (join r2 f2s)) (set! r2 1))
(loop (+ p (car ws)) (cdr ws)))
(else (set! p (min p b1))
(when verbose?
(display (list 'wheel n p f1 f1s r1 f2 f2s r2)) (newline))
(let loop ((more1? #t) (more2? #t) (more3? #t))
(cond ((and (not (enough? n p f1 f2)) (or more1? more2?))
(if more1?
(let ((f (rho r1 1 b2)))
(if (not f) (loop #f #t #t)
(begin (set! f1 (* f1 f)) (set! f1s (join f f1s))
(set! r1 (/ r1 f)) (loop #t #t#t))))
(let ((f (rho r2 1 b2)))
(if (not f) (loop #f #f #t)
(begin (set! f2 (* f2 f)) (set! f2s (join f f2s))
(set! r2 (/ r2 f)) (loop #f #t #t))))))
(more3?
(when verbose?
(display (list 'rho n p f1 f1s r1 f2 f2s r2)) (newline))
(loop #f #f #f))
((not (enough? n p f1 f2)) #f) ; failure to prove
((< r1 f1)
(when verbose? (display "fermat only") (newline))
(fermat-test? n f1s)) ; condition 1 only
((< r2 f2)
(when verbose? (display "lucas only") (newline))
(lucas-test? n f2s)) ; condition 2 only
(else
(when verbose? (display "fermat and lucas") (newline))
(and (fermat-test? n (cons r1 f1s)) ; 1 and 3
(lucas-test? n (cons r2 f2s)))))))))))) ; 2&4
(define rand #f)
(define randint #f)
(let ((two31 #x80000000) (a (make-vector 56 -1)) (fptr #f))
(define (mod-diff x y) (modulo (- x y) two31)) ; generic version
; (define (mod-diff x y) (logand (- x y) #x7FFFFFFF)) ; fast version
(define (flip-cycle)
(do ((ii 1 (+ ii 1)) (jj 32 (+ jj 1))) ((< 55 jj))
(vector-set! a ii (mod-diff (vector-ref a ii) (vector-ref a jj))))
(do ((ii 25 (+ ii 1)) (jj 1 (+ jj 1))) ((< 55 ii))
(vector-set! a ii (mod-diff (vector-ref a ii) (vector-ref a jj))))
(set! fptr 54) (vector-ref a 55))
(define (init-rand seed)
(let* ((seed (mod-diff seed 0)) (prev seed) (next 1))
(vector-set! a 55 prev)
(do ((i 21 (modulo (+ i 21) 55))) ((zero? i))
(vector-set! a i next) (set! next (mod-diff prev next))
(set! seed (+ (quotient seed 2) (if (odd? seed) #x40000000 0)))
(set! next (mod-diff next seed)) (set! prev (vector-ref a i)))
(flip-cycle) (flip-cycle) (flip-cycle) (flip-cycle) (flip-cycle)))
(define (next-rand)
(if (negative? (vector-ref a fptr)) (flip-cycle)
(let ((next (vector-ref a fptr))) (set! fptr (- fptr 1)) next)))
(define (unif-rand m)
(let ((t (- two31 (modulo two31 m))))
(let loop ((r (next-rand)))
(if (<= t r) (loop (next-rand)) (modulo r m)))))
(init-rand 19380110) ; happy birthday donald e knuth
(set! rand (lambda seed
(cond ((null? seed) (/ (next-rand) two31))
((eq? (car seed) 'get) (cons fptr (vector->list a)))
((eq? (car seed) 'set) (set! fptr (caadr seed))
(set! a (list->vector (cdadr seed))))
(else (/ (init-rand (modulo (numerator
(inexact->exact (car seed))) two31)) two31)))))
(set! randint (lambda args
(cond ((null? (cdr args))
(if (< (car args) two31) (unif-rand (car args))
(floor (* (next-rand) (car args)))))
((< (car args) (cadr args))
(let ((span (- (cadr args) (car args))))
(+ (car args)
(if (< span two31) (unif-rand span)
(floor (* (next-rand) span))))))
(else (let ((span (- (car args) (cadr args))))
(- (car args)
(if (< span two31) (unif-rand span)
(floor (* (next-rand) span))))))))))
(define (rand-prime n . args)
(let ((b (if (pair? args) (car args) 10)))
(let loop ((r (randint 1 b)) (n n))
(if (= n 1) (next-prime r)
(loop (+ (* r b) (randint b)) (- n 1))))))
(define verbose? #f)
(let ((n (rand-prime 42)))
(display "Proving primality of ")
(display n) (display ": ")
(display (bls-prime? n)) (newline))