MATH Jens Axel Søgaard Version 1.0 5th Aug 2007

--------------------------------------------------------------------
  _MATH_   Jens Axel Søgaard    Version 1.0  5th Aug 2007
--------------------------------------------------------------------

(require (planet "math.ss" ("soegaard" "math.plt" 1 1)))

This library contains number theory code that I have written over a
long period. The code began as an experiment. I grabbed a book on
number theory from the shelve ("Elementary Number Theory" by Gareth
A. Jones and J. Mary Jones) and began illustrating each definition and
each theorem with Scheme code. The first half of the surce code is
thus a well commented mix of definitions, theorems and code.

The second half contains more sophisticated algorithms mostly of from
the excellent book "Modern Computer Algebra" by Joachim von zur Gathen
and Jürgen Gerhard. The algorithms for factorizing large intgers come
from this book. 

Finally there are some definitions of special functions, mostly
inspired by the problems of the Euler Project.

Please write with bugs, suggestions, and improvements.


HISTORY
-------

Version 1.1 [7th Aug 2007]
 - ported Bradley Lucier's integer-root from Gambit
 - better algorithm for max-dividing-power
 - better algorithm for divisor-sum
 - fixed some documentation bugs 
 Thanks to Bradley Lucier and David Einsteing
 for code and comments.

Version 1.0 [5th Aug 2007]
 - first PLaneT release
 

_NUMBER THEORY_   
----------------

There are quite a number theoretic functions available.
A few highlights:

Prime tests and factorization works for quite big numbers.

  > (prime? 567348343)
  #f

  > (factorize 5673483433473648736487364)
  ((2 2) (3 2) (31 1) (4609093 1) (24487651 1) (45042553 1))

For modular arithmetic the form with-modulus is provided. It
evaluates it's body in an environment where + - * and ^ work
modulo the given integer. Furhermore inv computes the inverse.

  ; Euler disproved Fermat's assertion that all Fermat numbers were prime
  > (with-modulus 641
      (+ (^ 2 (^ 2 5)) 1))
  0

  ; The eleventh Fermat number is also composite.
  > (with-modulus 319489
      (+ (^ 2 (^ 2 11)) 1))
  0

  Enter (+ (^ 2 (^ 2 11)) 1) in the repl to see the number of digits.

The function bezout computes the coeffecient from the linear
combination for the greatest common divisor. 

 (bezout a b) = (list u v)   <=>   gcd(a,b) = au + bv

It also works for more numbers than two.

  > (bezout 2 3)
  (-1 1)

  > (bezout 6 15 35)
  (1 2 -1)
  > (+ (* 1 6) (* 2 15) (* -1 35))
  1
  > (gcd 6 5 35)
  1

The function solve-chinese solves the equations in the Chinese
Remainder Theorem.

      (solve-chinese (a1 ... ak) (n1 ... nk)) = x
  <=> x = a1 mod n1, ..., x=ak mod nk

  >(solve-chinese '(2 3 2) '(3 5 7))
  23


REFERENCE (In alphabetical order)
--------------------------------


_^_ 
  Alias for expt.

_divides?_ : integer integer -> integer
  (divides a b) <=> a divides b

_interval_ : integer integer -> integer
  (interval m n) = (list m m+1 ... n)

_as-power_ : natural -> (values integer integer)
  For an a>0 there exists an maximal r such that
    a = b^r  for an integer b
  (as-power a) return the values b and r

_bernoulli_ : natural -> fraction
  (bernoulli n) returns the n'th Bernoulli number
  See: <http://mathworld.wolfram.com/BernoulliNumber.html>
       <http://en.wikipedia.org/wiki/Bernoulli_number>
  Note: The implementation uses Ramajuan's improvement
        of the standard recurrence relation.

_bezout-binary_ : integer integer -> (list integer integer)
  If a and b are non-zero integers then, 
  there exists integers u and v such that
     gcd(a,b) = au + bv
  Bezout returns (list u v). 
  Note: u and v are not unique.

_bezout_ : integer ... -> (list integer ...)
       (bezout-binary a b c ...) = (list u v w ...)    
  <=>  gcd(a,b,c,...) = au + bv + cw + ...

_binomial_ : natural natural -> natural
  (binomial n m) returns the binomial coeffecient n choose k
 
