home *** CD-ROM | disk | FTP | other *** search
/ HAKERIS 11 / HAKERIS 11.ISO / linux / system / LinuxConsole 0.4 / linuxconsole0.4install-en.iso / guile0.4.lcm / share / guile / slib / factor.scm < prev    next >
Encoding:
Text File  |  2004-01-06  |  8.8 KB  |  246 lines
  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)
  246.