;;; number-theory.scm  --  Jens Axel Søgaard

(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)
  ;  if p is prime.
  
  ; 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))
  
  #;(define (max-dividing-power p n)
      (define (loop n-divided-by-p-to-e e)
        (if (divides? p n-divided-by-p-to-e)
            (loop (quotient n-divided-by-p-to-e p) (+ e 1))
            e))
      (if (= p 1)
          n
          (loop n 0)))
  
  (define (max-dividing-power2 p n)
    (define (find-start p-to-e e)
      ;(display (list 'fs 'p-to-e p-to-e  'e e)) (newline)
      ; p-to-e divides n  and  p-to-e = p^e
      (let ([p-to-e2 (square p-to-e)])
        (cond [(= p-to-e2 n) (* 2 e)]
              [(> p-to-e2 n) (find-power p-to-e e)]
              [(divides? p-to-e2 n) (if (divides? p (quotient n p-to-e2))
                                        (find-start p-to-e2 (* 2 e))
                                        (* 2 e))]
              [else (find-power p-to-e e)])))
    (define (find-power p-to-e e)
      ;(display (list 'fp 'p-to-e p-to-e  'e e)) (newline)
      ; p-to-e <= n < (square p-to-e)
      (+ e (max-dividing-power p (quotient n p-to-e))))
    (cond [(= p 1)              1]
          [(not (divides? p n)) 0]
          [else                 (find-start p 1)]))
  
  
  ; (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 (phi2 n)
    (* n (product
          (map (lambda (p) (- 1 (/ 1 p)))
               (prime-divisors n)))))
  
  (define (phi n)
    (let ((divisors (prime-divisors n)))
      (* (quotient n (product divisors))
         (product (map (lambda (p) (- p 1)) divisors)))))
  
  ; 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)
    ; TODO TODO
    (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 (and (< i ny) (not (= b 1)))
                                       (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 (integer-length (/ 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)))]
          [(simple-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 (integer-length 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 (integer-length 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/power-of-two a n)
    ; INPUT    n = 2^d
    ; OUTPUT   an n'th root of a, or #f
    (let loop ([d (- (integer-length n) 1)]
               [b a])
      (if (= d 0)
          b
          (let-values ([(s r) (integer-sqrt/remainder b)])
            (if (not (zero? r))
                #f
                (loop (- d 1) s))))))
  
  
  (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
      (cond
        [(= n 1) a]
        [(= n 2) (let-values ([(s r) (integer-sqrt/remainder a)])
                   (if (zero? r)
                       s
                       #f))]
        [else
         (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]))))))]))
  
  (define (integer-root/remainder a n)
    (let ([i (integer-root a n)])
      (values i (- a (expt i n)))))
  
  (define (integer-root x y)
    ; y'th root of x
    (cond 
      ((eq? x 0)
       0)
      ((eq? x 1)
       1)
      ((eq? y 1)
       x)
      ((eq? y 2)
       (integer-sqrt x))
      ((not (integer? y))
       1)
      (else
       (let ((length (integer-length x)))
         ;; (expt 2 (- length l 1)) <= x < (expt 2 length)
         (cond ((<= length y)
                1)
               ;; result is >= 2
               ((<= length (* 2 y))
                ;; result is < 4
                (if (< x (expt 3 y))
                    2
                    3))
               ((even? y)
                (integer-root (integer-sqrt x)
                              (quotient y 2)))
               (else
                (let* ((length/y/2 ;; length/y/2 >= 1 because (< (* 2 y) length)
                        (quotient (quotient (- length 1) y) 2)))
                  (let ((init-g
                         (let* ((top-bits
                                 (arithmetic-shift
                                  x
                                  (- (* length/y/2 y))))
                                (nth-root-top-bits
                                 (integer-root top-bits y)))
                           (arithmetic-shift
                            (+ nth-root-top-bits 1)
                            length/y/2))))
                    (let loop ((g init-g))
                      (let* ((a (expt g (- y 1)))
                             (b (* a y))
                             (c (* a (- y 1)))
                             (d (quotient (+ x (* g c)) b)))
                        (let ((diff (- d g)))
                          (cond ((not (negative? diff))
                                 g)
                                ((< diff -1)
                                 (loop d))
                                (else
                                 ;; once the difference is one, it's more
                                 ;; efficient to just decrement until g^y <= x
                                 (let loop ((g d))
                                   (if (not (< x (expt g y)))
                                       g
                                       (loop (- g 1)))))))))))))))))
  
  ; INPUT    a>0
  ; OUTPUT   (values b r)  where b^r=a and r maximal
  (define (simple-as-power a)
    ; Note: The simple version is used by pollard-factorize
    (let loop ([n (integer-length a)])
      (let-values ([(root rem) (integer-root/remainder a (add1 n))])
        (if (zero? rem)
            (values root (add1 n))
            (loop (sub1 n))))))
  
  (define (as-power a)
    (let ([r (apply gcd (map second (factorize a)))])
      (values (integer-root a r) r)))
  
  #;(define (prime-power? n)
      (let-values ([(b r) (as-power n)])
        (prime? b)))
  
  (define (perfect-power? a)
    (let-values ([(base n) (as-power a)])
      (and (> n 1) (> a 1))))
  
  (define (simple-perfect-power a)
    ; simple-perfect-power is used by pollard-fatorize
    (let-values ([(base n) (simple-as-power a)])
      (if (and (> n 1) (> a 1))
          (list base n)
          #f)))
  
  (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))))]))
  
  (define divisor-sum 
    (let ()
      (case-lambda 
        [(n)   (divisor-sum n 1)]
        [(n k) (let* ([f  (factorize n)]
                     [ps (map first f)]
                     [es (map second f)])
                 (define (divisor-sum0 p e)
                   (+ e 1))
                 (define (divisor-sum1 p e)
                   (let loop ([sum 1]
                              [n 0]
                              [p-to-n 1])
                     (cond [(= n e) sum]
                           [else (let ([t (* p p-to-n)])
                                   (loop (+ t sum) (+ n 1) t))])))
                 (define (divisor-sumk p e)
                   (let ([p-to-k (expt p k)])
                     (let loop ([sum 1]
                                [n 0]
                                [p-to-kn 1])
                       (cond [(= n e) sum]
                             [else (let ([t (* p-to-k p-to-kn)])
                                     (loop (+ t sum) (+ n 1) t))]))))
                 (product
                  (map (cond [(= k 0) divisor-sum0]
                             [(= k 1) divisor-sum1]
                             [else    divisor-sumk])
                       ps es)))])))
  
  (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 (:range 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))))
                                     <))])))
  
  
  ; unreasonable long time for integer-root:
  ; (x 40398581 y 8)
  
  )