summaryrefslogtreecommitdiffstats
path: root/factor.scm
diff options
context:
space:
mode:
authorBryan Newbold <bnewbold@robocracy.org>2017-02-20 00:05:28 -0800
committerBryan Newbold <bnewbold@robocracy.org>2017-02-20 00:05:28 -0800
commitbd9733926076885e3417b74de76e4c9c7bc56254 (patch)
tree2c99dced547d48407ad44cb0e45e31bb4d02ce43 /factor.scm
parentfa3f23105ddcf07c5900de47f19af43d1db1b597 (diff)
downloadslib-bd9733926076885e3417b74de76e4c9c7bc56254.tar.gz
slib-bd9733926076885e3417b74de76e4c9c7bc56254.zip
Import Upstream version 2c7upstream/2c7
Diffstat (limited to 'factor.scm')
-rw-r--r--factor.scm199
1 files changed, 145 insertions, 54 deletions
diff --git a/factor.scm b/factor.scm
index a5c7b44..f10f0d5 100644
--- a/factor.scm
+++ b/factor.scm
@@ -1,5 +1,5 @@
-;;;; "factor.scm", prime test and factorization for Scheme
-;;; Copyright (C) 1991, 1992, 1993 Aubrey Jaffer.
+;;;; "factor.scm" factorization, prime test and generation
+;;; Copyright (C) 1991, 1992, 1993, 1998 Aubrey Jaffer.
;
;Permission to copy this software, to redistribute it, and to use it
;for any purpose is granted, subject to the following restrictions and
@@ -17,8 +17,31 @@
;promotional, or sales literature without prior written consent in
;each case.
-(require 'random)
+(require 'common-list-functions)
(require 'modular)
+(require 'random)
+(require 'byte)
+
+;;@body
+;;@0 is the random-state (@pxref{Random Numbers}) used by these
+;;procedures. If you call these procedures from more than one thread
+;;(or from interrupt), @code{random} may complain about reentrant
+;;calls.
+(define prime:prngs
+ (make-random-state "repeatable seed for primes"))
+
+
+;;@emph{Note:} The prime test and generation procedures implement (or
+;;use) the Solovay-Strassen primality test. See
+;;
+;;@itemize @bullet
+;;@item Robert Solovay and Volker Strassen,
+;;@cite{A Fast Monte-Carlo Test for Primality},
+;;SIAM Journal on Computing, 1977, pp 84-85.
+;;@end itemize
+
+;;; Solovay-Strassen Prime Test
+;;; if n is prime, then J(a,n) is congruent mod n to a**((n-1)/2)
;;; (modulo p 16) is because we care only about the low order bits.
;;; The odd? tests are inline of (expt -1 ...)
@@ -35,41 +58,104 @@
(if (odd? (quotient (- (* qq qq) 1) 8))
(- (prime:jacobi-symbol (quotient p 2) q))
(prime:jacobi-symbol (quotient p 2) q))))))
+;;@args p q
+;;Returns the value (+1, @minus{}1, or 0) of the Jacobi-Symbol of
+;;exact non-negative integer @1 and exact positive odd integer @2.
+(define jacobi-symbol prime:jacobi-symbol)
-;;;; Solovay-Strassen Prime Test
-;;; if n is prime, then J(a,n) is congruent mod n to a**((n-1)/2)
-
-;;; See:
-;;; Robert Solovay and Volker Strassen,
-;;; "A Fast Monte-Carlo Test for Primality,"
-;;; SIAM Journal on Computing, 1977, pp 84-85.
+;;@body
+;;@0 the maxinum number of iterations of Solovay-Strassen that will
+;;be done to test a number for primality.
+(define prime:trials 30)
;;; checks if n is prime. Returns #f if not prime. #t if (probably) prime.
;;; probability of a mistake = (expt 2 (- prime:trials))
;;; choosing prime:trials=30 should be enough
-(define prime:trials 30)
-;;; prime:product is a product of small primes.
-(define prime:product
- (let ((p 210))
- (for-each (lambda (s)
- (set! s (string->number s))
- (set! p (or (and s (exact? s) s) p)))
- '("2310" "30030" "510510" "9699690" "223092870"
- "6469693230" "200560490130"))
- p))
+(define (Solovay-Strassen-prime? n)
+ (do ((i prime:trials (- i 1))
+ (a (+ 2 (random (- n 2) prime:prngs))
+ (+ 2 (random (- n 2) prime:prngs))))
+ ((not (and (positive? i)
+ (= (gcd a n) 1)
+ (= (modulo (prime:jacobi-symbol a n) n)
+ (modular:expt n a (quotient (- n 1) 2)))))
+ (if (positive? i) #f #t))))
+
+;;; prime:products are products of small primes.
+(define (primes-gcd? n comps)
+ (comlist:notevery (lambda (prd) (= 1 (gcd n prd))) comps))
+(define prime:prime-sqr 121)
+(define prime:products '(105))
+(define prime:sieve (bytes 0 0 1 1 0 1 0 1 0 0 0))
+(letrec ((lp (lambda (comp comps primes nexp)
+ (cond ((< comp (quotient most-positive-fixnum nexp))
+ (let ((ncomp (* nexp comp)))
+ (lp ncomp comps
+ (cons nexp primes)
+ (next-prime nexp (cons ncomp comps)))))
+ ((< (quotient comp nexp) (* nexp nexp))
+ (set! prime:prime-sqr (* nexp nexp))
+ (set! prime:sieve (make-bytes nexp 0))
+ (for-each (lambda (prime)
+ (byte-set! prime:sieve prime 1))
+ primes)
+ (set! prime:products (reverse (cons comp comps))))
+ (else
+ (lp nexp (cons comp comps)
+ (cons nexp primes)
+ (next-prime nexp (cons comp comps)))))))
+ (next-prime (lambda (nexp comps)
+ (set! comps (reverse comps))
+ (do ((nexp (+ 2 nexp) (+ 2 nexp)))
+ ((not (primes-gcd? nexp comps)) nexp)))))
+ (lp 3 '() '(2 3) 5))
(define (prime:prime? n)
(set! n (abs n))
- (cond ((<= n 36) (and (memv n '(2 3 5 7 11 13 17 19 23 29 31)) #t))
- ((= 1 (gcd n prime:product))
- (do ((i prime:trials (- i 1))
- (a (+ 1 (random (- n 1))) (+ 1 (random (- n 1)))))
- ((not (and (positive? i)
- (= (gcd a n) 1)
- (= (modulo (prime:jacobi-symbol a n) n)
- (modular:expt n a (quotient (- n 1) 2)))))
- (if (positive? i) #f #t))))
- (else #f)))
+ (cond ((< n (bytes-length prime:sieve)) (positive? (byte-ref prime:sieve n)))
+ ((even? n) #f)
+ ((primes-gcd? n prime:products) #f)
+ ((< n prime:prime-sqr) #t)
+ (else (Solovay-Strassen-prime? n))))
+;;@args n
+;;Returns @code{#f} if @1 is composite; @code{#t} if @1 is prime.
+;;There is a slight chance @code{(expt 2 (- prime:trials))} that a
+;;composite will return @code{#t}.
+(define prime? prime:prime?)
+(define probably-prime? prime:prime?) ;legacy
+
+(define (prime:prime< start)
+ (do ((nbr (+ -1 start) (+ -1 nbr)))
+ ((or (negative? nbr) (prime:prime? nbr))
+ (if (negative? nbr) #f nbr))))
+
+(define (prime:primes< start count)
+ (do ((cnt (+ -2 count) (+ -1 cnt))
+ (lst '() (cons prime lst))
+ (prime (prime:prime< start) (prime:prime< prime)))
+ ((or (not prime) (negative? cnt))
+ (if prime (cons prime lst) lst))))
+;;@args start count
+;;Returns a list of the first @2 prime numbers less than
+;;@1. If there are fewer than @var{count} prime numbers
+;;less than @var{start}, then the returned list will have fewer than
+;;@var{start} elements.
+(define primes< prime:primes<)
+
+(define (prime:prime> start)
+ (do ((nbr (+ 1 start) (+ 1 nbr)))
+ ((prime:prime? nbr) nbr)))
+
+(define (prime:primes> start count)
+ (set! start (max 0 start))
+ (do ((cnt (+ -2 count) (+ -1 cnt))
+ (lst '() (cons prime lst))
+ (prime (prime:prime> start) (prime:prime> prime)))
+ ((negative? cnt)
+ (reverse (cons prime lst)))))
+;;@args start count
+;;Returns a list of the first @2 prime numbers greater than @1.
+(define primes> prime:primes>)
;;;;Lankinen's recursive factoring algorithm:
;From: ld231782@longs.LANCE.ColoState.EDU (L. Detweiler)
@@ -81,7 +167,7 @@
; | f(u,v,2b,n/2) or f(u+b,v+b,2b,(n-u-v-b)/2) if n even
;Thm: f(1,1,2,(m-1)/2) = (p,q) iff pq=m for odd m.
-
+
;It may be illuminating to consider the relation of the Lankinen function in
;a `computational hierarchy' of other factoring functions.* Assumptions are
;made herein on the basis of conventional digital (binary) computers. Also,
@@ -89,7 +175,7 @@
;be factored is prime). However, all algorithms would probably perform to
;the same constant multiple of the given orders for complete composite
;factorizations.
-
+
;Thm: Eratosthenes' Sieve is very roughtly O(ln(n)/n) in time and
; O(n*log2(n)) in space.
;Pf: It works with all prime factors less than n (about ln(n)/n by the prime
@@ -101,7 +187,7 @@
;Pf: It tests all odd factors less than the square root of n (about
; sqrt(n)/2), with log2(n) time for each division. It requires only
; log2(n) space for the number and divisors.
-
+
;Thm: Lankinen's algorithm is O(sqrt(n)/2) in time and O((sqrt(n)/2)*log2(n))
; in space.
;Pf: The algorithm is easily modified to seach only for factors p<q for all
@@ -129,26 +215,31 @@
(or (prime:f (+ u b) v (+ b b) (quotient (- n v) 2))
(prime:f u (+ v b) (+ b b) (quotient (- n u) 2))))))
-(define (prime:factor m)
- (case m
- ((-1 0 1) (list m))
- (else
- (if (negative? m) (cons -1 (prime:factor (- m)))
- (let* ((s (gcd m prime:product))
- (r (quotient m s)))
- (if (even? s)
- (append
- (if (= 1 r) '() (prime:factor r))
- (cons 2 (let ((s/2 (quotient s 2)))
- (if (= s/2 1) '()
- (or (prime:f 1 1 2 (quotient (- s/2 1) 2))
- (list s/2))))))
- (if (= 1 s) (or (prime:f 1 1 2 (quotient (- m 1) 2)) (list m))
- (append
- (if (= 1 r) '()
- (or (prime:f 1 1 2 (quotient (- r 1) 2)) (list r)))
- (or (prime:f 1 1 2 (quotient (- s 1) 2)) (list s))))))))))
+(define (prime:fo m)
+ (let* ((s (gcd m (car prime:products)))
+ (r (quotient m s)))
+ (if (= 1 s)
+ (or (prime:f 1 1 2 (quotient (- m 1) 2)) (list m))
+ (append
+ (if (= 1 r) '()
+ (or (prime:f 1 1 2 (quotient (- r 1) 2)) (list r)))
+ (or (prime:f 1 1 2 (quotient (- s 1) 2)) (list s))))))
-(define jacobi-symbol prime:jacobi-symbol)
-(define prime? prime:prime?)
+(define (prime:fe m)
+ (if (even? m)
+ (cons 2 (prime:fe (quotient m 2)))
+ (if (eqv? 1 m)
+ '()
+ (prime:fo m))))
+
+(define (prime:factor k)
+ (case k
+ ((-1 0 1) (list k))
+ (else (if (negative? k)
+ (cons -1 (prime:fe (- k)))
+ (prime:fe k)))))
+;;@args k
+;;Returns a list of the prime factors of @1. The order of the
+;;factors is unspecified. In order to obtain a sorted list do
+;;@code{(sort! (factor @var{k}) <)}.
(define factor prime:factor)