path: root/src/algebra/numtheor.spad.pamphlet

\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra numtheor.spad}
\author{Martin Brock, Timothy Daly, Michael Monagan, Robert Sutor, 
Clifton J. Williamson}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package INTHEORY IntegerNumberTheoryFunctions}
\subsection{The inverse function}
The inverse function is derived from the 
{\bf Extended Euclidean Algorithm}.
If we divide one integer by another nonzero integer we get an integer
quotient plus a remainder which is, in general, a rational number. 
For instance, 
\[13/5 = 2 + 3/5\]
where 2 is the quotient and $3/5$ is the remainder.

If we multiply thru by the denominator of the remainder we get an answer
in integer terms which no longer involves division:
\[13 = 2(5) + 3\]

This gives a method of dividing integers. Specifically, if a and b are
positive integers, there exist unique non-negative integers q
and r so that
\[a = qb + r , {\rm\ where\ } 0 \le r < b\]
q is called the quotient and r the remainder.

The greatest common divisor of integers a and b, denoted by gcd(a,b),
is the largest integer that divides (without remainder) both a and
b. So, for example: gcd(15, 5) = 5, gcd(7, 9) = 1, gcd(12, 9) = 3,
gcd(81, 57) = 3.

The gcd of two integers can be found by repeated application of the
division algorithm, this is known as the Euclidean Algorithm. You
repeatedly divide the divisor by the remainder until the remainder is
0. The gcd is the last non-zero remainder in this algorithm. The
following example shows the algorithm.

Finding the gcd of 81 and 57 by the {\bf Euclidean Algorithm}:
\[
\begin{array}{rcl}
81 & = & 1(57) + 24\\
57 & = & 2(24) + 9\\
24 & = & 2(9) + 6\\
9 & = & 1(6) + 3\\
6 & = & 2(3) + 0
\end{array}
\]

So the greatest commmon divisor, the GCD(81,51)=3.

If the gcd(a, b) = r then there exist integers s and t so that
\[s(a) + t(b) = r\]

By back substitution in the steps in the Euclidean Algorithm, it is
possible to find these integers s and t. We shall do this with the
above example:

Starting with the next to last line, we have:
\[3 = 9 -1(6)\]
From the line before that, we see that $6 = 24 - 2(9)$, so:
\[3 = 9 - 1(24 - 2(9)) = 3(9) - 1(24)\]
From the line before that, we have $9 = 57 - 2(24)$, so:
\[3 = 3( 57 - 2(24)) - 1(24) = 3(57) - 7(24)\]
And, from the line before that $24 = 81 - 1(57)$, giving us:
\[3 = 3(57) - 7( 81 - 1(57)) = 10(57) -7(81)\]
So we have found s = -7 and t = 10.

The {\bf Extended Euclidean Algorithm} computes the GCD($a$,$b$) and 
the values for $s$ and $t$.

Suppose we were doing arithmetics modulo 26 and we needed to find the
inverse of a number mod 26. This turned out to be a difficult task
(and not always possible). We observed that a number $x$ had an inverse
mod 26 (i.e., a number $y$ so that $xy = 1 {\rm\ mod\ } 26$) 
if and only if $gcd(x,26) = 1$. 
In the general case the inverse
of $x$ exists if and only if $gcd(x, n) = 1$ and if it exists then
there exist integers $s$ and $t$ so that

\[sx + tn = 1\]

But this says that $sx = 1 + (-t)n$, or in other words, 
\[sx \equiv 1 {\rm\ mod\ } n\]  
So, s (reduced mod n if need be) is the inverse of $x {\rm\ mod\ }n$. 
The extended Euclidean algorithm calculates $s$ efficiently.

\subsubsection{Finding the inverse mod n}

We will number the steps of the Euclidean algorithm starting
with step 0. The quotient obtained at step $i$ will be denoted by $q_i$ 
and an auxillary number, $s_i$. For the first two steps, the value
of this number is given: $s_0 = 0$ and $s_1 = 1$. For the remainder of the
steps, we recursively calculate 
\[s_i = s_{i-2} - s_{i-1} q_{i-2} {\rm\ (mod\ n)}\] 

