[plt-scheme] Fast primality testing
Hi Joshua,
> I'm playing around with a problem from Al Zimmerman's Programming Contest,
> http://www.recmath.org/contest/description.php
> not really to try to win the contest, but just to have an excuse to
> teach myself some more scheme.
>
> One part I'm not really interested in teaching myself, though: is
> there a reasonably fast primality-checker in some library somewhere? No big requirements here: it just has to operate
to return true or
> false, pretty fast, for numbers up to 17 digits.
I have attached a file with various number theoretic functions.
Prime tests and factorization works for quite big numbers.
> (prime? 567348343)
#f
> (factorize 5673483433473648736487364)
((2 2) (3 2) (31 1) (4609093 1) (24487651 1) (45042553 1))
For modular arithmetic the form with-modulus is provided. It
evaluates it's body in an environment where + - * and ^ work
modulo the given integer. Furhermore inv computes the inverse.
; Euler disproved Fermat's assertion that all Fermat numbers were prime
> (with-modulus 641
(+ (^ 2 (^ 2 5)) 1))
0
; The eleventh Fermat number is also composite.
> (with-modulus 319489
(+ (^ 2 (^ 2 11)) 1))
0
Enter (+ (^ 2 (^ 2 11)) 1) in the repl to see the number of digits.
The function bezout computes the coeffecient from the linear
combination for the greatest common divisor.
(bezout a b) = (list u v) <=> gcd(a,b) = au + bv
It also works for more numbers than two.
> (bezout 2 3)
(-1 1)
> (bezout 6 15 35)
(1 2 -1)
> (+ (* 1 6) (* 2 15) (* -1 35))
1
> (gcd 6 5 35)
1
The function solve-chinese finds the coefficient in the Chinese
Remainder Theorem
(solve-chinese (a1 ... ak) (n1 ... nk)) = x
<=> x = a1 mod n1, ..., x=ak mod nk
>(solve-chinese (2 3 2) (3 5 7))
23
Here is list of the useful functions. See the code for theorems and
explanations.
;;; DIVISORS
divides? a b
coprime? a . bs) ; aka relatively-prime?
pairwise-coprime? ; aka pairwise-relatively-prime?
positive-divisors n
positive-divisors-of-n-below n m
;;; INTEGER ROOTS
is-nth-root a n candidate
integer-root a n
as-power a
perfect-power? a
perfect-power a
;;; MODULAR ARITHMETIC
bezout a . bs
solve-chinese as ns
inverse a n
syntax: (with-modulus n body ...)
evaluate body where + - * / ^ mod inv
behaves like in Z modulo N
;;; PRIMES
prime? : n
next-prime n
prev-prime n
next-primes n primes-wanted
prev-primes n primes-wanted
nth-prime n
factorize n
prime-divisors/exponents n
prime-divisors n
prime-exponents n
max-dividing-power n m
prime-power n
prime-power? n
;;; RANDOM
big-random n
random-integer from to)
;;; SPECIAL NUMBERS
cyclotomic p x
fermat n
mersenne p
;;; NUMBER THEORETIC FUNCTIONS
phi (Euler phi)
legendre a p
quadratic-residue? a n
;;; PRIMITIVE ROOTS
unit-group n
order g n
orders n
primitive-root? g n
exists-primitive-root? n
primitive-root? g n
primitive-root n
primitive-roots n
primitive-root n
;;; ADVANCED STUFF
pollard n
p-adic-newton-iteration phi Dphi p l g0 s0
;;; CONVENIENCES
interval m n
square z
--
Jens Axel Søgaard
;;; number-theory.scm -- Jens Axel Soegaard
(require (lib "27.ss" "srfi"))
;;; UTILIIES
; (interval m n) -> (list m m+1 ... n)
(define (interval m n)
(if (> m n)
'()
(cons m (interval (+ m 1) n))))
; square : number -> number
(define (square z)
(* z z))
(define (boolify x)
(if x #t #f))
;; big-random : N -> N
;; return random integer in the interval [0;n[
;(define (big-random n)
; (let ([l 30])
; (let ([bits-needed (ceiling (/ (log n) (log 2)))]
; [M (expt 2 l)])
; (let loop ([blocks (quotient bits-needed l)]
; [r (random (inexact->exact (expt 2 (remainder bits-needed l))))])
; (if (= blocks 0)
; (remainder r n)
; (loop (- blocks 1) (+ (* r M) (random M))))))))
(define big-random random-integer)
(define (random-integer-in-interval from to)
; return random integer from the half-open
; interval [from;to)
(+ from (random-integer (- to from))))
;;;;; START HERE
; THEOREM (Division algorithm)
; If a and b in Z with b<>0 then there exists
; unique integers q and r such that
; a = qb+r and 0<=r<|b|
; q is called the quotient and r the remainder.
; SCHEME
; q = (quotient a b)
; r = (remainder a b)
; DEF (Factor, divides)
; If a and b are integers, and b=qa for some q,
; then a divides b. In that case we write a|b,
; and a is called a factor of b.
; (divides? a b) <=> a|b
(define (divides? a b)
(= (remainder b a) 0))
; DEF (Common divisor, Greatest common divisor, gcd)
; If d|a and d|b then d is a common divisor of a and b.
; The greatest common divisior d satisfies
; 1) d|a and d|b
; 2) c|a and c|b => d|c
; Note, there is a unique positive greatest common divisor
; if a and b are bot zero.
; SCHEME
; d = (gcd a b) <=> d is a positive greatest common divisor
; d = (gcd a b c) <=> d is a positive gcd of a, b and c
; LEMMA
; a=qb+r => gcd(a,b)=gcd(b,r)
; ALGORITHM (Euclid)
; The lemma leads to Euclid's algorithm for computation of gcd.
(define (gcd-euclid a b)
(define (loop a b) ; a>=b>0
(let ([r (remainder a b)])
(if (= r 0)
b
(loop b r))))
(loop (max (abs a) (abs b))
(min (abs a) (abs b))))
; THEOREM (Bezout's identity)
; If a and b are integers (not both zero), then
; there exists integers u and v such that
; gcd(a,b) = au + bv
; Note: u and v are not unique
; (bezout-binary a b) = (list u v) <=> gcd(a,b) = au + bv
(define (bezout-binary a b)
(define (loop a b ua va ub vb) ; a>=b>0 , a = ua*a+ub*b, b = ub*a+ub*b
; (display (list a b ua va ub vb)) (newline)
(let ([r (remainder a b)]
[q (quotient a b)])
(if (= r 0)
(list ub vb)
(loop b r ub vb (- ua (* q ub)) (- va (* q vb))))))
(if (> a b)
(loop a b 1 0 0 1)
(loop b a 0 1 1 0)))
; This shows how bezout below was constructed
;(define (bezout/3 a b c)
; (let ([uv (bezout-binary a b)]
; [st (bezout-binary (gcd a b) c)])
; (display st)
; (let ([u (first uv)]
; [v (second uv)]
; [s (first st)]
; [t (second st)])
; (list (* s u) (* s v) t))))
; (bezout-binary a b c ...) -> (list u v w ...) <=> gcd(a,b,c,...) = au + bv + cw + ...
(define (bezout a . bs)
(if (null? bs)
(list 1)
(let ([uvs (apply bezout bs)]
[st (bezout-binary (apply gcd bs) a)])
(let ([s (first st)]
[t (second st)])
(cons t (map (lambda (u) (* s u))
uvs))))))
; DEF (Coprime, relatively prime)
; Two or more integers are called coprime, if their greatest common divisor is 1.
; a, b and c coprime <=> gcd(a,b,c)=1
(define (coprime? a . bs)
(= 1 (apply gcd (cons a bs))))
(define (pairwise-coprime? a . bs)
(cond
[(null? bs) #t]
[else (and (andmap (lambda (b)
(coprime? a b))
bs)
(apply pairwise-coprime? bs))]))
(define relatively-prime? coprime?)
(define pairwiae-relatively-prime? pairwise-coprime?)
; DEF (Prime)
; An integer p>1 is a prime if the only positive divisors of p are 1 and p.
;(define (prime? p)
; (and (> p 1)
; (= (length (positive-divisors p)) 2)))
(define (odd-prime? n)
(and (odd? n) (prime? n)))
; LEMMA An integer n>1 is composite <=> exists prime p s.t. p|n and p<=sqrt(n)
;(define (prime? p)
; (and (> p 1)
; (= 1 (length (positive-divisors-of-n-below p (inexact->exact (floor (sqrt p))))))))
;;;; TODO: USE THE FACTORIZE TO COMPUTE THESE
(define (positive-divisors n)
(positive-divisors-of-n-below n n))
(define (positive-divisors-of-n-below n m)
(cond
[(<= m 0) empty]
[(divides? m n) (cons m (positive-divisors-of-n-below n (- m 1)))]
[else (positive-divisors-of-n-below n (- m 1))]))
; DEF (Cyclotomic polynomial)
; The p-th cyclotomic polynomial is
; phi_p(x) = 1 + x + x^2 + ... + x^(p-1)
; THEOREM phi_p is irreducible if p is prime
; PROOF Apply Eisenstein's criterion
(define (cyclotomic p x)
(if (= p 1)
1
(+ 1 (* x (cyclotomic (- p 1) x)))))
; (cyclotomic 2 3) = 1+3 = 4
; THEOREM (Fundamental Theorem of Arithmetic)
; An integer n>1 has a factorisation into prime-powers
; e1 ek
; n = p1 ... pk , where pi is prime and ei>0
; (prime-divisors n) -> (list p1 p2 ... pk)
;(define (prime-divisors n)
; (reverse (filter prime? (positive-divisors n))))
; (prime-exponents n) -> (list e1 e2 ... ek)
;(define (prime-exponents n)
; (map (lambda (p)
; (max-dividing-power p n))
; (prime-divisors n)))
; (max-dividing-power p n) = m <=> p^m | n and p^(m+1) doesn't divide n
(define (max-dividing-power p n)
(define (loop p-to-e e)
(if (divides? p-to-e n)
(loop (* p p-to-e) (+ e 1))
(- e 1)))
(loop 1 0))
; (prime-divisors/exponents n) -> '((p1 e1) (p2 e2) ... (pk ek))
;(define (prime-divisors/exponents n)
; (map list
; (prime-divisors n)
; (prime-exponents n)))
;(define prime-divisors/exponents factorize)
; prime-power : N -> (union #f (list prime exponent))
; if n is a prime power, return list of prime and exponent in question,
; otherwise return #f
(define (prime-power n)
(let ([factorization (prime-divisors/exponents n)])
(if (= (length factorization) 1)
(first (prime-divisors/exponents n))
#f)))
(define (prime-power? n)
(boolify (prime-power n)))
(define (odd-prime-power? n)
(let ([p/e (prime-power n)])
(and p/e
(odd? (first p/e)))))
; (nth-prime n) = pn , where p1=2, p2=3, p3=5, ...
(define (nth-prime n)
(define (next c) (+ c 2))
(define (loop candidate m)
(if (prime? candidate)
(if (= m n)
candidate
(loop (next candidate) (+ m 1)))
(loop (next candidate) m)))
(if (= n 1)
2
(loop 3 2)))
; DEF (Fermat numbers)
; The n'th Fermat number is
; Fn = 2^(2^n) + 1
; I.e. F1=5, F2=17, F3=257, ...
; (fermat n) = Fn
(define (fermat n)
(+ 1 (expt 2 (expt 2 n))))
; DEF (Mersenne numbers)
; If p is a prime then 2^p-1 is called a Mersenne number.
(define (mersenne p)
(- (expt 2 p) 1))
; MODULAR ARITHMETIC
; THEOREM
; If gcd(a,n)=1 then there exist b such that
; ab=1 mod n
; The number b is called an inverse of a modulo n.
; (inverse a n) = b , where ab=1 mod n and b in {0,...,n-1}
(define (inverse a n)
(if (coprime? a n)
(modulo (first (bezout a n)) n)
#f))
; Within a (with-modulus n form1 ...) the return values of
; the arithmetival operations +, -, * and ^ are automatically
; reduced modulo n. Furthermore (mod x)=(modulo x n) and
; (inv x)=(inverse x n).
; Example: (with-modulus 3 (^ 2 4)) ==> 1
(define-syntax with-modulus (lambda (stx)
(syntax-case stx ()
[(with-modulus e form ...)
(with-syntax ([+ (datum->syntax-object (syntax with-modulus) '+)]
[- (datum->syntax-object (syntax with-modulus) '-)]
[* (datum->syntax-object (syntax with-modulus) '*)]
[^ (datum->syntax-object (syntax with-modulus) '^)]
[mod (datum->syntax-object (syntax with-modulus) 'mod)]
[inv (datum->syntax-object (syntax with-modulus) 'inv)])
(syntax (let* ([n e]
[mod (lambda (x) (modulo x n))]
[inv (lambda (x) (inverse x n))]
[+ (compose mod +)]
[- (compose mod -)]
[* (compose mod *)]
[square (lambda (x) (* x x))]
[^ (rec ^ (lambda (a b)
(cond
[(= b 0) 1]
[(even? b) (square (^ a (/ b 2)))]
[else (* a (^ a (sub1 b)))])))])
form ...)))])))
; THEOREM (The Chinese Remainder Theorem)
; Let n1,...,nk be positive integers with gcd(ni,nj)=1 whenever i<>j,
; and let a1,...,ak be any integers. Then the solutions to
; x=a1 mod n1, ..., x=ak mod nk
; has a single solution in {0,...,n-1}, where n=n1*...nk.
; Example : (solve-chinese '(2 3 2) '(3 5 7)) = 23
(define (solve-chinese as ns)
; the ns should be coprime
(let* ([n (apply * ns)]
[cs (map (lambda (ni) (/ n ni)) ns)]
[ds (map inverse cs ns)])
(modulo (apply + (map * as cs ds)) n)))
; THEOREM (Fermat's Little Theorem)
; If p is prime and a<>0 mod p then a^(p-1)=1 mod p
; THEOREM (Wilson's Theorem)
; n prime <=> (n-1)! = -1 mod n
; NB: prime-wilson? is very inefficient, so this is
; included only for completeness.
(define (prime-wilson? n)
(with-modulus n
(define (fac x)
(if (= 0 x)
1
(* x (fac (sub1 x)))))
(= (mod (fac (- n 1)))
(mod -1))))
;;; EULER'S FUNCTION
; DEFINITION (Euler's phi function)
; phi(n) is the number of integers a=1,2,... such that gcd(a,n)=1
; THEOREM
; If m and n are coprime then
; phi(mn) = phi(m) phi(n)
; THEOREM (Euler's phi function)
; If the prime power factorization of p is
; e1 ek
; n = p1 ... pk , where pi is prime and ei>0
; then
; k 1
; phi(n) = n * product (1 - ---- )
; i=1 pi
(define (phi n)
(* n (apply *
(map (lambda (p)
(- 1 (/ 1 p)))
(prime-divisors n)))))
;;; THE GROUP OF UNITS
; DEFINITION (Order)
; If G is a finite group with identity element e,
; then the order of g in G is the least k>0 such that
; g^k=e
; DEFINITION (Un)
; The group of units in Zn with respect to multiplication
; modulo n is called Un.
(define (unit-group n)
(filter (lambda (m)
(coprime? m n))
(interval 1 (- n 1))))
(define (order g n)
(if (not (coprime? g n))
(error "In (order g n) the g and n must me coprime")
(with-modulus n
(let loop ([k 1]
[a g])
(if (= a 1)
k
(loop (+ k 1) (* a g)))))))
(define (orders n)
(map (lambda (m) (order m n))
(unit-group n)))
; DEFINITION (Primitive Root)
; A generator g of Un is called a primitive root mod n.
; I.e. order(g)=phi(n) <=> g is primitive
(define (primitive-root? g n)
(if (not (coprime? g n))
(error "In (primitive-root? g n) the g and n must me coprime")
(= (order g n) (phi n))))
; THEOREM (Existence of primitive roots)
; Un is cyclic (i.e. have a primitive root)
; <=> n = 1, 2, 4, p^e, 2*p^e where p is an odd prime
(define (exists-primitive-root? n)
(cond
[(member n '(1 2 4)) #t]
[(odd? n) (odd-prime-power? n)]
[else (odd-prime-power? (quotient n 2))]))
; LEMMA
; a in Un is a primitive root
; <=> phi(n)/q
; a <> 1 in Un for all primes q dividing phi(n)
(define (primitive-root? g n)
(if (not (coprime? g n))
(error "In (primitive-root? g n) the g and n must me coprime")
(let ([phi-n (phi n)])
(with-modulus n
(andmap (lambda (x) x)
(map (lambda (q)
(not (= (^ g (/ phi-n q)) 1)))
(prime-divisors phi-n)))))))
; primitive-root : N -> Un
; return primitive root of n if one exists,
; otherwise return #f
(define (primitive-root n)
(and (exists-primitive-root? n)
(let* ([phi-n (phi n)]
[qs (prime-divisors phi-n)])
(define (primitive-root? g)
(with-modulus n
(andmap (lambda (x) x)
(map (lambda (q)
(not (= (^ g (/ phi-n q)) 1)))
qs))))
(let loop ([g 1])
(cond
[(= g n) #f]
[(not (coprime? g n)) (loop (+ g 1))]
[(primitive-root? g) g]
[else (loop (+ g 1))])))))
; LEMMA
; If Un has a primitive root, then it has phi(phi(n)) primitive roots
; primitive-roots : integer -> list
; return list of all primitive roots of Un
(define (primitive-roots n)
(if (not (exists-primitive-root? n))
empty
(let* ([phi-n (phi n)]
[qs (prime-divisors phi-n)])
(define (primitive-root? g)
(with-modulus n
(andmap (lambda (x) x)
(map (lambda (q)
(not (= (^ g (/ phi-n q)) 1)))
qs))))
(let loop ([g 1]
[roots empty])
(cond
[(= g n) (reverse roots)]
[(not (coprime? g n)) (loop (+ g 1) roots)]
[(primitive-root? g) (loop (+ g 1) (cons g roots))]
[else (loop (+ g 1) roots)])))))
(define (primitive-root n)
(and (exists-primitive-root? n)
(cond
; U_p^e , p odd
[(and (odd-prime-power? n)
(not (prime? n))) (let* ([p (first (prime-power n))]
[g (primitive-root p)])
(if (= (order g (* p p)) (phi (* p p)))
g
(modulo (+ g p) n)))]
; U_2p^e , p odd
[(and (even? n)
(odd-prime? (quotient n 2))) (let ([g (primitive-root (quotient n 2))])
(if (odd? g)
g
(modulo (+ g (quotient n 2)) n)))]
; General case
[else (let* ([phi-n (phi n)]
[qs (prime-divisors phi-n)])
(define (primitive-root? g)
(with-modulus n
(andmap (lambda (x) x)
(map (lambda (q)
(not (= (^ g (/ phi-n q)) 1)))
qs))))
(let loop ([g 1])
(cond
[(= g n) #f]
[(not (coprime? g n)) (loop (+ g 1))]
[(primitive-root? g) g]
[else (loop (+ g 1))])))])))
;;; QUADRATIC RESIDUES
; DEFINITION (Quadratic residue)
; a in Un is a quadratic residue,
; if exists s : a=s^2 (mod n)
; p is prime
(define (legendre a p)
(let ([l (with-modulus p
(^ a (/ (- p 1) 2)))])
(if (<= 0 l 1)
l
-1)))
(define (quadratic-residue? a n)
(let* ([ps (prime-divisors n)]
[odd-ps (if (= (first ps) 2)
(rest ps)
ps)])
(and (andmap (lambda (p)
(= (legendre a p) 1))
odd-ps)
(cond
[(divides? 8 n) (= (modulo a 8) 1)]
[(divides? 4 n) (= (modulo a 4) 1)]
[else #t]))))
;;; PRIMALITY TESTS
; THEOREM (Fermat's little theorem) [MCA,p.75]
; p
; p prime, a in Z => a = a mod p
; p-1
; (not p|a) => a = a mod p
; [MCA, p.507 - Fermat test]
(define (prime-fermat? n)
(cond
[(= n 2) #t]
[else (let* ([a (random-integer-in-interval 2 (- n 1))]
[b (with-modulus n
(^ a (- n 1)))])
(if (= b 1)
'possibly-prime
#f))]))
; The Fermat test answers 'possibly-prime, if n is a prime or a
; Carmichael number with gcd(a,n)=1.
; THEOREM 18.5 Strong pseudoprimality test
; If N is prime, then algorithm 18.5 returns "probably prime".
; If N is composite and not a Carmichael number, then the
; algorithm returns "composite" with probability at least 1/2.
; If N is a Carmichael number, the algorithm returns a
; proper divisor of N with probability at least 1/2.
; [MCA, p. 509 - Algorithm 18.5 - Strong pseudoprimality test]
(define (prime-strong-pseudo-single? n)
(cond
[(= n 2) #t]
[(= n 3) #t]
[(< n 2) #f]
[else (let* ([a (random-integer-in-interval 2 (- n 1))]
[g (gcd a n)])
(if (> g 1)
g
; 3. write n-1 = 2^ny * m , m odd
(let loop ([ny 0]
[m (- n 1)])
(if (even? m)
(loop (add1 ny) (/ m 2))
; 4. for i=1,...,ny do bi <- b_{i-1}^2 rem N
(let ([b (with-modulus n (^ a m))])
(if (= b 1)
'probably-prime
(let loop ([i 0]
[b b]
[b-old b])
(if (< i ny)
(loop (add1 i)
(with-modulus n (* b b))
b)
(if (= b 1)
(let ([g (gcd (+ b-old 1) n)])
(if (or (= g 1) (= g n))
'probably-prime
g))
'composite)))))))))]))
(define integer-log2
(let ([log2 (log 2)])
(lambda (n)
(/ (log n) log2))))
(define (big-log2 n)
(define (loop n l)
(if (<= n 1)
l
(loop (quotient n 2) (add1 l))))
(loop n 0))
(define prime-strong-pseudo-certainty 1/10000000)
(define prime-strong-pseudo-trials (inexact->exact (ceiling (integer-log2 (/ 1 prime-strong-pseudo-certainty)))))
; For large numbers we can repeat the above test to improve the probability
(define (prime-strong-pseudo/explanation n)
(if #f ;(< n 1000000)
(error "prime-strong-pseudo: n must be 'large'" n)
(let loop ([trials prime-strong-pseudo-trials]
[result (prime-strong-pseudo-single? n)])
(if (= trials 0)
'very-probably-prime
(cond
[(equal? result 'probably-prime) (loop (sub1 trials) (prime-strong-pseudo-single? n))]
[else result])))))
(define (prime-strong-pseudo? n)
(let ([explanation (prime-strong-pseudo/explanation n)])
(if (or (equal? explanation 'very-probably-prime)
(equal? explanation #t))
#t
#f)))
(define prime? prime-strong-pseudo?)
(define (next-prime n)
(cond
[(= n 2) 3]
[(even? n) (let ([n+1 (add1 n)])
(if (prime? n+1)
n+1
(next-prime n+1)))]
[else (let ([n+2 (+ n 2)])
(if (prime? n+2)
n+2
(next-prime n+2)))]))
(define (prev-prime n)
(cond
[(= n 3) 2]
[(even? n) (let ([n-1 (sub1 n)])
(if (prime? n-1)
n-1
(prev-prime n-1)))]
[else (let ([n-2 (- n 2)])
(if (prime? n-2)
n-2
(prev-prime n-2)))]))
(define (next-primes n primes-wanted)
(if (= primes-wanted 0)
'()
(let ([next (next-prime n)])
(cons next (next-primes next (sub1 primes-wanted))))))
(define (prev-primes n primes-wanted)
(if (= primes-wanted 0)
'()
(let ([prev (prev-prime n)])
(cons prev (prev-primes prev (sub1 primes-wanted))))))
;;; ALGORITHM 19.6 Floyd's cycle detection trick
; [MCA, p. 536]
; The function f : alpha -> alpha generates a sequence
; given x0 in alpa, by x0, x1=f(x0), x2=f(x1), ...
; Floyd-detect-cycle returns an index i>0 s.t. xi = x_2i
(define (floyd-detect-cycle f x0)
(do ([xi x0 (f xi)]
[yi x0 (f (f yi))]
[i 0 (add1 i)])
[(= xi yi) i]))
;;; ALGORITHM 19.8 Pollard's rho method
; INPUT N>=3 neither a prime nor a perfect power
; OUTPUT Either a proper divisor of N or #f
(define (pollard n)
; (display "pollard: n = ") (display n) (display " , ")
(let ([x0 (big-random n)])
; (display "x0 = ") (display x0) (newline)
(do ([xi x0 (remainder (+ (* xi xi) 1) n)]
[yi x0 (remainder (+ (square (+ (* yi yi) 1)) 1) n)]
[i 0 (add1 i)]
[g 1 (gcd (- xi yi) n)])
[(or (< 1 g n)
(> i (sqrt n))
) (if (< 1 g n)
g
#f)]
;(display (list i xi yi)) (newline)
)))
(define (pollard-factorize n)
; (display "pf: ") (display n) (newline)
(cond
[(prime? n) `((, n 1))]
[(even? n) `((2 1) ,@(pollard-factorize (quotient n 2)))]
[(divides? 3 n) `((3 1) ,@(pollard-factorize (quotient n 3)))]
[(perfect-power n) => (lambda (base-and-exp)
(cond
[(prime? (car base-and-exp)) (list base-and-exp)]
[else (map (lambda (b-and-e)
(list (car b-and-e) (* (cadr base-and-exp)
(cadr b-and-e))))
(pollard-factorize (car base-and-exp)))]))]
[else (let loop ([divisor (pollard n)])
(if divisor
(append (pollard-factorize divisor)
(pollard-factorize (quotient n divisor)))
(loop (pollard n))))]))
;(define (factorize n)
; (define (combine-same-base list-of-base-and-exponents)
; (match list-of-base-and-exponents
; [() '()]
; [((b e)) list-of-base-and-exponents]
; [((b1 e1) (b2 e2) . more) (if (= b1 b2)
; (combine-same-base (cons (list b1 (+ e1 e2))
; (cdr (cdr list-of-base-and-exponents))))
; (cons (car list-of-base-and-exponents)
; (combine-same-base (cdr list-of-base-and-exponents))))]))
; (combine-same-base
; (mergesort
; (pollard-factorize n)
; (lambda (x y)
; (match (list x y)
; [((prime-x exponent-x) (prime-y expnonent-y)) (<= prime-x prime-y)]
; [((prime-x exponent-x) prime-y) (<= prime-x prime-y)]
; [(prime-x (prime-y expnonent-y)) (<= prime-x prime-y)]
; [(prime-x prime-y) (<= prime-x prime-y)])))))
(define (factorize n)
(define (combine-same-base list-of-base-and-exponents)
(let ([l list-of-base-and-exponents])
(cond
[(null? l) '()]
[(null? (cdr l)) l]
[else (let ([b1 (first (first l))]
[e1 (second (first l))]
[b2 (first (second l))]
[e2 (second (second l))]
[more (cddr l)])
(if (= b1 b2)
(combine-same-base (cons (list b1 (+ e1 e2))
(cdr (cdr list-of-base-and-exponents))))
(cons (car list-of-base-and-exponents)
(combine-same-base (cdr list-of-base-and-exponents)))))])))
(combine-same-base
(mergesort
(pollard-factorize n)
(lambda (x y)
(let ([id-or-first (lambda (x)
(if (number? x)
x
(first x)))])
(<= (id-or-first x) (id-or-first y)))))))
(define prime-divisors/exponents factorize)
(define (prime-divisors n)
(map car (prime-divisors/exponents n)))
(define (prime-exponents n)
(map cadr (factorize n)))
;;; ALGORITHM 9.22 p-adic Newton Iteration [MCA, p.264]
; INPUT phi in Z[y] (represented a normal function f : Z -> Z)
; p in Z, l>0,
; g0 in Z with phi(g)=0 mod p, phi'(go) invertible mod p
; and a modular inverse s0 of phi'(g0) mod p
; OUTPUT
; g in R with phi(g)=0 mod p^l and g=g0 mod p
(define (p-adic-newton-iteration phi Dphi p l g0 s0)
(let ([r (inexact->exact (ceiling (integer-log2 l)))])
(let loop ([i 1]
[gi g0]
[si s0])
(cond
[(< i r) (let ([g_i+1 (modulo (- gi (* (phi gi) si)) (expt p (expt 2 i)))])
(loop (+ i 1)
g_i+1
(modulo (- (* 2 si) (* (Dphi g_i+1) si si)) (expt p (expt 2 i)))))]
[else (modulo (- gi (* (phi gi) si)) (expt p l))]))))
;(= (p-adic-newton-iteration (lambda (y) (- (* y y y y) 1))
; (lambda (y) (* 4 (* y y y)))
; 5
; 4
; 2
; 3)
; 182)
; Return candidate if it's the nth root of a, otherwise #f
(define (is-nth-root a n candidate)
(if (= (expt candidate n) a)
candidate
#f))
(define (integer-root/odd-odd a n)
; INPUT a odd, n odd
; OUTPUT The n'th root of a, if it's an integer, #f otherwise
(unless (and (odd? a) (odd? n))
(error "integer-root/odd-odd: Both a and n must be odd; given " a n))
; Newton iteration with phi(y)=y^n-a and initial guess g0=1
(let ([candidate
; Newton iteration with phi(y)=y^n-a and initial guess g0=1
(let* ([k (do ([k 1 (add1 k)])
[(> (expt 2 (* n k)) a) k])]
[r (inexact->exact (ceiling (integer-log2 k)))])
(let loop ([i 1]
[gi 1]
[si 1]
[ti 1])
; (display `((k ,k) (r ,r) (i ,i) (gi ,gi) (si ,si) (ti ,ti))) (newline)
(cond
[(< i r) (let* ([g_i+1 (modulo (- gi (* (- (* gi ti) a) si))
(expt 2 (expt 2 i)))]
[t_i+1 (modulo (expt g_i+1 (- n 1))
(expt 2 (expt 2 (+ i 1))))])
(loop (+ i 1)
g_i+1
(modulo (- (* 2 si) (* n t_i+1 si si))
(expt 2 (expt 2 i)))
t_i+1))]
[else (modulo (- gi (* (- (* gi ti) a) si))
(expt 2 (expt 2 i)))])))])
(is-nth-root a n candidate)))
(define (integer-root/power-of-two a n)
; INPUT n a power of 2
; gcd(6,a)=1
; OUTPUT
;
(let ([phi (lambda (y) (- (expt y n) a))]
[Dphi (lambda (y) (* n (expt y (- n 1))))])
(let ([candidate1 (p-adic-newton-iteration phi Dphi 3 11 1 (inverse (Dphi 1) 3))])
(if (= (expt candidate1 n) a)
candidate1
(let ([candidate2 (p-adic-newton-iteration phi Dphi 3 11 2 (inverse (Dphi 2) 3))])
(is-nth-root a n candidate2))))))
(define (integer-root-factor a n)
; factor a = 2^d 3^e b^r , where gcd(6,b)=1
(let* ([d (max-dividing-power 2 a)]
[e (max-dividing-power 3 a)]
[b-to-r (quotient a (* (expt 2 d) (expt 3 e)))]
; factor n = 2^f c , where gcd(2,c)=1
[f (max-dividing-power 2 n)]
[two-to-f (expt 2 f)]
[c (quotient n two-to-f)]
;
[b (integer-root/power-of-two
(integer-root/odd-odd b-to-r c)
two-to-f)]
[r (max-dividing-power b b-to-r)])
(list d e b r)))
(define (integer-root a n)
; factor a = 2^d 3^e b^r , where gcd(6,b)=1
(if (= n 1)
a
(let ([d (max-dividing-power 2 a)])
(if (not (divides? n d))
#f
(let ([e (max-dividing-power 3 a)])
(if (not (divides? n e))
#f
(let* ([b-to-r (quotient a (* (expt 2 d) (expt 3 e)))]
; factor n = 2^f c , where gcd(2,c)=1
[f (max-dividing-power 2 n)]
[two-to-f (expt 2 f)]
[c (quotient n two-to-f)])
;
(cond
[(integer-root/odd-odd b-to-r c)
=> (lambda (cth-root--of--b-to-r)
(cond
[(integer-root/power-of-two cth-root--of--b-to-r two-to-f)
=> (lambda (nth-root--of--b-to-r)
(* (expt 2 (quotient d n))
(expt 3 (quotient e n))
nth-root--of--b-to-r))]
[else #f]))]
[else #f]))))))))
; INPUT a>0
; OUTPUT (values b r) where b^r=a and r maximal
(define (as-power a)
(let loop ([n (+ (big-log2 a) 1)])
;(display n) (newline)
(let ([root (integer-root a (add1 n))])
;(display "> ") (display root) (newline)
(if root
(values root (add1 n))
(loop (sub1 n))))))
(define (perfect-power? a)
(let-values ([(base n) (as-power a)])
(and (> n 1) (> a 1))))
(define (perfect-power a)
(let-values ([(base n) (as-power a)])
(if (and (> n 1) (> a 1))
(list base n)
#f)))
![[logo]](/logo.png){kind=logo}