```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 ``` ```; sophie germain primes (define (inverse x m) (let loop ((x x) (a 0) (b m) (u 1)) (if (positive? x) (loop (modulo b x) u x (- a (* (quotient b x) u))) (if (= b 1) (modulo a m) (error 'inverse "must be coprime"))))) (define (drop n xs) (let loop ((n n) (xs xs)) (if (or (zero? n) (null? xs)) xs (loop (- n 1) (cdr xs))))) (define sort #f) (define merge #f) (let () (define dosort (lambda (pred? ls n) (if (= n 1) (list (car ls)) (let ((i (quotient n 2))) (domerge pred? (dosort pred? ls i) (dosort pred? (list-tail ls i) (- n i))))))) (define domerge (lambda (pred? l1 l2) (cond ((null? l1) l2) ((null? l2) l1) ((pred? (car l2) (car l1)) (cons (car l2) (domerge pred? l1 (cdr l2)))) (else (cons (car l1) (domerge pred? (cdr l1) l2)))))) (set! sort (lambda (pred? l) (if (null? l) l (dosort pred? l (length l))))) (set! merge (lambda (pred? l1 l2) (domerge pred? l1 l2)))) (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 (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 (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 (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 (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)))))))) (define (germain b f m) (define (p c) (- (* c 3003 (expt 10 b)) 1)) (define (q c) (- (* c 6006 (expt 10 b)) 1)) (define (r c) (- (* c 15015 (expt 10 (- b 1))) 1)) (define (p-start p) (inverse (* 3003 (expt 10 b)) p)) (define (q-start p) (inverse (* 6006 (expt 10 b)) p)) (define (r-start p) (inverse (* 15015 (expt 10 (- b 1))) p)) (let ((sieve (make-vector (+ m 1) #t)) ; ignore sieve[0] (ps (drop 6 (primes f)))) ; sieving primes ; sieve on p (do ((ps ps (cdr ps))) ((null? ps)) (do ((i (p-start (car ps)) (+ i (car ps)))) ((< m i)) (vector-set! sieve i #f))) ; sieve on q (do ((ps ps (cdr ps))) ((null? ps)) (do ((i (q-start (car ps)) (+ i (car ps)))) ((< m i)) (vector-set! sieve i #f))) ; sieve on r (do ((ps ps (cdr ps))) ((null? ps)) (do ((i (r-start (car ps)) (+ i (car ps)))) ((< m i)) (vector-set! sieve i #f))) ; collect and test primes (let ((gs (list))) ; sophie-germain primes (do ((c 1 (+ c 1))) ((< m c) (sort < gs)) (when (and (vector-ref sieve c) (prime? (p c))) (when (prime? (q c)) (set! gs (cons (p c) gs))) (when (prime? (r c)) (set! gs (cons (r c) gs)))))))) (for-each (lambda (p) (display p) (newline)) (germain 25 100000 100000)) ```
 ```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 208 209 210 ``` ```1831829999999999999999999999999 48873824999999999999999999999999 82912829999999999999999999999999 91351259999999999999999999999999 97747649999999999999999999999999 97747649999999999999999999999999 119339219999999999999999999999999 119339219999999999999999999999999 120585464999999999999999999999999 132417284999999999999999999999999 154609454999999999999999999999999 155705549999999999999999999999999 163258094999999999999999999999999 166291124999999999999999999999999 168468299999999999999999999999999 168468299999999999999999999999999 169729559999999999999999999999999 177747569999999999999999999999999 184849664999999999999999999999999 189309119999999999999999999999999 195495299999999999999999999999999 207086879999999999999999999999999 208888679999999999999999999999999 209309099999999999999999999999999 209309099999999999999999999999999 210194984999999999999999999999999 217176959999999999999999999999999 231366134999999999999999999999999 236411174999999999999999999999999 240164924999999999999999999999999 251411159999999999999999999999999 285795509999999999999999999999999 295254959999999999999999999999999 306561254999999999999999999999999 323963639999999999999999999999999 326035709999999999999999999999999 329459129999999999999999999999999 347266919999999999999999999999999 367266899999999999999999999999999 380029649999999999999999999999999 384504119999999999999999999999999 396290894999999999999999999999999 420389969999999999999999999999999 458933474999999999999999999999999 460525064999999999999999999999999 479864384999999999999999999999999 480329849999999999999999999999999 487897409999999999999999999999999 507957449999999999999999999999999 507957449999999999999999999999999 511696184999999999999999999999999 513873359999999999999999999999999 514023509999999999999999999999999 531320789999999999999999999999999 549488939999999999999999999999999 561275714999999999999999999999999 571696124999999999999999999999999 579894314999999999999999999999999 581711129999999999999999999999999 591065474999999999999999999999999 603317714999999999999999999999999 635629994999999999999999999999999 647581934999999999999999999999999 658918259999999999999999999999999 659233574999999999999999999999999 676065389999999999999999999999999 686876189999999999999999999999999 694533839999999999999999999999999 694533839999999999999999999999999 718677959999999999999999999999999 732446714999999999999999999999999 734488754999999999999999999999999 760059299999999999999999999999999 760059299999999999999999999999999 766440674999999999999999999999999 779188409999999999999999999999999 787581794999999999999999999999999 807611804999999999999999999999999 832942109999999999999999999999999 837431594999999999999999999999999 850269419999999999999999999999999 852221369999999999999999999999999 853002149999999999999999999999999 874068194999999999999999999999999 877056179999999999999999999999999 877491614999999999999999999999999 886650764999999999999999999999999 887176289999999999999999999999999 892671779999999999999999999999999 899443544999999999999999999999999 902671769999999999999999999999999 902671769999999999999999999999999 907041134999999999999999999999999 908797889999999999999999999999999 911260349999999999999999999999999 921050129999999999999999999999999 921560639999999999999999999999999 931560629999999999999999999999999 933347414999999999999999999999999 947626679999999999999999999999999 951935984999999999999999999999999 959728769999999999999999999999999 959728769999999999999999999999999 964713749999999999999999999999999 971095124999999999999999999999999 990329339999999999999999999999999 999668669999999999999999999999999 999668669999999999999999999999999 1025854829999999999999999999999999 1030914884999999999999999999999999 1039308269999999999999999999999999 1047971924999999999999999999999999 1049143094999999999999999999999999 1052866814999999999999999999999999 1062671609999999999999999999999999 1077146069999999999999999999999999 1078602524999999999999999999999999 1083947864999999999999999999999999 1084668584999999999999999999999999 1089398309999999999999999999999999 1104533429999999999999999999999999 1126350224999999999999999999999999 1132040909999999999999999999999999 1132040909999999999999999999999999 1138302164999999999999999999999999 1159788629999999999999999999999999 1177686509999999999999999999999999 1180884704999999999999999999999999 1181064884999999999999999999999999 1188977789999999999999999999999999 1189908719999999999999999999999999 1204142939999999999999999999999999 1230839609999999999999999999999999 1230839609999999999999999999999999 1239368129999999999999999999999999 1246154909999999999999999999999999 1277401124999999999999999999999999 1285389104999999999999999999999999 1296320024999999999999999999999999 1318467149999999999999999999999999 1323542219999999999999999999999999 1325719394999999999999999999999999 1332521189999999999999999999999999 1340058719999999999999999999999999 1340058719999999999999999999999999 1342386044999999999999999999999999 1342821479999999999999999999999999 1344353009999999999999999999999999 1346560214999999999999999999999999 1347025679999999999999999999999999 1347671324999999999999999999999999 1352130779999999999999999999999999 1395328934999999999999999999999999 1404908504999999999999999999999999 1425058634999999999999999999999999 1427566139999999999999999999999999 1438151714999999999999999999999999 1439848409999999999999999999999999 1447340894999999999999999999999999 1447400954999999999999999999999999 1470358889999999999999999999999999 1470358889999999999999999999999999 1473121649999999999999999999999999 1485914429999999999999999999999999 1485914429999999999999999999999999 1491154664999999999999999999999999 1492671179999999999999999999999999 1492671179999999999999999999999999 1493151659999999999999999999999999 1493181689999999999999999999999999 1502370869999999999999999999999999 1538346809999999999999999999999999 1558376819999999999999999999999999 1615223609999999999999999999999999 1616454839999999999999999999999999 1653752099999999999999999999999999 1677265589999999999999999999999999 1706004299999999999999999999999999 1719728009999999999999999999999999 1798887089999999999999999999999999 1814082269999999999999999999999999 1866694829999999999999999999999999 2013421409999999999999999999999999 2074172099999999999999999999999999 2090478389999999999999999999999999 2096003909999999999999999999999999 2105733629999999999999999999999999 2126874749999999999999999999999999 2264081819999999999999999999999999 2332249919999999999999999999999999 2355373019999999999999999999999999 2361769409999999999999999999999999 2408285879999999999999999999999999 2438676239999999999999999999999999 2443150709999999999999999999999999 2459456999999999999999999999999999 2467715249999999999999999999999999 2469306839999999999999999999999999 2554802249999999999999999999999999 2595582989999999999999999999999999 2597354759999999999999999999999999 2629276649999999999999999999999999 2791768979999999999999999999999999 2809817009999999999999999999999999 2838315479999999999999999999999999 2866273409999999999999999999999999 2894681789999999999999999999999999 2900807909999999999999999999999999 2916273359999999999999999999999999 2940717779999999999999999999999999 ```