Import Upstream version 2a6upstream/2a6

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)