The algorithm starts by "dividing" $n$ by $x$. 
If the last non-zero remainder occurs at step $k$, 
then if this remainder is 1, $x$ has an inverse and it is $s_{k+2}$. 
If the remainder is not 1, then $x$ does not have an inverse. 

Find the inverse of 15 mod 26.
\[
\begin{array}{crcll}
Step 0: &26 &=& 1(15) + 11 &s_0 = 0\\
Step 1: &15 &=& 1(11) + 4  &s_1 = 1\\
Step 2: &11 &=& 2(4) + 3   
&s_2 = 0 - 1( 1) {\rm\ mod\ } 26 = 25\\
Step 3: &4  &=& 1(3) + 1   
&s_3 = 1 - 25( 1) {\rm\ mod\ } 26 = -24 {\rm\ mod\ } 26 = 2\\
Step 4: &3  &=& 3(1) + 0   
&s_4 = 25 - 2( 2) {\rm\ mod\ } 26 = 21\\
&&& &s_5 = 2 - 21( 1) {\rm\ mod\ } 26 = -19 {\rm\ mod\ } 26 = 7
\end{array}
\]
Notice that $15(7) = 105 = 1 + 4(26) \equiv 1 ({\rm\ mod\ } 26)$. 

Using the half extended Euclidean algorithm we compute 1/a mod b.
<<inverse(a,b)>>=
  inverse(a,b) ==
    borg:I:=b
    c1:I := 1
    d1:I := 0
    while b ~= 0 repeat
      q:I := a quo b
      r:I := a-q*b
      (a,b):=(b,r)
      (c1,d1):=(d1,c1-q*d1)
    a ~= 1 => error("moduli are not relatively prime")
    positiveRemainder(c1,borg)
  
@
<<inverse : (I,I) -> I>>=
  inverse : (I,I) -> I
@
Since this algorithm in local we need to reproduce it in the 
input file for testing purposes.
<<TESTinverse>>=
)clear completely 

inverse:(INT,INT)->INT

inverse(a,b) ==
  borg:INT:=b
  c1:INT := 1
  d1:INT := 0
  while b ~= 0 repeat
    q := a quo b
    r := a-q*b
    print [a, "=", q, "*(", b, ")+", r]
    (a,b):=(b,r)
    (c1,d1):=(d1,c1-q*d1)
  a ~= 1 => error("moduli are not relatively prime")
  positiveRemainder(c1,borg)

if ((inverse(26,15)*26)::IntegerMod(15) ~= 1) then print "DALY BUG"
if ((inverse(15,26)*15)::IntegerMod(26) ~= 1) then print "DALY BUG"

@
\subsection{The Chinese Remainder Algorithm}
\subsubsection{Chinese Remainder Theorem}
Let $m_1$,$m_2$,\ldots,$m_r$ be positive integers that are pairwise
relatively prime. Let $x_1$,$x_2$,\ldots,$x_r$ be integers with 
$0 \le x_i < m_i$. Then, there is exactly one $x$ in the interval
\[0 \le x < m_1 \cdot m_2 \cdots m_r\]
that satisfies the remainder equations
\[{\rm\ irem\ }(x,m_i) = x_i,\ \ \ i=1,2,\ldots,r\]
where {\bf irem} is the positive integer remainder function.
\subsubsection{Chinese Remainder Example}
Let $x_1 = 4$, $m_1 = 5$, $x_2 = 2$, $m_2 = 3$. We know that 
\[{\rm\ irem\ }(x,m_1) = x_1\]
\[{\rm\ irem\ }(x,m_2) = x_2\]
where $0 \le x_1 < m_1$ and $0 \le x_2 < m_2$. By the extended
Euclidean Algorithm there are integers $c$ and $d$ such that
\[c m_1 + d m_2 = 1\].
\noindent
In this case we are looking for an integer such that
\[{\rm\ irem\ }(x,5) = 4\]
\[{\rm\ irem\ }(x,3) = 2\]

