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
|
;; "primes.scm", test and generate prime numbers.
; Written by Michael H Coffin (mhc@edsdrd.eds.com)
;
; This code is in the public domain.
;Date: Thu, 23 Feb 1995 07:47:49 +0500
;From: mhc@edsdrd.eds.com (Michael H Coffin)
;;
;; Test numbers for primality using Rabin-Miller Monte-Carlo
;; primality test.
;;
;; Public functions:
;;
;; (primes start count . iter)
;;
;; (probably-prime? p . iter)
;;
;;
;; Please contact the author if you have problems or suggestions:
;;
;; Mike Coffin
;; 1196 Whispering Knoll
;; Rochester Hills, Mi. 48306
;;
;; mhc@edsdrd.eds.com
;;
(require 'random)
;; The default number of times to perform the Rabin-Miller test. The
;; probability of a composite number passing the Rabin-Miller test for
;; primality with this many random numbers is at most
;; 1/(4^primes:iterations). The default yields about 1e-9.
;;
(define primes:iter 15)
;; Is n probably prime?
;;
(define (primes:probably-prime? n . iter)
(let ((iter (if (null? iter) primes:iter (car iter))))
(primes:prob-pr? n iter)))
;; Return a list of the first `number' odd probable primes less
;; than `start'.
(define (primes:primes< start number . iter)
(let ((iter (if (null? iter) primes:iter (car iter))))
(do ((candidate (if (odd? start) start (- start 1))
(- candidate 2))
(count 0)
(result '())
)
((or (< candidate 3) (>= count number)) result)
(if (primes:prob-pr? candidate iter)
(begin
(set! count (1+ count))
(set! result (cons candidate result)))
))))
(define (primes:primes> start number . iter)
(let ((iter (if (null? iter) primes:iter (car iter))))
(do ((candidate (if (odd? start) start (+ 1 start))
(+ 2 candidate))
(count 0)
(result '())
)
((= count number) (reverse result))
(if (primes:prob-pr? candidate iter)
(begin
(set! count (1+ count))
(set! result (cons candidate result)))
))))
;; Is n probably prime? First we check for divisibility by small
;; primes; if it passes that, and it's less than the maximum small
;; prime squared, we try Rabin-Miller.
;;
(define (primes:prob-pr? n count)
(and (not (primes:dbsp? n))
(or (< n (* primes:max-small-prime primes:max-small-prime))
(primes:rm-prime? n count))))
;; Is `n' Divisible By a Small Prime?
;;
(define primes:dbsp?
(let ((sqrt (cond ((provided? 'inexact) sqrt)
(else (require 'root) integer-sqrt))))
(lambda (n)
(let ((limit (min (sqrt n) primes:max-small-prime))
(divisible #f)
)
(do ((i 0 (1+ i)))
((let* ((divisor (vector-ref primes:small-primes i)))
(set! divisible (= (modulo n divisor) 0))
(or divisible (> divisor limit)))
divisible)
)))))
;; Does `n' pass the R.-M. primality test for `m' random numbers?
;;
(define (primes:rm-prime? n m)
(do ((i 0 (1+ i))
(x (+ 2 (random (- n 2) primes:prngs))))
((or (= i m) (primes:rm-composite? n x))
(= i m))))
;; Does `x' prove `n' composite using Rabin-Miller?
;;
(define (primes:rm-composite? n x)
(let ((f (primes:extract2s (- n 1))))
(primes:rm-comp? n (cdr f) (car f) x)))
;; Is `n' (where n-1 = 2^k * q) proven composite by `x'?
;;
(define (primes:rm-comp? n q k x)
(let ((y (primes:expt-mod x q n)))
(if (= y 1)
#f
(let loop ((j 0) (y y))
(cond ((= j k) #t)
((= y (- n 1)) #f)
((= y 1) #t)
(else (loop (1+ j) (primes:expt-mod y 2 n)))
)))))
;; Extract factors of 2; that is, factor x as 2^k * q
;; and return (k . q)
;;
(define (primes:extract2s x)
(do ((k 0 (1+ k))
(q x (quotient q 2)))
((odd? q) (cons k q))
))
;; Raise `base' to the power `exp' modulo `modulus' Could use the
;; modulo package, but we only need this function (and besides, this
;; implementation is quite a bit faster).
;;
(define (primes:expt-mod base exp modulus)
(do ((y 1)
(k exp (quotient k 2))
(z base (modulo (* z z) modulus)))
((= k 0) y)
(if (odd? k)
(set! y (modulo (* y z) modulus)))
))
(define primes:prngs
(make-random-state "repeatable seed for primes"))
;; This table seems big enough so that making it larger really
;; doesn't have much effect.
;;
(define primes:max-small-prime 997)
(define primes:small-primes
'#( 2 3 5 7 11 13 17 19 23 29
31 37 41 43 47 53 59 61 67 71
73 79 83 89 97 101 103 107 109 113
127 131 137 139 149 151 157 163 167 173
179 181 191 193 197 199 211 223 227 229
233 239 241 251 257 263 269 271 277 281
283 293 307 311 313 317 331 337 347 349
353 359 367 373 379 383 389 397 401 409
419 421 431 433 439 443 449 457 461 463
467 479 487 491 499 503 509 521 523 541
547 557 563 569 571 577 587 593 599 601
607 613 617 619 631 641 643 647 653 659
661 673 677 683 691 701 709 719 727 733
739 743 751 757 761 769 773 787 797 809
811 821 823 827 829 839 853 857 859 863
877 881 883 887 907 911 919 929 937 941
947 953 967 971 977 983 991 997 ))
(define primes< primes:primes<)
(define primes> primes:primes>)
(define probably-prime? primes:probably-prime?)
(provide 'primes)
|