_coprime?_ integer ... -> boolean
      (coprime a b ...) 
  <=> (gcd a b ...) = 1

_digits_ : integer -> (list digit ...)
  (digits n) returns a list of the digits of n
  Example: (digits 123) = (list 1 2 3).

_digits->number_ : (list digit ...) -> number
  (digits->number (list d1 ...)) returns
  the number whose digits are d1,...

_diviors_ : naural -> (list natural ...)
  (diviors n) = (list d1 ... dk),
  where d1=1, ... dk=n is all positive divisors of n.

_divisor-sum_ : natural [integer] -> number
  (divisor-sum n) returns the sum of all divisors of n
  (divosr-sum n k) returns the sum of all kth powers of divisors of n

_euler-number_ : natural natural -> natural
  (euler-number n m) calculates the Euler number <n,k>
  Note: The implementation uses the standard recurrence
        <n,k> = (k+1) <n-1,k> + (n-k) <n-1,k-1>

_euler-phi_ : natural -> natural
  See phi.

_exists-primitive-root?_ : natural -> boolean
  Theorem      Un is cyclic (i.e. a primitive root exists)
           <=> n = 1, 2, 4, p^e or 2*p^e where p is an odd prime
  (exists-primitive-root? n) returns true, 
  if a pritive root exists in Un. 

_factorial_ : natural -> natural
  (factorial n) = (* 1 2 ... n)
  Note: For large n the implementation uses product splitting.

_factorize_ : natural -> (list (list prime natural) ...)
  The Fundamental Theorem of Arithmetic states
  that n>1 can be written uniquly as
    n = p1^e1 ... pk^ek where ei>0 and pi is prime

  (factorize n) returns (list (list p1 e1) ... (list pk ek))
  Note: This function is also called factorize.

_farey_ : natural -> (list fraction ...)
  (farey n) returns the sorted list of all fraction 
  in the interval 0 to 1, whose numerator is smaller or equal to n.
  Example: (farey 3) = (list 0 1/3 1/2 2/3 1)

_fermat_ : natural -> integer
  (fermat n) = Fn, where Fn is the n'th Fermat number:
      Fn = 2^(2^n) + 1
   I.e. F1=5, F2=17, F3=257, ...

_fibonacci_ : natural -> natural
  Returns the n'th Fibonacci number.

_gcd_ : integer integer ... -> integer
  (gcd a b ...) returns a positive greatest divisor of a, b, ...

_heptagonal_ : natural -> natural
  Returns the n'th heptagonal number

_heptagonal?_ : natural -> boolean
  (hexagonal? n) <=> n is a heptagonal number

_hexagonal_ : natural -> natural
  Returns the n'th hexagonal number

_hexagonal?_ : natural -> boolean
  (hexagonal? n) <=> n is a hexagonal number

_inverse_ : integer natural -> integer
  If gcd(a,n)=1 then there exists b such that
    ab=1 mod n
  The number b is called the inverse of a modulo n.
  (inverse a n) returns b such that ab=1 mod n and b in {0,...,n-1}