The algorithm we use is to first 
compute the positive integer
remainder of $x_1$ and $m_1$ to get a new $x_1$:
\[
\begin{array}{rcl}
x_1 & = & {\rm\ positiveRemainder\ }(x_1,m_1)\\
4 & = & {\rm\ positiveRemainder\ }(4,5)
\end{array}
\]
Next compute the positive integer
remainder of $x_2$ and $m_2$ to get a new $x_2$:
\[
\begin{array}{rcl}
x_2 & = & {\rm\ positiveRemainder\ }(x_2,m_2)\\
2 & = & {\rm\ positiveRemainder\ }(2,3)
\end{array}
\]
Then we compute
\[x_1 + m_1 \cdot {\rm\ positiveRemainder\ }
(((x_2-x_1)\cdot{\rm inverse}(m_1,m_2)),m_2)\]
or
\[4+5*{\rm\ positiveRemainder\ }(((2-4)*{\rm\ inverse\ }(5,3)),3)\]
or
\[4+5*{\rm\ positiveRemainder\ }(-2*2),3)\]
or
\[4+5*2\]
or
\[14\]
<<chineseRemainder(x1,m1,x2,m2)>>=
  chineseRemainder(x1,m1,x2,m2) ==
    m1 < 0 or m2 < 0 => error "moduli must be positive"
    x1 := positiveRemainder(x1,m1)
    x2 := positiveRemainder(x2,m2)
    x1 + m1 * positiveRemainder(((x2-x1) * inverse(m1,m2)),m2)

@
This function has a restricted signature which only allows for
computing the chinese remainder of two numbers and two moduli.
<<chineseRemainder: (I,I,I,I) -> I>>=
  chineseRemainder: (I,I,I,I) -> I
    ++ \spad{chineseRemainder(x1,m1,x2,m2)} returns w, where w is such that
    ++ \spad{w = x1 mod m1} and \spad{w = x2 mod m2}. Note: \spad{m1} and
    ++ \spad{m2} must be relatively prime.
@
We test the particular example. The result should be 14.
<<TESTchineseRemainder>>=
)clear all
x1:=4
m1:=5
x2:=2
m2:=3
result:=chineseRemainder(x1,m1,x2,m2)
expected:=14
if (result-expected ~=0) then print "DALY BUG"

@
<<package INTHEORY IntegerNumberTheoryFunctions>>=
)abbrev package INTHEORY IntegerNumberTheoryFunctions
++ Author: Michael Monagan, Martin Brock, Robert Sutor, Timothy Daly
++ Date Created: June 1987
++ Date Last Updated:
++ Basic Operations:
++ Related Domains:
++ Also See:
++ AMS Classifications:
++ Keywords: number theory, integer
++ Examples:
++ References: Knuth, The Art of Computer Programming Vol.2
++ Description:
++ This package provides various number theoretic functions on the integers.
IntegerNumberTheoryFunctions(): Exports == Implementation where
 I ==> Integer
 RN ==> Fraction I
 SUP ==> SparseUnivariatePolynomial
 NNI ==> NonNegativeInteger

 Exports ==> with
  bernoulli : I -> RN
    ++ \spad{bernoulli(n)} returns the nth Bernoulli number.
    ++ this is \spad{B(n,0)}, where \spad{B(n,x)} is the \spad{n}th Bernoulli
    ++ polynomial.
