; baillie-wagstaff pseudoprimality test
(define (isqrt n)
(if (not (and (positive? n) (integer? n)))
(error 'isqrt "must be positive integer")
(let loop ((x n))
(let ((y (quotient (+ x (quotient n x)) 2)))
(if (< y x) (loop y) x)))))
(define (square? n)
(let ((n2 (isqrt n)))
(= (* n2 n2) n)))
(define (expm b e m)
(define (m* x y) (modulo (* x y) m))
(cond ((zero? e) 1)
((even? e) (expm (m* b b) (/ e 2) m))
(else (m* b (expm (m* b b) (/ (- e 1) 2) m)))))
(define (jacobi a m)
(if (not (integer? a)) (error 'jacobi "must be integer")
(if (not (and (integer? m) (positive? m) (odd? m)))
(error 'jacobi "modulus must be odd positive integer")
(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 (primes n)
(let ((sieve (make-vector n #t)))
(let loop ((p 2) (ps (list)))
(cond ((= n p) (reverse ps))
((vector-ref sieve p)
(do ((i p (+ i p))) ((<= n i))
(vector-set! sieve i #f))
(loop (+ p 1) (cons p ps)))
(else (loop (+ p 1) ps))))))
(define (strong-pseudoprime? n a)
(let loop ((r 0) (s (- n 1)))
(if (even? s) (loop (+ r 1) (/ s 2))
(if (= (expm a s n) 1) #t
(let loop ((r r) (s s))
(cond ((zero? r) #f)
((= (expm a s n) (- n 1)) #t)
(else (loop (- r 1) (* s 2)))))))))
(define (selfridge n)
(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 (chain n u v u2 v2 d q m)
(let loop ((u u) (v v) (u2 u2) (v2 v2) (q q) (qkd q) (m m))
(if (zero? m) (values u v qkd)
(let* ((u2 (modulo (* u2 v2) n))
(v2 (modulo (- (* v2 v2) (* 2 q)) n))
(q (modulo (* q q) n)))
(if (odd? m)
(let* ((t1 (* u2 v)) (t2 (* u v2))
(t3 (* v2 v)) (t4 (* u2 u d))
(u (+ t1 t2)) (v (+ t3 t4))
(qkd (modulo (* qkd q) n)))
(loop (modulo (quotient (if (odd? u) (+ u n) u) 2) n)
(modulo (quotient (if (odd? v) (+ v n) v) 2) n)
u2 v2 q qkd (quotient m 2)))
(loop u v u2 v2 q qkd (quotient m 2)))))))
(define (standard-lucas-pseudoprime? n)
; assumes odd positive integer not a square
(call-with-values
(lambda () (selfridge n))
(lambda (d p q)
(if (zero? p) (= n d)
(call-with-values
(lambda () (chain n 0 2 1 p d q (/ (+ n 1) 2)))
(lambda (u v qkd) (zero? u)))))))
(define (powers-of-two n)
(let loop ((s 0) (n n))
(if (odd? n) (values s n)
(loop (+ s 1) (/ n 2)))))
(define (strong-lucas-pseudoprime? n)
; assumes odd positive integer not a square
(call-with-values
(lambda () (selfridge n))
(lambda (d p q)
(if (zero? p) (= n d)
(call-with-values
(lambda () (powers-of-two (+ n 1)))
(lambda (s t)
(call-with-values
(lambda () (chain n 1 p 1 p d q (quotient t 2)))
(lambda (u v qkd)
(if (or (zero? u) (zero? v)) #t
(let loop ((r 1) (v v) (qkd qkd))
(if (= r s) #f
(let* ((v (modulo (- (* v v) (* 2 qkd)) n))
(qkd (modulo (* qkd qkd) n)))
(if (zero? v) #t (loop (+ r 1) v qkd))))))))))))))
(define prime?
(let ((ps (primes 100)))
(lambda (n)
(if (not (integer? n)) (error 'prime? "must be integer"))
(if (or (< n 2) (square? n)) #f
(let loop ((ps ps))
(if (pair? ps)
(if (zero? (modulo n (car ps))) (= n (car ps)) (loop (cdr ps)))
(and (strong-pseudoprime? n 2)
(strong-pseudoprime? n 3)
(strong-lucas-pseudoprime? n))))))))
(do ((ns (list 79 83 89 111 323 5777 3825123056546413051 (- (expt 2 89) 1)) (cdr ns))) ((null? ns))
(let ((n (car ns)))
(display n) (display " ")
(display (strong-pseudoprime? n 2)) (display " ")
(display (strong-pseudoprime? n 3)) (display " ")
(display (standard-lucas-pseudoprime? n)) (display " ")
(display (strong-lucas-pseudoprime? n)) (display " ")
(display (prime? n)) (newline)))
(display ";;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;") (newline)
(define (lucas p q m n)
(define (mod n) (if (zero? m) n (modulo n m)))
(let bits ((n n) (bs (list)))
(if (positive? n) (bits (quotient n 2) (cons (modulo n 2) bs))
(let loop ((un 0) (uu 1) (vn 2) (vv p) (qn 1) (bs bs))
(cond ((null? bs) (values un vn qn))
((odd? (car bs)) ; 1-bit
(loop (mod (- (* uu vn) qn)) ; u(2n+1)
(mod (* uu vv)) ; u(2n+2)
(mod (- (* vv vn) (* p qn))) ; v(2n+1)
(mod (- (* vv vv) (* 2 q qn))) ; v(2n+2)
(mod (* qn qn q)) (cdr bs)))
(else ; 0-bit
(loop (mod (* un vn)) ; u(2n)
(mod (- (* uu vn) qn)) ; u(2n+1)
(mod (- (* vn vn) (* 2 qn))) ; v(2n)
(mod (- (* vv vn) (* p qn))) ; v(2n+1)
(mod (* qn qn)) (cdr bs))))))))
(define (standard-lucas-pseudoprime? n)
; assumes odd positive integer not a square
(call-with-values
(lambda () (selfridge n))
(lambda (d p q)
(if (zero? p) (= n d)
(call-with-values
; (lambda () (chain n 0 2 1 p d q (/ (+ n 1) 2)))
(lambda () (lucas p q n (/ (+ n 1) 2)))
(lambda (u v qkd) (zero? u)))))))
(define (powers-of-two n)
(let loop ((s 0) (n n))
(if (odd? n) (values s n)
(loop (+ s 1) (/ n 2)))))
(define (strong-lucas-pseudoprime? n)
; assumes odd positive integer not a square
(call-with-values
(lambda () (selfridge n))
(lambda (d p q)
(if (zero? p) (= n d)
(call-with-values
(lambda () (powers-of-two (+ n 1)))
(lambda (s t)
(call-with-values
; (lambda () (chain n 1 p 1 p d q (quotient t 2)))
(lambda () (lucas p q n (quotient t 2)))
(lambda (u v qkd)
(if (or (zero? u) (zero? v)) #t
(let loop ((r 1) (v v) (qkd qkd))
(if (= r s) #f
(let* ((v (modulo (- (* v v) (* 2 qkd)) n))
(qkd (modulo (* qkd qkd) n)))
(if (zero? v) #t (loop (+ r 1) v qkd))))))))))))))
(define prime?
(let ((ps (primes 100)))
(lambda (n)
(if (not (integer? n)) (error 'prime? "must be integer"))
(if (or (< n 2) (square? n)) #f
(let loop ((ps ps))
(if (pair? ps)
(if (zero? (modulo n (car ps))) (= n (car ps)) (loop (cdr ps)))
(and (strong-pseudoprime? n 2)
(strong-pseudoprime? n 3)
(strong-lucas-pseudoprime? n))))))))
(do ((ns (list 79 83 89 111 323 5777 3825123056546413051 (- (expt 2 89) 1)) (cdr ns))) ((null? ns))
(let ((n (car ns)))
(display n) (display " ")
(display (strong-pseudoprime? n 2)) (display " ")
(display (strong-pseudoprime? n 3)) (display " ")
(display (standard-lucas-pseudoprime? n)) (display " ")
(display (strong-lucas-pseudoprime? n)) (display " ")
(display (prime? n)) (newline)))