```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 ``` ```; aks primality prover (define (split n xs) (let loop ((n n) (xs xs) (zs '())) (if (or (zero? n) (null? xs)) (values (reverse zs) xs) (loop (- n 1) (cdr xs) (cons (car xs) zs))))) (define (make-list n x) (let loop ((n n) (xs '())) (if (zero? n) xs (loop (- n 1) (cons x xs))))) (define (square x) (* x x)) (define (log2 n) (/ (log n) (log 2))) (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 (ilog b n) (let loop1 ((lo 0) (b^lo 1) (hi 1) (b^hi b)) (if (< b^hi n) (loop1 hi b^hi (* hi 2) (* b^hi b^hi)) (let loop2 ((lo lo) (b^lo b^lo) (hi hi) (b^hi b^hi)) (if (<= (- hi lo) 1) (if (= b^hi n) hi lo) (let* ((mid (quotient (+ lo hi) 2)) (b^mid (* b^lo (expt b (- mid lo))))) (cond ((< n b^mid) (loop2 lo b^lo mid b^mid)) ((< b^mid n) (loop2 mid b^mid hi b^hi)) (else mid)))))))) (define (td-prime? n) (if (even? n) (= n 2) (let loop ((d 3)) (cond ((< n (* d d)) #t) ((zero? (modulo n d)) #f) (else (loop (+ d 2))))))) (define (uniq-factors n) (define (uniq-cons x xs) (if (null? xs) (list x) (if (= x (car xs)) xs (cons x xs)))) (let twos ((n n) (fs '())) (if (even? n) (twos (/ n 2) (uniq-cons 2 fs)) (if (= n 1) fs (let odds ((n n) (d 3) (fs fs)) (cond ((< n (* d d)) (reverse (uniq-cons n fs))) ((zero? (modulo n d)) (odds (/ n d) d (uniq-cons d fs))) (else (odds n (+ d 2) fs)))))))) (define (prime-power? n) (if (even? n) (if (= (expt 2 (ilog 2 n)) n) 2 #f) (let loop ((a 2) (n n)) (let* ((b (expm a n n)) (p (gcd (- b a) n))) (cond ((or (= p 1) (= a b)) #f) ((prime? p) (if (= (expt p (ilog p n)) n) p #f)) (else (loop (+ a 1) p))))))) ; The order of n modulo r, ; written as ord_r(n), ; is the smallest number k ; such that n^k = 1 (mod r). (define (ord r n) (do ((k 2 (+ k 1))) ((= (expm n k r) 1) k))) ; Find the smallest r such ; that ord_r(n) > 4(log2 n)^2. (define (compute-r n) (let ((target (* 4 (square (log2 n))))) (do ((r 3 (+ r 1))) ((and (td-prime? r) (not (= r n)) (< target (ord r n))) r)))) (define (phi n) (let loop ((fs (uniq-factors n)) (t n)) (if (null? fs) t (loop (cdr fs) (* t (- 1 (/ (car fs)))))))) (define (poly-mult-mod xs ys r n) (define (times x) (lambda (y) (* x y))) (define (plus xs ys) (let loop ((xs xs) (ys ys) (zs (list))) (cond ((null? xs) (reverse (append (reverse ys) zs))) ((null? ys) (reverse (append (reverse xs) zs))) (else (loop (cdr xs) (cdr ys) (cons (+ (car xs) (car ys)) zs)))))) (define (mod-poly xs) (let-values (((hs ts) (split r (reverse xs)))) (reverse (plus hs ts)))) (define (mod-n x) (modulo x n)) (let ((xs (reverse xs))) (let loop ((xs (cdr xs)) (zs (map (times (car xs)) ys))) (if (null? xs) (map mod-n (mod-poly zs)) (loop (cdr xs) (plus (cons 0 zs) (map (times (car xs)) ys))))))) (define (poly-power-mod bs e r n) (let loop ((bs bs) (e e) (rs (list 1))) (if (zero? e) rs (loop (poly-mult-mod bs bs r n) (quotient e 2) (if (even? e) rs (poly-mult-mod rs bs r n)))))) (define (binomial-test? a r n) (not (equal? (poly-power-mod (list 1 a) n r n) (append (list 1) (make-list (- r 1) 0) (list a))))) (define (aks-prime? n) (if (prime-power? n) #f (let* ((r (compute-r n)) (phi-r (phi r)) (sqrt-phi-r (sqrt phi-r)) (log2-n (log2 n)) (sqrt-phi-r-log2-n (* sqrt-phi-r log2-n))) (let loop ((i 2)) (if (< i sqrt-phi-r-log2-n) (if (and (not (= n i)) (zero? (modulo n i))) #f (loop (+ i 1))) (let loop ((a 1)) (if (< a sqrt-phi-r-log2-n) (if (binomial-test? a r n) (loop (+ a 1)) #f) #t))))))) (display (compute-r 89)) (newline) (display (aks-prime? 89)) (newline) ```
 ```1 2 ``` ```191 #t ```