# [plt-scheme] Fast primality testing

 From: Jens Axel Søgaard (jensaxel at soegaard.net) Date: Thu Jul 7 03:13:09 EDT 2005 Previous message: [plt-scheme] fast primality testing Next message: [plt-scheme] Fast primality testing Messages sorted by: [date] [thread] [subject] [author]

```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)))

```

 Posted on the users mailing list. Previous message: [plt-scheme] fast primality testing Next message: [plt-scheme] Fast primality testing Messages sorted by: [date] [thread] [subject] [author]