summaryrefslogtreecommitdiffstats
path: root/primes.scm
diff options
context:
space:
mode:
Diffstat (limited to 'primes.scm')
-rw-r--r--primes.scm181
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)