Initial population.

1 files changed, 736 insertions, 0 deletions

diff --git a/src/algebra/numtheor.spad.pamphlet b/src/algebra/numtheor.spad.pamphlet
new file mode 100644
index 00000000..b5cf7404
--- /dev/null
+++ b/src/algebra/numtheor.spad.pamphlet
@@ -0,0 +1,736 @@
+\documentclass{article}
+\usepackage{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
+    for k in 0.. while even? b repeat b := b quo 2
+    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
+        for k in 0.. until odd? a repeat a := a quo 2
+        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 := concat(1,nil())
+    for f in factors factor n repeat
+      newList : List Integer := nil()
+      for k in 0..f.exponent repeat
+        pow := f.factor ** k
+        for m in oldList repeat
+          newList := concat(pow * m,newList)
+      oldList := newList
+    sort(#1 < #2,newList)
+
+  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}