diff options
Diffstat (limited to 'primes.scm')
| -rw-r--r-- | primes.scm | 181 |
1 files changed, 181 insertions, 0 deletions
diff --git a/primes.scm b/primes.scm new file mode 100644 index 0000000..a27b240 --- /dev/null +++ b/primes.scm @@ -0,0 +1,181 @@ +;; "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? n) + (let ((limit (min (sqrt n) primes:max-small-prime)) + (divisible #f) + ) + (do ((i 0 (1+ i))) + ((let* ((divisor (array-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))))) + ((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))) + )) + +;; 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) |