_legendre_ : integer integer -> integer
  The Legendre symbol of a and p is defined as:
                    /  1 if gcd(a,p)=1 and a is a square modulo p
  (legendre a p) = {  -1 if gcd(a,p)=1 and a is not a square modulo p
                    \  0 if gcd(a,p)<> 1
  Note: p must be prime
  Note: See also quadractic-residue?

_lucas_ : natural -> natural
  Returns the n'th Lucas number

_max-dividing-power_ : natural integer -> natural
      (max-dividing-power p n) = m  
  <=> p^m divides n  and  p^(m+1) doesn't divide n
  Note: In Mathematica this one is called IntegerExponent

_mediant_ : fraction fraction -> number
  If x=a/b and y=c/d then the mediant of x and y is (a+c)/(b+d).

_mersenne_ : natural -> natural
  (mersenne n) = 2^n-1
  Note: If p is a prime, then 2^p-1 is a Mersenne number

_moebius-mu_ : natural -> {-1,0,1}
                   / 1 if n is a product of an even number of primes
  moebius-mu(n) = { -1 if n is a product of an odd number of primes
                   \ 0 if n has multiple prime factors

_natural?_ : number -> boolean
  (natural? n) = (and (integer? n) (positive? n))

_next-prime_ : integer -> prime
  (next-prime n) returns the smallest prime larger than n.
  Note: See also prev-prime.

_next-primes_ : integer natural -> (list prime ...)
  (next-primes n k) returns a list of the k smallest primes greater
  than n.

_number-append_ : natural natural -> natural
  (number-append n1 n2) returns the number whose
  digits are the digits from n1 and n2.
  Example: (number-append 123 456) = 123456

_number-length_ : natural -> natural
  Returns the number of digits.

_number-reverse_ : natural -> natural
  (number-reverse n) returns the number r,
  whose digts are the same as n, but in reversed order.
  Example: (number-reverse 123) = 321

_nth-prime_ : natural -> prime
  Let p1=2, p2=3, p3=5, ... be the series of primes,
  then (nth-prime n) returns pn.

_octagonal_ : natural -> natural
  Returns the n'th octagonal number

_octagonal?_ : natural -> boolean
  (hexagonal? n) <=> n is a octagonal number

_odd-prime?_ integer -> boolean
      (odd-prime? n)
  <=> (and (odd? n) (prime? n))

_odd-prime-power?_ : natural -> boolean
      (odd-prime-power n)
  <=> (and (odd n) (prime-power? n))

_order_ : natural natural -> natural
  If G is a finite group with identity e, then the order of g in G is
  the least k>0 such that g^k=e. 
  (order g n) returns the order of g in Un, where Un is the group of
  units in Zn with respect to multiplication modulo n. 
  Note: g and n must be coprime.

_orders_ : natural -> (list natural)
  (orders n) = (list (order g1 n) ... (order gk n),
  where g1,...gk is the elements in (unit-group n).

_pairwise-coprime?_ :  integer ... -> boolean
      (pairwise-coprime? x1 ...)
  <=> for all i,j with i<>j :  (gcd xi xj)=1

_palindromic?_ : natural -> boolean
  A number n is palindromic if 
   (digits n) and (reverse (digits n)) are equal

_partitions_ : natural -> natural
  (parititions n) returns the number of partitions of n.

_pentagonal_ : natural -> natural
  Returns the n't pentagonal number.

_pentagonal?_ : natural -> boolean
  (pentagonal? n) <=> n is a pentagonal number

_perfect-power_ : natural -> (list base exponent) or #f
  If a can be written as a=b^r with r>1 maximal,
  then (perfect-power a) returns (list b r).
  Otherwise #f is returned.

_perfect-power?_ : natural -> boolean
      (perfect-power a) 
  <=> a can be written as a power b^r with r>1

_perfect-square_ : natural -> natural or #f
  If n=s^2 for an integer s, the (perfect-square n) 
  returns s, otherwise #f is returned.

_phi_ : natural -> natural
  Euler's phi function of 
             e1     ek
       n = p1 ... pk     , where pi is prime and ei>0
  is
                     k          1
     phi(n) = n * product (1 - ---- )
                    i=1         pi
  Note: This function is also called totient.

_powers-of_ : number natural -> (list integer ...)
  (powers-of a n) = (list a^0 a^1 ... a^n)

_prev-prime_ : integer -> #f or prime
  (prev-prime n) returns the greatest prime smaller than n.
  Unless no such prime exists, in which case #f is returned.

_next-primes_integer natural -> (list prime ...)
  (prev-primes n k) returns a list of the k greatest primes smaller
  than n.

_prime?_ : natural -> boolean  
      (prime? n) where n>1
  <=> 1 and n is the only positive divisors of n

_prime-diviors_ : natural -> (list prime ...)
  The Fundamental Theorem of Arithmetic states
  that n>1 can be written uniquly as
    n = p1^e1 ... pk^ek where ei>0 and pi is prime

  (prime-divisors n) returns (list p1 ... pk),
  that is, it returns the list of all primes dividing n

_prime-divisors/exponents_ : natural -> (list (list natural natural) ...)
  The Fundamental Theorem of Arithmetic states
  that n>1 can be written uniquly as
    n = p1^e1 ... pk^ek where ei>0 and pi is prime

  (prime-divisors/exponents n) returns (list (list p1 e1) ... (list pk ek))
  Note: This function is also called factorize.

_prime-exponents_ : natural -> (list natural ...)
  The Fundamental Theorem of Arithmetic states
  that n>1 can be written uniquly as
    n = p1^e1 ... pk^ek where ei>0 and pi is prime

  (prime-exponents n) returns (list e1 ... ek)

_prime-power_ : natural -> #f or (list prime exponent)
  If n=p^k for p prime and k>0, then
  (prime-power n) returns (list p k),
  otherwise #f is returned.

_prime-power?_ : natural -> boolean
      (prime-power? n)
  <=> n is a prime power

_primitive-root?_ : natural natural -> boolean
  A generator g of Un is called a primitive root mod n. Thus
       (primitive-root? g n)
  <=>  g is a primitive of Un
  <=>  g is a generator of Un 
  <=>  order(g) = phi(n)
  Note: The implementation uses the following observation:
         a in Un is a primitive root
    <=>   phi(n)/q
         a         <> 1  in Un for all primes q dividing phi(n)

_primitive-root_ : natural -> #f or natural
  (primitive-root n) returns a primitive root of Un, if
  such an root exists, otherwise #f is returned.

_primitive-roots_ : integer -> (list integer ...)
  (primitive-roots n) returns a list of all primitive roots of Un.
  Note: See primitive-root? for definition of a primitive root.

_positive-divisors_ : integer -> (list integer ...)
  (positive-divisors n) returns list of 
  all positive divisors of n 

_positive-divisors-of-n-below_ : integer integer -> (list integer ...)
  (positive-divisors n m) returns list of 
  all positive divisors of n less or equal to m

_product_ : (list number) -> number
  (product (list x1 ... xn)) = (* x1 ... xn)

_quadratic-residue?_ : natural natural -> boolean
  Definition:      a is a quadratic residue mod n
               <=> exists s such that a=s^2 (mod n)
  (quadratic-residue? a n) <=> a us a quadratic residue mod n

_quotient_ : integer integer -> integer
  If a and b in Z with b non-zero, then
  there exists q and r such that 
     a = qb+r   and 0<=r<=|b|
  q is returned by quotient

_random-integer_ : integer integer -> integer
  (random-integer from to) returns a random
  integer in the semi-open interval [from;to)

_remainder_ : integer integer -> integer
  If a and b in Z with b non-zero, then
  there exists q and r such that 
     a = qb+r   and 0<=r<=|b|
  r is returned by remainder

_second-degree-solutions_ : number number number -> (list number ...)
  (second-degree-solutions a b c) returns the list of real
  solutions to ax^2+bx+c=0.

_second-degree-integer-solutions_ : number number number -> (list integer ...)
  (second-degree-integer-solutions a b c) returns the list of integer
  solutions to ax^2+bx+c=0.

_second-degree-natural-solutions_ : number number number -> (list natural ...)
  (second-degree-natural-solutions a b c) returns the list of natural
  solutions to ax^2+bx+c=0.

_solve-chinese_ : (list integer ...) (list integer ...) -> integer
  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.
  (solve-chinese (list a1 ...) (list n1 ...) returns x.
  Example : (solve-chinese '(2 3 2) '(3 5 7)) = 23

_square_ : number -> number
  (square z) = z*z

_square?_ : natural -> boolean
  (square? n) <=> n is a perfect square

_sum_ : (list number) -> number
  (sum (list x1 ... xn)) = (+ x1 ... xn)

_tangent-number_ : natural -> natural
  (tangent-number n) returns the n'th tangent number
  See: http://mathworld.wolfram.com/TangentNumber.html
  Note: The implementation uses the method from 
  "Concrete Mathematics" p. 287

_totient_ : natural -> natural
  See phi.

_triangle_ : natural -> integer
  Returns the n't triangle number
  (triangle n) =  (quotient (* n (+ n 1)) 2)

_triangle?_ : natural -> booleand
  (triangle? n) <=> n is a triangle number

_unit-group_ : natural -> (list natural ...)
  The group of units in Zn with respect to multiplication modulo n is
  called Un. 
  (unit-group n) returns the list of units with respect to
  multiplication modulo n. This unit-group returns the elements of Un.

(with-modulus expr1 expr2 ...) : syntax
  The with-modulus form makes it easy to modular calculations. First
  the expresssion expr1 is evaluated, the result must be a natural
  number n. The body expressions expr2 ... is now evaluated in an
  environment, where the arithmetical operations +, -, * and ^
  (exponation) are automatically reduced modulo n.
  Example: (with-modulus 3 (^ 2 4)) evaluates to 1.