number-theory.ss
```;;; number-theory.scm  --  Jens Axel Søgaard  --  5th aug 2007

(module number-theory mzscheme
(provide (all-defined))

(require (lib "include.ss"))
(require (lib "27.ss" "srfi"))
(require (only (lib "etc.ss") compose rec)
(only (lib "list.ss") filter mergesort) )
(define first car)
(define second cadr)
(define rest cdr)
(define empty '())

(define ^ expt)

; inteval : integer integer -> (list integer)
;   return (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))

;; product : (list number) -> number
;(define (product xs)
;  ; (prod (list x1 x2 x3 x4 ...))
;  ;  = (list (* x1 x2) (* x3 x4) ...)
;  (define (prod xs)
;    (cond
;      [(null? xs)        '()]
;      [(null? (cdr xs))  xs]
;      [else             (cons (* (car xs) (cadr xs))
;                              (prod (cddr xs)))]))
;  (if (null? xs)
;      1
;      (let loop ([xs xs])
;        (if (null? (cdr xs))
;            (car xs)
;            (loop (prod xs))))))

(define (product xs)
(apply * xs))

(define (sum xs)
(apply + xs))

(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 pairwise-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
;   In Mathematica this one is called IntegerExponent
(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 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  (product ns)]
[cs (map (lambda (ni) (/ n ni)) ns)]
[ds (map inverse cs ns)])
(modulo (sum (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 (product
(map (lambda (p)
(- 1 (/ 1 p)))
(prime-divisors n)))))

; Euler's phi function is also known as the totient.
(define totient phi)

;;; 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 there exists an s such that a=s^2 (mod n)
;   The number s os called a squre root of a modulo 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)))

(require (lib "42.ss" "srfi"))

(define *SMALL-PRIME-LIMIT* 1000000)

(define prime?
; TODO: make N settable
(let ()
(define N *SMALL-PRIME-LIMIT*)  ; NOTE: this affects
(define ps (make-vector (+ N 1) #t))
(define ! vector-set!)
(! ps 0 #f)
(! ps 1 #f)
(do-ec (: n 2 (+ N 1))
(if (vector-ref ps n))
(: m (+ n n) (+ N 1) n)
(! ps m #f))
(lambda (n)
(if (< n N)
(vector-ref ps n)
(prime-strong-pseudo? n)))))

(define (next-prime n)
(cond
[(< n 2)   2]
[(= 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]
[(< n 3)   #f]
[(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)])
(if next
(cons next (next-primes next (sub1 primes-wanted)))
'()))))

(define (prev-primes n primes-wanted)
(if (= primes-wanted 0)
'()
(let ([prev (prev-prime n)])
(if prev
(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)
(if (< n 10000)  ; todo: find best value
(factorize-small n)
(cond
[(= n 1)              '()]
[(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 (base-and-exponenent<? x y)
(let ([id-or-first
(lambda (x)
(if (number? x)
x
(first x)))])
(<= (id-or-first x) (id-or-first y))))

(define (factorize n)
(if (< n *SMALL-PRIME-LIMIT*)   ; NOTE: Do measurement of best cut
(factorize-small n)
(factorize-large n)))

(require (only (lib "1.ss" "srfi") find-tail))

(define (small-prime-factors-over n p)
; p prime
(cond
[(< n p)         '()]
[(= n p)         (list (list p 1))]
[(prime? n)      (list (list n 1))]
[(divides? p n)  (let ([m (max-dividing-power p n)])
(cons (list p m)
(small-prime-factors-over
(quotient n (expt p m))
(next-prime p))))]
[else            (small-prime-factors-over n (next-prime p))]))

(define (factorize-small n)
; fast for small n, but works correctly
; for large n too
(small-prime-factors-over n 2))

(define (combine-same-base list-of-base-and-exponents)
; list-of-base-and-exponents must be sorted
(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)))))])))

(define (factorize-large n)
(combine-same-base
(mergesort
(pollard-factorize n)
base-and-exponenent<?)))

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

; powers-of : number natural -> number
;   returns a list of numbers: a^0, ..., a^n
(define (powers-of a n)
(let loop ([i    0]
[a^i  1])
(if (<= i n)
(cons a^i (loop (+ i 1) (* a^i a)))
'())))

(define (factorization->divisors f)
(cond
[(null? f) '(1)]
[else      (let ([p (first (first f))]
[n (second (first f))]
[g (rest f)])
; f = p^n * g
(let ([divisors-of-g (factorization->divisors g)])
(apply append
(map (lambda (p^i)
(map (lambda (d) (* p^i d))
divisors-of-g))
(powers-of p n)))))]))

(define (divisors n)
(factorization->divisors (factorize n)))

;;; Integer Functions

(require (only (lib "1.ss" "srfi") every))

(define (one? n)
(= n 1))

(define euler-phi phi)

; moebius-mu : natural -> {-1,0-1}
;   mu(n) =  1  if n is a product of an even number of primes
;         = -1  if n is a product of an odd number of primes
;         =  0  if n has a multiple prime factor
(define (moebius-mu n)
(let* ([f          (factorize n)]
[primes     (map first f)]
[exponents  (map second f)])
(if (every one? exponents)
(if (even? (length primes))
1
-1)
0)))

; divisor-sum : natural integer -> number
;   returns the sum of the kth power of all divisors of n
(define divisor-sum
(case-lambda
[(n)   (divisor-sum n 1)]
[(n k) (let loop ([sum 0] [ds (divisors n)])
(if (null? ds)
sum
(loop (+ sum (expt (car ds) k))
(cdr ds))))]))

(require (lib "42.ss" "srfi"))

;;; Bernoulli Numbers

; bernoulli : natutal -> fraction
;   compute the n'th Bernoulli number
;   <http://mathworld.wolfram.com/BernoulliNumber.html>
;   <http://en.wikipedia.org/wiki/Bernoulli_number>
(define (bernoulli n)
; Implementation note:
;   - uses Ramanujan's improvement of the standard recurrence relation
;     of the Bernoulli numbers <http://en.wikipedia.org/wiki/Bernoulli_number>
;   - memoizes previous computations
;   - avoids an explicit call to compute the binomials
;   TODO: Memoize the computed numbers, also from  call to call
;         (simple change, the numbers are already memoized internally,
;          just switch to evectors)
(let ([b (make-vector (max (+ n 1) 2) #f)])
(vector-set! b 0 1)
(vector-set! b 1 -1/2)
(let ()
(define (next-binom old x k)
; calculate binom(x,k-6) from the old binom(x,k)
(* old
(/ (* k (- k 1) (- k 2) (- k 3) (- k 4) (- k 5) )
(* (- x (- k 1)) (- x (- k 2)) (- x (- k 3)) (- x (- k 4)) (- x (- k 5)) (- x (- k 6))))))
(define (A m M)
(if (< M 1)
0
(sum-ec (:parallel
(: j 1 (+ M 1))
(:do ([bin (binomial (+ m 3) (- m 6))]) #t ((next-binom bin (+ m 3) (- m (* 6 j))))))
(* bin
(bern (- m (* 6 j)))))))
(define (bern n)
(let ([bn (vector-ref b n)])
(if bn
bn
(let ([bn
(cond
[(odd? n)               0]
[(= (remainder n 6) 0)  (/   (-  (/ (+ n 3)  3)  (A n (/ n       6)))   (binomial (+ n 3) n))]
[(= (remainder n 6) 2)  (/   (-  (/ (+ n 3)  3)  (A n (/ (- n 2) 6)))   (binomial (+ n 3) n))]
[(= (remainder n 6) 4)  (/   (-  (/ (+ n 3) -6)  (A n (/ (- n 4) 6)))   (binomial (+ n 3) n))])])
(vector-set! b n bn)
bn))))
(bern n))))

; tangent-number : natural -> natural
;  compute the n'th tangent number
;  <http://mathworld.wolfram.com/TangentNumber.html>
(define (tangent-number n)
; Implementation note:
;   See "Concrete Mathematics" p 287 for the method
(let ([T (make-vector (+ n 2) 0)])
; T[x] = x
(vector-set! T 1 1)
(do-ec (: k (+ n 1))
(begin
; differentiate T[x]
(do-ec (: i (+ k 1))
(vector-set! T i
(* (+ i 1) (vector-ref T (+ i 1)))))
; multiply T[x] with 1+x^2
(do-ec (: i (+ k 1) 1 -1)
(vector-set! T i
(+ (vector-ref T i)
(vector-ref T (- i 2)))))))
(vector-ref T 0)))

; bernoulli/tangent : natural -> fraction
;   compute the n'th Bernoulli number
(define (bernoulli/tangent n)
; Implementation note:
;   The Bernoulli number is computed on the basis
;   of the corresponding tangent number. Since
;   the tangent number computation uses bignum
;   integers only, this may or may not be faster
;   for a one-of computation than the (bernoulli n)
;   function.
(cond
[(= n 0)   1]
[(= n 1)  -1/2]
[(= n 2)   1/6]
[(odd? n)  0]
[else      (let ([m (/ n 2)])
(* (tangent-number (- n 1))
(/ n
(let ([4^m (expt 4 m)])
(* 4^m (- 4^m 1))))
(if (even? m)
-1
1)))]))

; binomial : natural natural -> natural
;  compute the binomial coeffecient n choose k
(define (binomial n k)
; ;  <http://www.swox.com/gmp/manual/Binomial-Coefficients-Algorithm.html>
; TODO: Check range of n and k
(cond
[(= k 0)            1]
[(= k 1)            n]
[(= k 2)            (/ (* n (- n 1)) 2)]
[(> k (/ n 2))      (binomial n (- n k))]
[else               (* (+ n (- k) 1)
(product-ec (:range i 2 (+ k 1))
(/ (+ n (- k) i)
i)))]))

; factorial : natural -> natural
;  computes n * (n-1) * ... 1
(define (factorial n)
(define (fact-small n)
(if (= n 1)
1
(* n (factorial (- n 1)))))
(define (fact-big n)
; Implementation note:
;   uses product splitting
;   TODO: further improvements, collect twos
;         see: <http://www.swox.com/gmp/manual/Factorial-Algorithm.html#Factorial-Algorithm>
(let ([n/2 (quotient n 2)])
(product (append
(list-ec (:range i 1 n/2) i)
(reverse (list-ec (:range i n/2 (+ n 1)) i))))))
(cond
[(= n 0)   1]
[(< n 400) (fact-small n)]
[else      (fact-big n)]))

; euler-number : N0 N0 -> N0
;   computes the Euler number <n,k>
(define (euler-number n k)
; Implementation note:
;   Uses standard recurrence : <n,k> = (k+1) <n-1,k> + (n-k) <n-1,k-1>
;   Time: =(n^2)
(if (= k 0)
1
(let ([E (make-vector (max (+ k 1) (+ n 1)) 0)])
(vector-set! E 0 1) ; <0,0> = 1
(do-ec (: i 1 (+ n 1))
(do-ec (:range j (- i 1) 0 -1)
(vector-set! E j
(+ (* (+ j 1) (vector-ref E j))
(* (- i j) (vector-ref E (- j 1)))))))
(vector-ref E k))))

#;(list-ec (: n 10) (list-ec (: k 10) (euler-number n k)))

(define (perfect-square n)
(let ([sqrt-n (integer-sqrt n)])
(if (= (* sqrt-n sqrt-n) n)
sqrt-n
#f)))

(define (second-degree-solutions a b c)
; return list of solutions to a a x^2 + b x + c = 0
(let ([d (- (* b b) (* 4 a c))])
(cond
[(< d 0)
'()]
[(= d 0)
(list (/ b (* -2 a)))]
[else
(let ([sqrt-d (sqrt d)])
(list (/ (- (- b) sqrt-d) (* 2 a))
(/ (+ (- b) sqrt-d) (* 2 a))))])))

(define (second-degree-integer-solutions a b c)
; return list of integer solutions to a x^2 + b x + c = 0
(filter integer? (second-degree-solutions a b c)))

(define (natural? n)
(and (integer? n)
(positive? n)))

(define (second-degree-natural-solutions a b c)
; return list of integer solutions to a x^2 + b x + c = 0
(filter natural? (second-degree-solutions a b c)))

(define (triangle n)
(quotient (* n (+ n 1)) 2))

(define (triangle? n)
(not (null? (second-degree-natural-solutions 1/2 1/2 (- n)))))

(define (square? n)
(perfect-square n))

(define (pentagonal n)
(quotient (* n (- (* 3 n) 1)) 2))

(define (pentagonal? n)
(not (null? (second-degree-natural-solutions 3/2 -1/2 (- n)))))

(define (hexagonal n)
(* n (- (* 2 n) 1)))

(define (hexagonal? n)
(not (null? (second-degree-natural-solutions 2 -1 (- n)))))

(define (heptagonal n)
(quotient (* n (- (* 5 n) 3)) 2))

(define (heptagonal? n)
(not (null? (second-degree-natural-solutions 5/2 -3/2 (- n)))))

(define (octagonal n)
(* n (- (* 3 n) 2)))

(define (octagonal? n)
(not (null? (second-degree-natural-solutions 3 -2 (- n)))))

(require (planet "memoize.ss" ("dherman" "memoize.plt" 2 1)))

(define (partitions n)
(define p
(memo-lambda* (n)
(cond
[(< n 0) 0]
[(= n 0) 1]
[else   (sum-ec (:range k 1 (add1 (inexact->exact (floor (/ (+ 1.0 (sqrt (+ 1.0 (* 24.0 n)))) 6.0)))))
(* (if (odd? k) 1 -1)
(+ (p (- n (quotient (* k (- (* 3 k) 1)) 2)))
(p (- n (quotient (* k (+ (* 3 k) 1)) 2))))))])))
(p n))

(define (digits n)
(define (d x)
(if (< x 10)
(list x)
(cons (remainder x 10)
(d (quotient x 10)))))
(unless (integer? n)
(error 'digits "expected an integer, got: n"))
(reverse (d (if (negative? n) (- n) n))))

(define (digits->number ds)
(define (rev-digits->number ds)
(if (null? ds)
0
(+ (car ds)
(* 10 (rev-digits->number (cdr ds))))))
(rev-digits->number (reverse ds)))

(define (number-append x y)
(unless (and (natural? x) (natural? y))
(error 'number-append "expected two natural numbers, got: " x y))
(string->number
(string-append
(number->string x)
(number->string y))))

(define (number-reverse n)
(digits->number
(reverse
(digits n))))

(define (palindromic? n)
(let ([ds (digits n)])
(equal? ds (reverse ds))))

(define (number-length n)
(unless (and (integer? n) (not (negative? n)))
(error 'number-lengh "expected non-negative integer, got: " n))
(let ([len (inexact->exact (floor (/ (log n) (log 10))))])
(let loop ([len len])
(if (< (- (expt 10 len) 1) n (expt 10 (+ len 1)))
(+ len 1)
(loop (+ len 1))))))

;;; Code for Fibonacci by Ian Glover / SICP ex. 1.19.

(define (generator a b p q count)
(cond ((= count 0) b)
((even? count)
(generator a
b
(+ (* p p) (* q q))
(+ (* 2 p q) (* q q))
(quotient count 2)))
(else (generator (+ (* b q) (* a q) (* a p))
(+ (* b p) (* a q))
p
q
(- count 1)))))

(define (fibonacci n)
(generator 1 0 0 1 n))

(define (lucas n)
(generator 1 3 0 1 n))

(define (fibonacci-mod n mod)
(define (generator a b p q count)
(cond ((= count 0) b)
((even? count)
(generator a
b
(remainder (+ (* p p) (* q q)) mod)
(remainder (+ (* 2 p q) (* q q)) mod)
(quotient count 2)))
(else (generator (remainder (+ (* b q) (* a q) (* a p)) mod)
(remainder (+ (* b p) (* a q)) mod)
p
q
(- count 1)))))
(generator 1 0 0 1 n))

(require (planet "memoize.ss" ("dherman" "memoize.plt")))

(define (mediant x y)
(/ (+ (numerator x) (numerator y))
(+ (denominator x) (denominator y))))

(require (only (lib "list.ss") sort))
(define farey
(memo-lambda* (d)
(define (delete-last xs)
(reverse! (cdr (reverse xs))))
(define (successive? x y)
(= (- (* (denominator x) (numerator y))
(* (denominator y) (numerator x)))
1))
(cond
[(= d 1) (list 0 1)]
[else    (let ([fs (farey (- d 1))])
(sort (append fs
(filter {lambda (x) (<= (denominator x) d)}
(map mediant
(delete-last fs) (rest fs))))
<))])))

)```