[plt-scheme] Fast primality testing

From: Jens Axel Søgaard (jensaxel at soegaard.net)
Date: Thu Jul 7 03:13:09 EDT 2005

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.