summaryrefslogtreecommitdiffstats
path: root/factor.scm
blob: f10f0d589c42341e58e30eaa4f5c91c5a951e2c9 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
;;;; "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
;understandings.
;
;1.  Any copy made of this software must include this copyright notice
;in full.
;
;2.  I have made no warrantee or representation that the operation of
;this software will be error-free, and I am under no obligation to
;provide any services, by way of maintenance, update, or otherwise.
;
;3.  In conjunction with products arising from the use of this
;material, there shall be no use of my name in any advertising,
;promotional, or sales literature without prior written consent in
;each case.

(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 ...)

(define (prime:jacobi-symbol p q)
  (cond ((zero? p) 0)
	((= 1 p) 1)
	((odd? p)
	 (if (odd? (quotient (* (- (modulo p 16) 1) (- q 1)) 4))
	     (- (prime:jacobi-symbol (modulo q p) p))
	     (prime:jacobi-symbol (modulo q p) p)))
	(else
	 (let ((qq (modulo q 16)))
	   (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)

;;@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 (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 (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)

;                  |  undefined if n<0,
;                  |  (u,v) if n=0,
;Let f(u,v,b,n) := | [otherwise]
;                  |  f(u+b,v,2b,(n-v)/2) or f(u,v+b,2b,(n-u)/2) if n odd
;                  |  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,
;complexity orders are given for the worst case scenarios (when the number to
;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
;    number thm), requiring an array of size proportional to n with log2(n)
;    space for each entry.

;Thm: `Odd factors' is O((sqrt(n)/2)*log2(n)) in time and O(log2(n)) in
;     space.
;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
;    pq=m.  Then the recursive call tree forms a geometric progression
;    starting at one, and doubling until reaching sqrt(n)/2, or a length of
;    log2(sqrt(n)/2).  From the formula for a geometric progression, there is
;    a total of about 2^log2(sqrt(n)/2) = sqrt(n)/2 calls.  Assuming that
;    addition, subtraction, comparison, and multiplication/division by two
;    occur in constant time, this implies O(sqrt(n)/2) time and a
;    O((sqrt(n)/2)*log2(n)) requirement of stack space.

(define (prime:f u v b n)
  (if (<= n 0)
      (cond ((negative? n) #f)
	    ((= u 1) #f)
	    ((= v 1) #f)
	    ; Do both of these factors need to be factored?
	    (else (append (or (prime:f 1 1 2 (quotient (- u 1) 2))
			      (list u))
			  (or (prime:f 1 1 2 (quotient (- v 1) 2))
			      (list v)))))
      (if (even? n)
	  (or (prime:f u v (+ b b) (quotient n 2))
	      (prime:f (+ u b) (+ v b) (+ b b) (quotient (- n (+ u v b)) 2)))
	  (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: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 (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)