```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 ``` ```; 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))) ```
 ```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 ``` ```79 #t #t #t #t #t 83 #t #t #t #t #t 89 #t #t #t #t #t 111 #f #f #f #f #f 323 #f #f #t #f #f 5777 #f #f #t #t #f 3825123056546413051 #t #t #f #f #f 618970019642690137449562111 #t #t #t #t #t ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 79 #t #t #f #f #t 83 #t #t #f #f #t 89 #t #t #t #f #t 111 #f #f #f #f #f 323 #f #f #t #f #f 5777 #f #f #t #f #f 3825123056546413051 #t #t #f #f #f 618970019642690137449562111 #t #t #t #t #t ```