diff options
Diffstat (limited to 'primes.scm')
| -rw-r--r-- | primes.scm | 187 |
1 files changed, 0 insertions, 187 deletions
diff --git a/primes.scm b/primes.scm deleted file mode 100644 index 672e899..0000000 --- a/primes.scm +++ /dev/null @@ -1,187 +0,0 @@ -;; "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) |