module/slib/factor.scm

   1 ;;;; "factor.scm" factorization, prime test and generation
   2 ;;; Copyright (C) 1991, 1992, 1993, 1998 Aubrey Jaffer.
   3 ;
   4 ;Permission to copy this software, to redistribute it, and to use it
   5 ;for any purpose is granted, subject to the following restrictions and
   6 ;understandings.
   7 ;
   8 ;1.  Any copy made of this software must include this copyright notice
   9 ;in full.
  10 ;
  11 ;2.  I have made no warrantee or representation that the operation of
  12 ;this software will be error-free, and I am under no obligation to
  13 ;provide any services, by way of maintenance, update, or otherwise.
  14 ;
  15 ;3.  In conjunction with products arising from the use of this
  16 ;material, there shall be no use of my name in any advertising,
  17 ;promotional, or sales literature without prior written consent in
  18 ;each case.
  19
  20 (require 'common-list-functions)
  21 (require 'modular)
  22 (require 'random)
  23 (require 'byte)
  24
  25 ;;@body
  26 ;;@0 is the random-state (@pxref{Random Numbers}) used by these
  27 ;;procedures.  If you call these procedures from more than one thread
  28 ;;(or from interrupt), @code{random} may complain about reentrant
  29 ;;calls.
  30 (define prime:prngs
  31   (make-random-state "repeatable seed for primes"))
  32
  33
  34 ;;@emph{Note:} The prime test and generation procedures implement (or
  35 ;;use) the Solovay-Strassen primality test. See
  36 ;;
  37 ;;@itemize @bullet
  38 ;;@item Robert Solovay and Volker Strassen,
  39 ;;@cite{A Fast Monte-Carlo Test for Primality},
  40 ;;SIAM Journal on Computing, 1977, pp 84-85.
  41 ;;@end itemize
  42
  43 ;;; Solovay-Strassen Prime Test
  44 ;;;   if n is prime, then J(a,n) is congruent mod n to a**((n-1)/2)
  45
  46 ;;; (modulo p 16) is because we care only about the low order bits.
  47 ;;; The odd? tests are inline of (expt -1 ...)
  48
  49 (define (prime:jacobi-symbol p q)
  50   (cond ((zero? p) 0)
  51         ((= 1 p) 1)
  52         ((odd? p)
  53          (if (odd? (quotient (* (- (modulo p 16) 1) (- q 1)) 4))
  54              (- (prime:jacobi-symbol (modulo q p) p))
  55              (prime:jacobi-symbol (modulo q p) p)))
  56         (else
  57          (let ((qq (modulo q 16)))
  58            (if (odd? (quotient (- (* qq qq) 1) 8))
  59                (- (prime:jacobi-symbol (quotient p 2) q))
  60                (prime:jacobi-symbol (quotient p 2) q))))))
  61 ;;@args p q
  62 ;;Returns the value (+1, @minus{}1, or 0) of the Jacobi-Symbol of
  63 ;;exact non-negative integer @1 and exact positive odd integer @2.
  64 (define jacobi-symbol prime:jacobi-symbol)
  65
  66 ;;@body
  67 ;;@0 the maxinum number of iterations of Solovay-Strassen that will
  68 ;;be done to test a number for primality.
  69 (define prime:trials 30)
  70
  71 ;;; checks if n is prime.  Returns #f if not prime. #t if (probably) prime.
  72 ;;;   probability of a mistake = (expt 2 (- prime:trials))
  73 ;;;     choosing prime:trials=30 should be enough
  74 (define (Solovay-Strassen-prime? n)
  75   (do ((i prime:trials (- i 1))
  76        (a (+ 2 (random (- n 2) prime:prngs))
  77           (+ 2 (random (- n 2) prime:prngs))))
  78       ((not (and (positive? i)
  79                  (= (gcd a n) 1)
  80                  (= (modulo (prime:jacobi-symbol a n) n)
  81                     (modular:expt n a (quotient (- n 1) 2)))))
  82        (if (positive? i) #f #t))))
  83
  84 ;;; prime:products are products of small primes.
  85 (define (primes-gcd? n comps)
  86   (comlist:notevery (lambda (prd) (= 1 (gcd n prd))) comps))
  87 (define prime:prime-sqr 121)
  88 (define prime:products '(105))
  89 (define prime:sieve (bytes 0 0 1 1 0 1 0 1 0 0 0))
  90 (letrec ((lp (lambda (comp comps primes nexp)
  91                (cond ((< comp (quotient most-positive-fixnum nexp))
  92                       (let ((ncomp (* nexp comp)))
  93                         (lp ncomp comps
  94                             (cons nexp primes)
  95                             (next-prime nexp (cons ncomp comps)))))
  96                      ((< (quotient comp nexp) (* nexp nexp))
  97                       (set! prime:prime-sqr (* nexp nexp))
  98                       (set! prime:sieve (make-bytes nexp 0))
  99                       (for-each (lambda (prime)
 100                                   (byte-set! prime:sieve prime 1))
 101                                 primes)
 102                       (set! prime:products (reverse (cons comp comps))))
 103                      (else
 104                       (lp nexp (cons comp comps)
 105                           (cons nexp primes)
 106                           (next-prime nexp (cons comp comps)))))))
 107          (next-prime (lambda (nexp comps)
 108                        (set! comps (reverse comps))
 109                        (do ((nexp (+ 2 nexp) (+ 2 nexp)))
 110                            ((not (primes-gcd? nexp comps)) nexp)))))
 111   (lp 3 '() '(2 3) 5))
 112
 113 (define (prime:prime? n)
 114   (set! n (abs n))
 115   (cond ((< n (bytes-length prime:sieve)) (positive? (byte-ref prime:sieve n)))
 116         ((even? n) #f)
 117         ((primes-gcd? n prime:products) #f)
 118         ((< n prime:prime-sqr) #t)
 119         (else (Solovay-Strassen-prime? n))))
 120 ;;@args n
 121 ;;Returns @code{#f} if @1 is composite; @code{#t} if @1 is prime.
 122 ;;There is a slight chance @code{(expt 2 (- prime:trials))} that a
 123 ;;composite will return @code{#t}.
 124 (define prime? prime:prime?)
 125 (define probably-prime? prime:prime?)   ;legacy
 126
 127 (define (prime:prime< start)
 128   (do ((nbr (+ -1 start) (+ -1 nbr)))
 129       ((or (negative? nbr) (prime:prime? nbr))
 130        (if (negative? nbr) #f nbr))))
 131
 132 (define (prime:primes< start count)
 133   (do ((cnt (+ -2 count) (+ -1 cnt))
 134        (lst '() (cons prime lst))
 135        (prime (prime:prime< start) (prime:prime< prime)))
 136       ((or (not prime) (negative? cnt))
 137        (if prime (cons prime lst) lst))))
 138 ;;@args start count
 139 ;;Returns a list of the first @2 prime numbers less than
 140 ;;@1.  If there are fewer than @var{count} prime numbers
 141 ;;less than @var{start}, then the returned list will have fewer than
 142 ;;@var{start} elements.
 143 (define primes< prime:primes<)
 144
 145 (define (prime:prime> start)
 146   (do ((nbr (+ 1 start) (+ 1 nbr)))
 147       ((prime:prime? nbr) nbr)))
 148
 149 (define (prime:primes> start count)
 150   (set! start (max 0 start))
 151   (do ((cnt (+ -2 count) (+ -1 cnt))
 152        (lst '() (cons prime lst))
 153        (prime (prime:prime> start) (prime:prime> prime)))
 154       ((negative? cnt)
 155        (reverse (cons prime lst)))))
 156 ;;@args start count
 157 ;;Returns a list of the first @2 prime numbers greater than @1.
 158 (define primes> prime:primes>)
 159
 160 ;;;;Lankinen's recursive factoring algorithm:
 161 ;From: ld231782@longs.LANCE.ColoState.EDU (L. Detweiler)
 162
 163 ;                  |  undefined if n<0,
 164 ;                  |  (u,v) if n=0,
 165 ;Let f(u,v,b,n) := | [otherwise]
 166 ;                  |  f(u+b,v,2b,(n-v)/2) or f(u,v+b,2b,(n-u)/2) if n odd
 167 ;                  |  f(u,v,2b,n/2) or f(u+b,v+b,2b,(n-u-v-b)/2) if n even
 168
 169 ;Thm: f(1,1,2,(m-1)/2) = (p,q) iff pq=m for odd m.
 170
 171 ;It may be illuminating to consider the relation of the Lankinen function in
 172 ;a `computational hierarchy' of other factoring functions.*  Assumptions are
 173 ;made herein on the basis of conventional digital (binary) computers.  Also,
 174 ;complexity orders are given for the worst case scenarios (when the number to
 175 ;be factored is prime).  However, all algorithms would probably perform to
 176 ;the same constant multiple of the given orders for complete composite
 177 ;factorizations.
 178
 179 ;Thm: Eratosthenes' Sieve is very roughtly O(ln(n)/n) in time and
 180 ;     O(n*log2(n)) in space.
 181 ;Pf: It works with all prime factors less than n (about ln(n)/n by the prime
 182 ;    number thm), requiring an array of size proportional to n with log2(n)
 183 ;    space for each entry.
 184
 185 ;Thm: `Odd factors' is O((sqrt(n)/2)*log2(n)) in time and O(log2(n)) in
 186 ;     space.
 187 ;Pf: It tests all odd factors less than the square root of n (about
 188 ;    sqrt(n)/2), with log2(n) time for each division.  It requires only
 189 ;    log2(n) space for the number and divisors.
 190
 191 ;Thm: Lankinen's algorithm is O(sqrt(n)/2) in time and O((sqrt(n)/2)*log2(n))
 192 ;     in space.
 193 ;Pf: The algorithm is easily modified to seach only for factors p<q for all
 194 ;    pq=m.  Then the recursive call tree forms a geometric progression
 195 ;    starting at one, and doubling until reaching sqrt(n)/2, or a length of
 196 ;    log2(sqrt(n)/2).  From the formula for a geometric progression, there is
 197 ;    a total of about 2^log2(sqrt(n)/2) = sqrt(n)/2 calls.  Assuming that
 198 ;    addition, subtraction, comparison, and multiplication/division by two
 199 ;    occur in constant time, this implies O(sqrt(n)/2) time and a
 200 ;    O((sqrt(n)/2)*log2(n)) requirement of stack space.
 201
 202 (define (prime:f u v b n)
 203   (if (<= n 0)
 204       (cond ((negative? n) #f)
 205             ((= u 1) #f)
 206             ((= v 1) #f)
 207             ; Do both of these factors need to be factored?
 208             (else (append (or (prime:f 1 1 2 (quotient (- u 1) 2))
 209                               (list u))
 210                           (or (prime:f 1 1 2 (quotient (- v 1) 2))
 211                               (list v)))))
 212       (if (even? n)
 213           (or (prime:f u v (+ b b) (quotient n 2))
 214               (prime:f (+ u b) (+ v b) (+ b b) (quotient (- n (+ u v b)) 2)))
 215           (or (prime:f (+ u b) v (+ b b) (quotient (- n v) 2))
 216               (prime:f u (+ v b) (+ b b) (quotient (- n u) 2))))))
 217
 218 (define (prime:fo m)
 219   (let* ((s (gcd m (car prime:products)))
 220          (r (quotient m s)))
 221     (if (= 1 s)
 222         (or (prime:f 1 1 2 (quotient (- m 1) 2)) (list m))
 223         (append
 224          (if (= 1 r) '()
 225              (or (prime:f 1 1 2 (quotient (- r 1) 2)) (list r)))
 226          (or (prime:f 1 1 2 (quotient (- s 1) 2)) (list s))))))
 227
 228 (define (prime:fe m)
 229   (if (even? m)
 230       (cons 2 (prime:fe (quotient m 2)))
 231       (if (eqv? 1 m)
 232           '()
 233           (prime:fo m))))
 234
 235 (define (prime:factor k)
 236   (case k
 237     ((-1 0 1) (list k))
 238     (else (if (negative? k)
 239               (cons -1 (prime:fe (- k)))
 240               (prime:fe k)))))
 241 ;;@args k
 242 ;;Returns a list of the prime factors of @1.  The order of the
 243 ;;factors is unspecified.  In order to obtain a sorted list do
 244 ;;@code{(sort! (factor @var{k}) <)}.
 245 (define factor prime:factor)