Project:
 ```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 ``` ```; cornacchia's algorithm ; 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) (define (times p q) (modulo (* p q) m)) (let loop ((b b) (e e) (x 1)) (if (zero? e) x (loop (times b b) (quotient e 2) (if (odd? e) (times b x) 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 (member n ps) #t (do ((ps ps (cdr ps))) ((or (null? ps) (not (spsp? n (car ps)))) (null? ps))))))) ; expm b e m -- modular exponentiation: b^e (mod m) (define (expm b e m) (define (times p q) (modulo (* p q) m)) (let loop ((b b) (e e) (x 1)) (if (zero? e) x (loop (times b b) (quotient e 2) (if (odd? e) (times b x) x))))) ; jacobi a m -- jacobi symbol (define (jacobi a m) (let loop ((a (modulo a m)) (m m) (t 1)) (if (= a 0) (if (= m 1) t 0) (let twos ((a a) (t t)) (if (even? a) (twos (quotient a 2) (if (member (modulo m 8) '(3 5)) (- t) t)) (loop (modulo m a) a (if (= 3 (modulo a 4) (modulo m 4)) (- t) t))))))) ; legendre a n -- legendre symbol (define (legendre a n) (if (and (< 2 n) (prime? n)) (jacobi a n) (error 'legendre "modulus must be odd prime"))) ; mod-sqrt a p -- modular square root: x such that x^2 == a (mod p) (define (mod-sqrt a p) ; tonelli-shanks algorithm (define (both n) (list n (- p n))) (when (not (and (odd? p) (prime? p))) (error 'mod-sqrt "modulus must be an odd prime")) (when (not (= (legendre a p) 1)) (error 'mod-sqrt "must be a quadratic residual")) (let ((a (modulo a p))) (case (modulo p 8) ((3 7) (both (expm a (/ (+ p 1) 4) p))) ((5) (let* ((x (expm a (/ (+ p 3) 8) p)) (c (expm x 2 p))) (if (= a c) (both x) (both (modulo (* x (expm 2 (/ (- p 1) 4) p)) p))))) (else (let* ((d (let loop ((d 2)) (if (= (jacobi d p) -1) d (loop (+ d 1))))) (s (let loop ((p (- p 1)) (s 0)) (if (odd? p) s (loop (quotient p 2) (+ s 1))))) (t (quotient (- p 1) (expt 2 s))) (big-a (expm a t p)) (big-d (expm d t p)) (m (let loop ((i 0) (m 0)) (cond ((= i s) m) ((= (- p 1) (expm (* big-a (expm big-d m p)) (expt 2 (- s 1 i)) p)) (loop (+ i 1) (+ m (expt 2 i)))) (else (loop (+ i 1) m)))))) (both (modulo (* (expm a (/ (+ t 1) 2) p) (expm big-d (/ m 2) p)) p))))))) ; isqrt n -- largest integer x such that x*x <= n (define (isqrt n) (let loop ((x n) (y (quotient (+ n 1) 2))) (if (<= 0 (- y x) 1) x (loop y (quotient (+ y (quotient n y)) 2))))) ; divides? d n -- #t if n/d is integral, else #f (define (divides? d n) (zero? (modulo n d))) ; square? n -- sqrt(n) if n is a perfect square, else #f (define square? (let ((q11 (make-vector 11 #f)) (q63 (make-vector 63 #f)) (q64 (make-vector 64 #f)) (q65 (make-vector 65 #f))) (do ((k 0 (+ k 1))) ((< 5 k)) (vector-set! q11 (modulo (* k k) 11) #t)) (do ((k 0 (+ k 1))) ((< 31 k)) (vector-set! q63 (modulo (* k k) 63) #t)) (do ((k 0 (+ k 1))) ((< 31 k)) (vector-set! q64 (modulo (* k k) 64) #t)) (do ((k 0 (+ k 1))) ((< 32 k)) (vector-set! q65 (modulo (* k k) 65) #t)) (lambda (n) (if (not (integer? n)) (error 'square? "must be integer") (if (< n 1) #f (if (not (vector-ref q64 (modulo n 64))) #f (let ((r (modulo n 45045))) (if (not (vector-ref q63 (modulo r 63))) #f (if (not (vector-ref q65 (modulo r 65))) #f (if (not (vector-ref q11 (modulo r 11))) #f (let ((q (isqrt n))) (if (= (* q q) n) q #f)))))))))))) (define-syntax fold-of (syntax-rules (range in is) ((_ "z" f b e) (set! b (f b e))) ((_ "z" f b e (v range fst pst stp) c ...) (let* ((x fst) (p pst) (s stp) (le? (if (positive? s) <= >=))) (do ((v x (+ v s))) ((le? p v) b) (fold-of "z" f b e c ...)))) ((_ "z" f b e (v range fst pst) c ...) (let* ((x fst) (p pst) (s (if (< x p) 1 -1))) (fold-of "z" f b e (v range x p s) c ...))) ((_ "z" f b e (v range pst) c ...) (fold-of "z" f b e (v range 0 pst) c ...)) ((_ "z" f b e (x in xs) c ...) (do ((t xs (cdr t))) ((null? t) b) (let ((x (car t))) (fold-of "z" f b e c ...)))) ((_ "z" f b e (x is y) c ...) (let ((x y)) (fold-of "z" f b e c ...))) ((_ "z" f b e p? c ...) (if p? (fold-of "z" f b e c ...))) ((_ f i e c ...) (let ((b i)) (fold-of "z" f b e c ...))))) (define-syntax list-of (syntax-rules () ((_ arg ...) (reverse (fold-of (lambda (d a) (cons a d)) '() arg ...))))) (define (cornacchia d p) (if (negative? (jacobi (- d) p)) #f (let loop ((a p) (b (apply max (mod-sqrt (- d) p)))) (if (< p (* b b)) (loop b (modulo a b)) (let* ((n (- p (* b b))) (c (/ n d))) (cond ((not (divides? d n)) #f) ((not (square? c)) #f) (else (list b (isqrt c))))))))) (display (cornacchia 4 1733)) (newline) (display (cornacchia 58 3031201)) (newline) (display (list-of (cons p (cornacchia 1 p)) (p in (primes 1000)) (= (modulo p 4) 1))) (newline) ```
 ```1 2 3 ``` ```(17 19) (127 228) ((5 2 1) (13 3 2) (17 4 1) (29 5 2) (37 6 1) (41 5 4) (53 7 2) (61 6 5) (73 8 3) (89 8 5) (97 9 4) (101 10 1) (109 10 3) (113 8 7) (137 11 4) (149 10 7) (157 11 6) (173 13 2) (181 10 9) (193 12 7) (197 14 1) (229 15 2) (233 13 8) (241 15 4) (257 16 1) (269 13 10) (277 14 9) (281 16 5) (293 17 2) (313 13 12) (317 14 11) (337 16 9) (349 18 5) (353 17 8) (373 18 7) (389 17 10) (397 19 6) (401 20 1) (409 20 3) (421 15 14) (433 17 12) (449 20 7) (457 21 4) (461 19 10) (509 22 5) (521 20 11) (541 21 10) (557 19 14) (569 20 13) (577 24 1) (593 23 8) (601 24 5) (613 18 17) (617 19 16) (641 25 4) (653 22 13) (661 25 6) (673 23 12) (677 26 1) (701 26 5) (709 22 15) (733 27 2) (757 26 9) (761 20 19) (769 25 12) (773 22 17) (797 26 11) (809 28 5) (821 25 14) (829 27 10) (853 23 18) (857 29 4) (877 29 6) (881 25 16) (929 23 20) (937 24 19) (941 29 10) (953 28 13) (977 31 4) (997 31 6)) ```