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
|
;;;; "factor.scm", prime test and factorization for Scheme
;;; Copyright (C) 1991, 1992, 1993 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 'random)
(require 'modular)
;;; (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))))))
;;;; 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.
;;; 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 (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)))
;;;;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: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 jacobi-symbol prime:jacobi-symbol)
(define prime? prime:prime?)
(define factor prime:factor)
|