<<chineseRemainder: (I,I,I,I) -> I>>
  divisors : I -> List I
    ++ \spad{divisors(n)} returns a list of the divisors of n.
  euler : I -> I
    ++ \spad{euler(n)} returns the \spad{n}th Euler number. This is
    ++ \spad{2^n E(n,1/2)}, where \spad{E(n,x)} is the nth Euler polynomial.
  eulerPhi : I -> I
    ++ \spad{eulerPhi(n)} returns the number of integers between 1 and n
    ++ (including 1) which are relatively prime to n. This is the Euler phi
    ++ function \spad{\phi(n)} is also called the totient function.
  fibonacci : I -> I
    ++ \spad{fibonacci(n)} returns the nth Fibonacci number. the Fibonacci
    ++ numbers \spad{F[n]} are defined by \spad{F[0] = F[1] = 1} and
    ++ \spad{F[n] = F[n-1] + F[n-2]}.
    ++ The algorithm has running time \spad{O(log(n)^3)}.
    ++ Reference: Knuth, The Art of Computer Programming
    ++ Vol 2, Semi-Numerical Algorithms.
  harmonic : I -> RN
    ++ \spad{harmonic(n)} returns the nth harmonic number. This is
    ++ \spad{H[n] = sum(1/k,k=1..n)}.
  jacobi : (I,I) -> I
    ++ \spad{jacobi(a,b)} returns the Jacobi symbol \spad{J(a/b)}.
    ++ When b is odd, \spad{J(a/b) = product(L(a/p) for p in factor b )}.
    ++ Note: by convention, 0 is returned if \spad{gcd(a,b) ~= 1}.
    ++ Iterative \spad{O(log(b)^2)} version coded by Michael Monagan June 1987.
  legendre : (I,I) -> I
    ++ \spad{legendre(a,p)} returns the Legendre symbol \spad{L(a/p)}.
    ++ \spad{L(a/p) = (-1)**((p-1)/2) mod p} (p prime), which is 0 if \spad{a}
    ++ is 0, 1 if \spad{a} is a quadratic residue \spad{mod p} and -1 otherwise.
    ++ Note: because the primality test is expensive,
    ++ if it is known that p is prime then use \spad{jacobi(a,p)}.
  moebiusMu : I -> I
    ++ \spad{moebiusMu(n)} returns the Moebius function \spad{mu(n)}.
    ++ \spad{mu(n)} is either -1,0 or 1 as follows:
    ++ \spad{mu(n) = 0} if n is divisible by a square > 1,
    ++ \spad{mu(n) = (-1)^k} if n is square-free and has k distinct
    ++ prime divisors.
  numberOfDivisors: I -> I
    ++ \spad{numberOfDivisors(n)} returns the number of integers between 1 and n
    ++ (inclusive) which divide n. The number of divisors of n is often
    ++ denoted by \spad{tau(n)}.
  sumOfDivisors : I -> I
    ++ \spad{sumOfDivisors(n)} returns the sum of the integers between 1 and n
    ++ (inclusive) which divide n. The sum of the divisors of n is often
    ++ denoted by \spad{sigma(n)}.
  sumOfKthPowerDivisors: (I,NNI) -> I
    ++ \spad{sumOfKthPowerDivisors(n,k)} returns the sum of the \spad{k}th
    ++ powers of the integers between 1 and n (inclusive) which divide n.
    ++ the sum of the \spad{k}th powers of the divisors of n is often denoted
    ++ by \spad{sigma_k(n)}.
 Implementation ==> add
  import IntegerPrimesPackage(I)

  -- we store the euler and bernoulli numbers computed so far in
  -- a Vector because they are computed from an n-term recurrence
  E: IndexedFlexibleArray(I,0)   := new(1, 1)
  B: IndexedFlexibleArray(RN,0)  := new(1, 1)
  H: Record(Hn:I,Hv:RN) := [1, 1]

  harmonic n ==
    s:I; h:RN
    n < 0 => error("harmonic not defined for negative integers")
    if n >= H.Hn then (s,h) := H else (s := 0; h := 0)
    for k in s+1..n repeat h := h + 1/k
    H.Hn := n
    H.Hv := h
    h

  fibonacci n ==
    n = 0 => 0
    n < 0 => (odd? n => 1; -1) * fibonacci(-n)
    f1, f2 : I
    (f1,f2) := (0,1)
    for k in length(n)-2 .. 0 by -1 repeat
      t := f2**2
      (f1,f2) := (t+f1**2,t+2*f1*f2)
      if bit?(n,k) then (f1,f2) := (f2,f1+f2)
    f2

  euler n ==
    n < 0 => error "euler not defined for negative integers"
    odd? n => 0
    l := (#E) :: I
    n < l => E(n)
    concat!(E, new((n+1-l)::NNI, 0)$IndexedFlexibleArray(I,0))
    for i in 1 .. l by 2 repeat E(i) := 0
    -- compute E(i) i = l+2,l+4,...,n given E(j) j = 0,2,...,i-2
    t,e : I
    for i in l+1 .. n by 2 repeat
      t := e := 1
      for j in 2 .. i-2 by 2 repeat
        t := (t*(i-j+1)*(i-j+2)) quo (j*(j-1))
        e := e + t*E(j)
      E(i) := -e
    E(n)

  bernoulli n ==
    n < 0 => error "bernoulli not defined for negative integers"
    odd? n =>
      n = 1 => -1/2
      0
    l := (#B) :: I
    n < l => B(n)
    concat!(B, new((n+1-l)::NNI, 0)$IndexedFlexibleArray(RN,0))
    for i in 1 .. l by 2 repeat B(i) := 0
    -- compute B(i) i = l+2,l+4,...,n given B(j) j = 0,2,...,i-2
    for i in l+1 .. n by 2 repeat
      t:I := 1
      b := (1-i)/2
      for j in 2 .. i-2 by 2 repeat
        t := (t*(i-j+2)*(i-j+3)) quo (j*(j-1))
        b := b + (t::RN) * B(j)
      B(i) := -b/((i+1)::RN)
    B(n)

<<inverse : (I,I) -> I>>
<<inverse(a,b)>>
<<chineseRemainder(x1,m1,x2,m2)>>

  jacobi(a,b) ==
    -- Revised by Clifton Williamson January 1989.
    -- Previous version returned incorrect answers when b was even.
    -- The formula J(a/b) = product ( L(a/p) for p in factor b) is only
    -- valid when b is odd (the Legendre symbol L(a/p) is not defined
    -- for p = 2).  When b is even, the Jacobi symbol J(a/b) is only
    -- defined for a = 0 or 1 (mod 4).  When a = 1 (mod 8),
    -- J(a/2) = +1 and when a = 5 (mod 8), we define J(a/2) = -1.
    -- Extending by multiplicativity, we have J(a/b) for even b and
    -- appropriate a.
    -- We also define J(a/1) = 1.
    -- The point of this is the following: if d is the discriminant of
    -- a quadratic field K and chi is the quadratic character for K,
    -- then J(d/n) = chi(n) for n > 0.
    -- Reference: Hecke, Vorlesungen ueber die Theorie der Algebraischen
    -- Zahlen.
    if b < 0 then b := -b
    b = 0 => error "second argument of jacobi may not be 0"
    b = 1 => 1
    even? b and positiveRemainder(a,4) > 1 =>
      error "J(a/b) not defined for b even and a = 2 or 3 (mod 4)"
    even? b and even? a => 0
    k : Integer := 0
    while even? b repeat
      b := b quo 2
      k := k + 1
    j:I := (odd? k and positiveRemainder(a,8) = 5 => -1; 1)
    b = 1 => j
    a := positiveRemainder(a,b)
    -- assertion: 0 < a < b and odd? b
    while a > 1 repeat
      if odd? a then
        -- J(a/b) = J(b/a) (-1) ** (a-1)/2 (b-1)/2
        if a rem 4 = 3 and b rem 4 = 3 then j := -j
        (a,b) := (b rem a,a)
      else
        -- J(2*a/b) = J(a/b) (-1) (b**2-1)/8
        k := 0
        until odd? a repeat
          a := a quo 2
          k := k + 1
        if odd? k and (b+2) rem 8 > 4 then j := -j
    a = 0 => 0
    j

  legendre(a,p) ==
    prime? p => jacobi(a,p)
    error "characteristic of legendre must be prime"

  eulerPhi n ==
    n = 0 => 0
    r : RN := 1
    for entry in factors factor n repeat
      r := ((entry.factor - 1) /$RN entry.factor) * r
    numer(n * r)

  divisors n ==
    oldList : List Integer := [1]
    for f in factors factor n repeat
      newList := oldList
      for k in 1..f.exponent repeat
        pow := f.factor ** k
        for m in oldList repeat
          newList := concat(pow * m,newList)
      oldList := newList
    sort(#1 < #2,oldList)

  numberOfDivisors n ==
    n = 0 => 0
    */[1+entry.exponent for entry in factors factor n]

  sumOfDivisors n ==
    n = 0 => 0
    r : RN := */[(entry.factor**(entry.exponent::NNI + 1)-1)/
      (entry.factor-1) for entry in factors factor n]
    numer r

  sumOfKthPowerDivisors(n,k) ==
    n = 0 => 0
    r : RN := */[(entry.factor**(k*entry.exponent::NNI+k)-1)/
      (entry.factor**k-1) for entry in factors factor n]
    numer r

  moebiusMu n ==
    n = 1 => 1
    t := factor n
    for k in factors t repeat
      k.exponent > 1 => return 0
    odd? numberOfFactors t => -1
    1

@
\subsection{TEST INTHEORY}
<<TEST INTHEORY>>=
<<TESTchineseRemainder>>
<<TESTinverse>>
@
\section{package PNTHEORY PolynomialNumberTheoryFunctions}
<<package PNTHEORY PolynomialNumberTheoryFunctions>>=
)abbrev package PNTHEORY PolynomialNumberTheoryFunctions
++ Author: Michael Monagan, Clifton J. Williamson
++ Date Created: June 1987
++ Date Last Updated: 10 November 1996 (Claude Quitte)
++ Basic Operations:
++ Related Domains:
++ Also See:
++ AMS Classifications:
++ Keywords: polynomial, number theory
++ Examples:
++ References: Knuth, The Art of Computer Programming Vol.2
++ Description:
++ This package provides various polynomial number theoretic functions
++ over the integers.
PolynomialNumberTheoryFunctions(): Exports == Implementation where
 I ==> Integer
 RN ==> Fraction I
 SUP ==> SparseUnivariatePolynomial
 NNI ==> NonNegativeInteger

 Exports ==> with
  bernoulli : I -> SUP RN
    ++ bernoulli(n) returns the nth Bernoulli polynomial \spad{B[n](x)}.
    ++ Note: Bernoulli polynomials denoted \spad{B(n,x)} computed by solving the
    ++ differential equation  \spad{differentiate(B(n,x),x) = n B(n-1,x)} where
    ++ \spad{B(0,x) = 1} and initial condition comes from \spad{B(n) = B(n,0)}.
  chebyshevT: I -> SUP I
    ++ chebyshevT(n) returns the nth Chebyshev polynomial \spad{T[n](x)}.
    ++ Note: Chebyshev polynomials of the first kind, denoted \spad{T[n](x)},
    ++ computed from the two term recurrence.  The generating function
    ++ \spad{(1-t*x)/(1-2*t*x+t**2) = sum(T[n](x)*t**n, n=0..infinity)}.
  chebyshevU: I -> SUP I
    ++ chebyshevU(n) returns the nth Chebyshev polynomial \spad{U[n](x)}.
    ++ Note: Chebyshev polynomials of the second kind, denoted \spad{U[n](x)},
    ++ computed from the two term recurrence.  The generating function
    ++ \spad{1/(1-2*t*x+t**2) = sum(T[n](x)*t**n, n=0..infinity)}.
  cyclotomic: I -> SUP I
    ++ cyclotomic(n) returns the nth cyclotomic polynomial \spad{phi[n](x)}.
    ++ Note: \spad{phi[n](x)} is the factor of \spad{x**n - 1} whose roots
    ++ are the primitive nth roots of unity.
  euler     : I -> SUP RN
    ++ euler(n) returns the nth Euler polynomial \spad{E[n](x)}.
    ++ Note: Euler polynomials denoted \spad{E(n,x)} computed by solving the
    ++ differential equation  \spad{differentiate(E(n,x),x) = n E(n-1,x)} where
    ++ \spad{E(0,x) = 1} and initial condition comes from \spad{E(n) = 2**n E(n,1/2)}.
  fixedDivisor: SUP I -> I
    ++ fixedDivisor(a) for \spad{a(x)} in \spad{Z[x]} is the largest integer
    ++ f such that f divides \spad{a(x=k)} for all integers k.
    ++ Note: fixed divisor of \spad{a} is
    ++ \spad{reduce(gcd,[a(x=k) for k in 0..degree(a)])}.
  hermite   : I -> SUP I
    ++ hermite(n) returns the nth Hermite polynomial \spad{H[n](x)}.
    ++ Note: Hermite polynomials, denoted \spad{H[n](x)}, are computed from
    ++ the two term recurrence.  The generating function is:
    ++ \spad{exp(2*t*x-t**2) = sum(H[n](x)*t**n/n!, n=0..infinity)}.
  laguerre  : I -> SUP I
    ++ laguerre(n) returns the nth Laguerre polynomial \spad{L[n](x)}.
    ++ Note: Laguerre polynomials, denoted \spad{L[n](x)}, are computed from
    ++ the two term recurrence.  The generating function is:
    ++ \spad{exp(x*t/(t-1))/(1-t) = sum(L[n](x)*t**n/n!, n=0..infinity)}.
  legendre  : I -> SUP RN
    ++ legendre(n) returns the nth Legendre polynomial \spad{P[n](x)}.
    ++ Note: Legendre polynomials, denoted \spad{P[n](x)}, are computed from
    ++ the two term recurrence.  The generating function is:
    ++ \spad{1/sqrt(1-2*t*x+t**2) = sum(P[n](x)*t**n, n=0..infinity)}.
 Implementation ==> add
  import IntegerPrimesPackage(I)

  x := monomial(1,1)$SUP(I)
  y := monomial(1,1)$SUP(RN)

  -- For functions computed via a fixed term recurrence we record
  -- previous values so that the next value can be computed directly

  E : Record(En:I, Ev:SUP(RN)) := [0,1]
  B : Record( Bn:I, Bv:SUP(RN) ) := [0,1]
  H : Record( Hn:I, H1:SUP(I), H2:SUP(I) ) := [0,1,x]
  L : Record( Ln:I, L1:SUP(I), L2:SUP(I) ) := [0,1,x]
  P : Record( Pn:I, P1:SUP(RN), P2:SUP(RN) ) := [0,1,y]
  CT : Record( Tn:I, T1:SUP(I), T2:SUP(I) ) := [0,1,x]
  U : Record( Un:I, U1:SUP(I), U2:SUP(I) ) := [0,1,0]

  MonicQuotient: (SUP(I),SUP(I)) -> SUP(I)
  MonicQuotient (a,b) ==
    leadingCoefficient(b) ~= 1 => error "divisor must be monic"
    b = 1 => a
    da := degree a
    db := degree b            -- assertion: degree b > 0
    q:SUP(I) := 0
    while da >= db repeat
      t := monomial(leadingCoefficient a, (da-db)::NNI)
      a := a - b * t
      q := q + t
      da := degree a
    q

  cyclotomic n ==
    --++ cyclotomic polynomial denoted phi[n](x)
    p:I; q:I; r:I; s:I; m:NNI; c:SUP(I); t:SUP(I)
    n < 0 => error "cyclotomic not defined for negative integers"
    n = 0 => x
    k := n; s := p := 1
    c := x - 1
    while k > 1 repeat
      p := nextPrime p
      (q,r) := divide(k, p)
      if r = 0 then
        while r = 0 repeat (k := q; (q,r) := divide(k,p))
        t := multiplyExponents(c,p::NNI)
        c := MonicQuotient(t,c)
        s := s * p
    m := (n quo s) :: NNI
    multiplyExponents(c,m)

  euler n ==
    p : SUP(RN); t : SUP(RN); c : RN; s : I
    n < 0 => error "euler not defined for negative integers"
    if n < E.En then (s,p) := (0$I,1$SUP(RN)) else (s,p) := E
    -- (s,p) := if n < E.En then (0,1) else E
    for i in s+1 .. n repeat
      t := (i::RN) * integrate p
      c := euler(i)$IntegerNumberTheoryFunctions / 2**(i::NNI) - t(1/2)
      p := t + c::SUP(RN)
    E.En := n
    E.Ev := p
    p

  bernoulli n ==
    p : SUP RN; t : SUP RN; c : RN; s : I
    n < 0 => error "bernoulli not defined for negative integers"
    if n < B.Bn then (s,p) := (0$I,1$SUP(RN)) else (s,p) := B
    -- (s,p) := if n < B.Bn then (0,1) else B
    for i in s+1 .. n repeat
      t := (i::RN) * integrate p
      c := bernoulli(i)$IntegerNumberTheoryFunctions
      p := t + c::SUP(RN)
    B.Bn := n
    B.Bv := p
    p

  fixedDivisor a ==
    g:I; d:NNI; SUP(I)
    d := degree a
    g := coefficient(a, minimumDegree a)
    for k in 1..d while g > 1 repeat g := gcd(g,a k)
    g

  hermite n ==
    s : I; p : SUP(I); q : SUP(I)
    n < 0 => error "hermite not defined for negative integers"
    -- (s,p,q) := if n < H.Hn then (0,1,x) else H
    if n < H.Hn then (s := 0; p := 1; q := x) else (s,p,q) := H
    for k in s+1 .. n repeat (p,q) := (2*x*p-2*(k-1)*q,p)
    H.Hn := n
    H.H1 := p
    H.H2 := q
    p

  legendre n ==
    s:I; t:I; p:SUP(RN); q:SUP(RN)
    n < 0 => error "legendre not defined for negative integers"
    -- (s,p,q) := if n < P.Pn then (0,1,y) else P
    if n < P.Pn then (s := 0; p := 1; q := y) else (s,p,q) := P
    for k in s+1 .. n repeat
      t := k-1
      (p,q) := ((k+t)$I/k*y*p - t/k*q,p)
    P.Pn := n
    P.P1 := p
    P.P2 := q
    p

  laguerre n ==
    k:I; s:I; t:I; p:SUP(I); q:SUP(I)
    n < 0 => error "laguerre not defined for negative integers"
    -- (s,p,q) := if n < L.Ln then (0,1,x) else L
    if n < L.Ln then (s := 0; p := 1; q := x) else (s,p,q) := L
    for k in s+1 .. n repeat
      t := k-1
      (p,q) := ((((k+t)$I)::SUP(I)-x)*p-t**2*q,p)
    L.Ln := n
    L.L1 := p
    L.L2 := q
    p

  chebyshevT n ==
    s : I; p : SUP(I); q : SUP(I)
    n < 0 => error "chebyshevT not defined for negative integers"
    -- (s,p,q) := if n < CT.Tn then (0,1,x) else CT
    if n < CT.Tn then (s := 0; p := 1; q := x) else (s,p,q) := CT
    for k in s+1 .. n repeat (p,q) := ((2*x*p - q),p)
    CT.Tn := n
    CT.T1 := p
    CT.T2 := q
    p

  chebyshevU n ==
    s : I; p : SUP(I); q : SUP(I)
    n < 0 => error "chebyshevU not defined for negative integers"
    if n < U.Un then (s := 0; p := 1; q := 0) else (s,p,q) := U
    for k in s+1 .. n repeat (p,q) := ((2*x*p - q),p)
    U.Un := n
    U.U1 := p
    U.U2 := q
    p

@
\section{License}
<<license>>=
--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
@
<<*>>=
<<license>>

<<package INTHEORY IntegerNumberTheoryFunctions>>
<<package PNTHEORY PolynomialNumberTheoryFunctions>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} Cohen, Joel S., 
{\sl Computer Algebra and Symbolic Computation}
{\sl Mathematical Methods},
A.K. Peters, Ltd, Natick, MA. USA (2003)
ISBN 1-56881-159-4
\bibitem{2} Geddes, Keith O., Czapor, Stephen R., Labahn, George
{\sl Algorithms for Computer Algebra}
Kluwer Academic Publishers
ISBN 0-7923-9259-0
\end{thebibliography}
\end{document}