Initial population.

1 files changed, 4986 insertions, 0 deletions

diff --git a/src/hyper/pages/ug08.ht b/src/hyper/pages/ug08.ht
new file mode 100644
index 00000000..f14e0289
--- /dev/null
+++ b/src/hyper/pages/ug08.ht
@@ -0,0 +1,4986 @@
+% Copyright The Numerical Algorithms Group Limited 1992-94. All rights reserved.
+% !! DO NOT MODIFY THIS FILE BY HAND !! Created by ht.awk.
+\texht{\setcounter{chapter}{7}}{} % Chapter 8
+
+%
+\newcommand{\ugProblemTitle}{Advanced Problem Solving}
+\newcommand{\ugProblemNumber}{8.}
+%
+% =====================================================================
+\begin{page}{ugProblemPage}{8. Advanced Problem Solving}
+% =====================================================================
+\beginscroll
+
+In this chapter we describe techniques useful in solving advanced problems
+with \Language{}.
+
+\beginmenu
+    \menudownlink{{8.1. Numeric Functions}}{ugProblemNumericPage}
+    \menudownlink{{8.2. Polynomial Factorization}}{ugProblemFactorPage}
+    \menudownlink{{8.3. Manipulating Symbolic Roots of a Polynomial}}{ugProblemSymRootPage}
+    \menudownlink{{8.4. Computation of Eigenvalues and Eigenvectors}}{ugProblemEigenPage}
+    \menudownlink{{8.5. Solution of Linear and Polynomial Equations}}{ugProblemLinPolEqnPage}
+    \menudownlink{{8.6. Limits}}{ugProblemLimitsPage}
+    \menudownlink{{8.7. Laplace Transforms}}{ugProblemLaplacePage}
+    \menudownlink{{8.8. Integration}}{ugProblemIntegrationPage}
+    \menudownlink{{8.9. Working with Power Series}}{ugProblemSeriesPage}
+    \menudownlink{{8.10. Solution of Differential Equations}}{ugProblemDEQPage}
+    \menudownlink{{8.11. Finite Fields}}{ugProblemFinitePage}
+    \menudownlink{{8.12. Primary Decomposition of Ideals}}{ugProblemIdealPage}
+    \menudownlink{{8.13. Computation of Galois Groups}}{ugProblemGaloisPage}
+    \menudownlink{{8.14. Non-Associative Algebras and Modelling Genetic Laws}}{ugProblemGeneticPage}
+\endmenu
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugProblemNumericTitle}{Numeric Functions}
+\newcommand{\ugProblemNumericNumber}{8.1.}
+%
+% =====================================================================
+\begin{page}{ugProblemNumericPage}{8.1. Numeric Functions}
+% =====================================================================
+\beginscroll
+%
+\Language{} provides two basic floating-point types: \axiomType{Float} and
+\axiomType{DoubleFloat}.  This section describes how to use numerical
+%-% \HDindex{function!numeric}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+operations defined on these types and the related complex types.
+%-% \HDindex{numeric operations}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+%
+As we mentioned in \downlink{``\ugIntroTitle''}{ugIntroPage} in Chapter \ugIntroNumber\ignore{ugIntro}, the \axiomType{Float} type is a software
+implementation of floating-point numbers in which the exponent and the
+%-% \HDindex{floating-point number}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+significand may have any number of digits.
+%-% \HDindex{number!floating-point}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+See \downlink{`Float'}{FloatXmpPage}\ignore{Float} for detailed information about this domain.
+The \axiomType{DoubleFloat} (see \downlink{`DoubleFloat'}{DoubleFloatXmpPage}\ignore{DoubleFloat}) is usually a hardware
+implementation of floating point numbers, corresponding to machine double
+precision.
+The types \axiomType{Complex Float} and \axiomType{Complex DoubleFloat} are
+%-% \HDindex{floating-point number!complex}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+the corresponding software implementations of complex floating-point numbers.
+%-% \HDindex{complex!floating-point number}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+In this section the term {\it floating-point type}  means any of these
+%-% \HDindex{number!complex floating-point}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+four types.
+%
+The floating-point types implement the basic elementary functions.
+These include (where \axiomSyntax{\$} means
+\axiomType{DoubleFloat},
+\axiomType{Float},
+\axiomType{Complex DoubleFloat}, or
+\axiomType{Complex Float}):
+
+\noindent
+\axiomFun{exp},  \axiomFun{log}: \axiom{\$ -> \$} \newline
+\axiomFun{sin},  \axiomFun{cos}, \axiomFun{tan}, \axiomFun{cot}, \axiomFun{sec}, \axiomFun{csc}: \axiom{\$ -> \$} \newline
+\axiomFun{sin},  \axiomFun{cos}, \axiomFun{tan}, \axiomFun{cot}, \axiomFun{sec}, \axiomFun{csc}: \axiom{\$ -> \$}  \newline
+\axiomFun{asin}, \axiomFun{acos}, \axiomFun{atan}, \axiomFun{acot}, \axiomFun{asec}, \axiomFun{acsc}: \axiom{\$ -> \$}  \newline
+\axiomFun{sinh},  \axiomFun{cosh}, \axiomFun{tanh}, \axiomFun{coth}, \axiomFun{sech}, \axiomFun{csch}: \axiom{\$ -> \$}  \newline
+\axiomFun{asinh}, \axiomFun{acosh}, \axiomFun{atanh}, \axiomFun{acoth}, \axiomFun{asech}, \axiomFun{acsch}: \axiom{\$ -> \$}  \newline
+\axiomFun{pi}: \axiom{() -> \$}  \newline
+\axiomFun{sqrt}: \axiom{\$ -> \$} \newline
+\axiomFun{nthRoot}: \axiom{(\$, Integer) -> \$}  \newline
+\axiomFunFrom{**}{Float}: \axiom{(\$, Fraction Integer) -> \$} \newline
+\axiomFunFrom{**}{Float}: \axiom{(\$,\$) -> \$}  \newline
+
+The handling of roots depends on whether the floating-point type
+%-% \HDindex{root!numeric approximation}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+is real or complex: for the real floating-point types,
+\axiomType{DoubleFloat} and \axiomType{Float}, if a real root exists
+the one with the same sign as the radicand is returned; for the
+complex floating-point types, the principal value is returned.
+%-% \HDindex{principal value}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+Also, for real floating-point types the inverse functions
+produce errors if the results are not real.
+This includes cases such as \axiom{asin(1.2)}, \axiom{log(-3.2)},
+\axiom{sqrt(-1.1)}.
+%
+\xtc{
+The default floating-point type is \axiomType{Float} so to evaluate
+functions using \axiomType{Float} or \axiomType{Complex Float}, just use
+normal decimal notation.
+}{
+\spadpaste{exp(3.1)}
+}
+\xtc{
+}{
+\spadpaste{exp(3.1 + 4.5 * \%i)}
+}
+\xtc{
+To evaluate functions using \axiomType{DoubleFloat}
+or \axiomType{Complex DoubleFloat},
+a declaration or conversion is required.
+}{
+\spadpaste{r: DFLOAT := 3.1; t: DFLOAT := 4.5; exp(r + t*\%i)}
+}
+\xtc{
+}{
+\spadpaste{exp(3.1::DFLOAT + 4.5::DFLOAT * \%i)}
+}
+%
+A number of special functions are provided by the package
+\axiomType{DoubleFloatSpecialFunctions} for the machine-precision
+%-% \HDindex{special functions}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+floating-point types.
+%-% \HDexptypeindex{DoubleFloatSpecialFunctions}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+The special functions provided are listed below, where \axiom{F} stands for
+the types \axiomType{DoubleFloat} and \axiomType{Complex DoubleFloat}.
+The real versions of the functions yield an error if the result is not real.
+%-% \HDindex{function!special}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+
+\noindent
+\axiomFun{Gamma}: \axiom{F -> F}\hfill\newline
+\axiom{Gamma(z)} is the Euler gamma function,
+%-% \HDindex{function!Gamma}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+   \texht{$\Gamma(z)$}{\axiom{Gamma(z)}},
+   defined by
+%-% \HDindex{Euler!gamma function}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+   \texht{\narrowDisplay{\Gamma(z) = \int_{0}^{\infty} t^{z-1} e^{-t} dt.}%
+   }{
+\newline
+\centerline{{\axiom{Gamma(z) = integrate(t**(z-1)*exp(-t), t=0..\%infinity).}}}
+   }
+
+\noindent
+\axiomFun{Beta}: \axiom{F -> F}\hfill\newline
+   \axiom{Beta(u, v)} is the Euler Beta function,
+%-% \HDindex{function!Euler Beta}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+   \texht{$B(u,v)$}{\axiom{B(u,v)}}, defined by
+%-% \HDindex{Euler!Beta function}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+   \texht{\narrowDisplay{B(u,v) = \int_{0}^{1} t^{u-1} (1-t)^{v-1} dt.}%
+   }{
+\newline
+\centerline{{      \axiom{Beta(u,v) = integrate(t**(u-1)*(1-t)**(b-1), t=0..1).}}}
+   }
+   This is related to \texht{$\Gamma(z)$}{\axiom{Gamma(z)}} by
+   \texht{\narrowDisplay{B(u,v) = \frac{\Gamma(u) \Gamma(v)}{\Gamma(u + v)}.}%
+   }{
+\newline
+\centerline{{      \axiom{Beta(u,v) = Gamma(u)*Gamma(v)/Gamma(u + v).}}}
+   }%
+
+\noindent
+\axiomFun{logGamma}: \axiom{F -> F}\hfill\newline
+   \axiom{logGamma(z)} is the natural logarithm of
+\texht{$\Gamma(z)$}{\axiom{Gamma(z)}}.
+   This can often be computed even if \texht{$\Gamma(z)$}{\axiom{Gamma(z)}}
+cannot.
+%
+
+\noindent
+\axiomFun{digamma}: \axiom{F -> F}\hfill\newline
+   \axiom{digamma(z)}, also called \axiom{psi(z)},
+%-% \HDindex{psi @ $\psi$}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+is the function \texht{$\psi(z)$}{\axiom{psi(z)}},
+%-% \HDindex{function!digamma}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+   defined by \texht{\narrowDisplay{\psi(z) = \Gamma'(z)/\Gamma(z).}%
+   }{
+\newline
+\centerline{{      \axiom{psi(z) = Gamma'(z)/Gamma(z).}}}
+  }%
+
+\noindent
+\axiomFun{polygamma}: \axiom{(NonNegativeInteger, F) -> F}\hfill\newline
+   \axiom{polygamma(n, z)} is the \eth{\axiom{n}} derivative of
+%-% \HDindex{function!polygamma}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+   \texht{$\psi(z)$}{\axiom{digamma(z)}}\texht{, written $\psi^{(n)}(z)$}{}.
+
+\noindent
+\axiomFun{besselJ}: \axiom{(F,F) -> F}\hfill\newline
+   \axiom{besselJ(v,z)} is the Bessel function of the first kind,
+%-% \HDindex{function!Bessel}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+   \texht{$J_\nu (z)$}{\axiom{J(v,z)}}.
+   This function satisfies the differential equation
+   \texht{\narrowDisplay{z^2 w''(z) + z w'(z) + (z^2-\nu^2)w(z) = 0.}%
+   }{
+\newline
+\centerline{{      \axiom{z**2*w''(z) + z*w'(z) + (z**2-v**2)*w(z) = 0.}}}
+  }%
+
+\noindent
+\axiomFun{besselY}: \axiom{(F,F) -> F}\hfill\newline
+   \axiom{besselY(v,z)} is the Bessel function of the second kind,
+%-% \HDindex{function!Bessel}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+   \texht{$Y_\nu (z)$}{\axiom{Y(v,z)}}.
+   This function satisfies the same differential equation as
+   \axiomFun{besselJ}.
+   The implementation simply uses the relation
+   \texht{\narrowDisplay{Y_\nu (z) = \frac{J_\nu (z) \cos(\nu \pi) - J_{-\nu} (z)}{\sin(\nu \pi)}.}%
+   }{
+\newline
+\centerline{{      \axiom{Y(v,z) = (J(v,z)*cos(v*\%pi) - J(-v,z))/sin(v*\%pi).}}}
+   }%
+
+\noindent
+\axiomFun{besselI}: \axiom{(F,F) -> F}\hfill\newline
+   \axiom{besselI(v,z)} is the modified Bessel function of the first kind,
+%-% \HDindex{function!Bessel}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+   \texht{$I_\nu (z)$}{\axiom{I(v,z)}}.
+   This function satisfies the differential equation
+   \texht{\narrowDisplay{z^2 w''(z) + z w'(z) - (z^2+\nu^2)w(z) = 0.}%
+   }{
+\newline
+\centerline{{      \axiom{z**2 * w''(z) + z * w'(z) - (z**2+v**2)*w(z) = 0.}}}
+   }%
+
+\noindent
+\axiomFun{besselK}: \axiom{(F,F) -> F}\hfill\newline
+   \axiom{besselK(v,z)} is the modified Bessel function of the second kind,
+%-% \HDindex{function!Bessel}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+   \texht{$K_\nu (z)$}{\axiom{K(v,z)}}.
+   This function satisfies the same differential equation as \axiomFun{besselI}.
+%-% \HDindex{Bessel function}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+   The implementation simply uses the relation
+   \texht{\narrowDisplay{K_\nu (z) = \pi \frac{I_{-\nu} (z) - I_{\nu} (z)}{2 \sin(\nu \pi)}.}%
+   }{
+\newline
+\centerline{{      \axiom{K(v,z) = \%pi*(I(v,z) - I(-v,z))/(2*sin(v*\%pi)).}}}
+   }
+
+\noindent
+\axiomFun{airyAi}: \axiom{F -> F}\hfill\newline
+   \axiom{airyAi(z)} is the Airy function \texht{$Ai(z)$}{\axiom{Ai(z)}}.
+%-% \HDindex{function!Airy Ai}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+   This function satisfies the differential equation
+   \texht{$w''(z) - z w(z) = 0.$}{\axiom{w''(z) - z * w(z) = 0.}}
+   The implementation simply uses the relation
+   \texht{\narrowDisplay{Ai(-z) = \frac{1}{3}\sqrt{z} ( J_{-1/3} (\frac{2}{3}z^{3/2}) + J_{1/3} (\frac{2}{3}z^{3/2}) ).}%
+   }{
+\newline
+\centerline{{      \axiom{Ai(-z) = 1/3 * sqrt(z) * (J(-1/3, 2/3*z**(3/2)) + J(1/3, 2/3*z**(3/2)) ).}}}
+   }%
+
+\noindent
+\axiomFun{airyBi}: \axiom{F -> F}\hfill\newline
+   \axiom{airyBi(z)} is the Airy function \texht{$Bi(z)$}{\axiom{Bi(z)}}.
+%-% \HDindex{function!Airy Bi}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+   This function satisfies the same differential equation as \axiomFun{airyAi}.
+%-% \HDindex{Airy function}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+   The implementation simply uses the relation
+   \texht{\narrowDisplay{Bi(-z) = \frac{1}{3}\sqrt{3 z} ( J_{-1/3} (\frac{2}{3}z^{3/2}) - J_{1/3} (\frac{2}{3}z^{3/2}) ).}%
+   }{
+\newline
+\centerline{{      \axiom{Bi(-z) = 1/3 *sqrt(3*z) * (J(-1/3, 2/3*z**(3/2)) - J(1/3, 2/3*z**(3/2)) ).}}}
+   }
+
+\noindent
+\axiomFun{hypergeometric0F1}: \axiom{(F,F) -> F}\hfill\newline
+   \axiom{hypergeometric0F1(c,z)} is the hypergeometric function
+%-% \HDindex{function!hypergeometric}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+   \texht{${}_0 F_1 ( ; c; z)$}{\axiom{0F1(; c; z)}}.%
+
+\xtc{
+The above special functions are defined only for small floating-point types.
+If you give \axiomType{Float} arguments, they are converted to
+\axiomType{DoubleFloat} by \Language{}.
+}{
+\spadpaste{Gamma(0.5)**2}
+}
+\xtc{
+}{
+\spadpaste{a := 2.1; b := 1.1; besselI(a + \%i*b, b*a + 1)}
+}
+%
+A number of additional operations may be used to compute numerical values.
+These are special polynomial functions that can be evaluated for values in
+any commutative ring \axiom{R}, and in particular for values in any
+floating-point type.
+The following operations are provided by the package
+\axiomType{OrthogonalPolynomialFunctions}:
+%-% \HDexptypeindex{OrthogonalPolynomialFunctions}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+
+\noindent
+\axiomFun{chebyshevT}: \axiom{(NonNegativeInteger, R) -> R}\hbox{}\hfill\newline
+   \axiom{chebyshevT(n,z)} is the \eth{\axiom{n}} Chebyshev polynomial of the first
+   kind, \texht{$T_n (z)$}{\axiom{T[n](z)}}.  These are defined by
+   \texht{\narrowDisplay{\frac{1-t z}{1-2 t z+t^2} = \sum_{n=0}^{\infty} T_n (z) t^n.}%
+   }{
+\newline
+\centerline{{   \axiom{(1-t*z)/(1-2*t*z+t**2) = sum(T[n](z) *t**n, n = 0..).}}}
+   }%
+
+\noindent
+\axiomFun{chebyshevU}: \axiom{(NonNegativeInteger, R) -> R}\hbox{}\hfill\newline
+   \axiom{chebyshevU(n,z)} is the \eth{\axiom{n}} Chebyshev polynomial of the second
+   kind, \texht{$U_n (z)$}{\axiom{U[n](z)}}. These are defined by
+   \texht{\narrowDisplay{\frac{1}{1-2 t z+t^2} = \sum_{n=0}^{\infty} U_n (z) t^n.}%
+   }{
+\newline
+\centerline{{   \axiom{1/(1-2*t*z+t**2) = sum(U[n](z) *t**n, n = 0..).}}}
+   }%
+
+\noindent
+\axiomFun{hermiteH}:   \axiom{(NonNegativeInteger, R) -> R}\hbox{}\hfill\newline
+   \axiom{hermiteH(n,z)} is the \eth{\axiom{n}} Hermite polynomial,
+   \texht{$H_n (z)$}{\axiom{H[n](z)}}.
+   These are defined by
+   \texht{\narrowDisplay{e^{2 t z - t^2} = \sum_{n=0}^{\infty} H_n (z) \frac{t^n}{n!}.}%
+   }{
+\newline
+\centerline{{   \axiom{exp(2*t*z-t**2) = sum(H[n](z)*t**n/n!, n = 0..).}}}
+   }%
+
+\noindent
+\axiomFun{laguerreL}:  \axiom{(NonNegativeInteger, R) -> R}\hbox{}\hfill\newline
+   \axiom{laguerreL(n,z)} is the \eth{\axiom{n}} Laguerre polynomial,
+   \texht{$L_n (z)$}{\axiom{L[n](z)}}.
+   These are defined by
+   \texht{\narrowDisplay{\frac{e^{-\frac{t z}{1-t}}}{1-t} = \sum_{n=0}^{\infty} L_n (z) \frac{t^n}{n!}.}%
+   }{
+\newline
+\centerline{{   \axiom{exp(-t*z/(1-t))/(1-t) = sum(L[n](z)*t**n/n!, n = 0..).}}}
+   }%
+
+\noindent
+\axiomFun{laguerreL}:  \axiom{(NonNegativeInteger, NonNegativeInteger, R) -> R}\hbox{}\hfill\newline
+   \axiom{laguerreL(m,n,z)} is the associated Laguerre polynomial,
+   \texht{$L^m_n (z)$}{\axiom{L<m>[n](z)}}.
+   This is the \eth{\axiom{m}} derivative of \texht{$L_n (z)$}{\axiom{L[n](z)}}.
+
+\noindent
+\axiomFun{legendreP}:  \axiom{(NonNegativeInteger, R) -> R}\hbox{}\hfill\newline
+   \axiom{legendreP(n,z)} is the \eth{\axiom{n}} Legendre polynomial,
+   \texht{$P_n (z)$}{\axiom{P[n](z)}}.  These are defined by
+   \texht{\narrowDisplay{\frac{1}{\sqrt{1-2 t z+t^2}} = \sum_{n=0}^{\infty} P_n (z) t^n.}%
+   }{
+\newline
+\centerline{{   \axiom{1/sqrt(1-2*z*t+t**2) = sum(P[n](z)*t**n, n = 0..).}}}
+   }%
+
+%
+\xtc{
+These operations require non-negative integers for the indices, but otherwise
+the argument can be given as desired.
+}{
+\spadpaste{[chebyshevT(i, z) for i in 0..5]}
+}
+\xtc{
+The expression \axiom{chebyshevT(n,z)} evaluates to the \eth{\axiom{n}} Chebyshev
+%-% \HDindex{polynomial!Chebyshev!of the first kind}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+polynomial of the first kind.
+}{
+\spadpaste{chebyshevT(3, 5.0 + 6.0*\%i)}
+}
+\xtc{
+}{
+\spadpaste{chebyshevT(3, 5.0::DoubleFloat)}
+}
+\xtc{
+The expression \axiom{chebyshevU(n,z)} evaluates to the \eth{\axiom{n}} Chebyshev
+%-% \HDindex{polynomial!Chebyshev!of the second kind}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+polynomial of the second kind.
+}{
+\spadpaste{[chebyshevU(i, z) for i in 0..5]}
+}
+\xtc{
+}{
+\spadpaste{chebyshevU(3, 0.2)}
+}
+\xtc{
+The expression \axiom{hermiteH(n,z)} evaluates to the \eth{\axiom{n}} Hermite
+%-% \HDindex{polynomial!Hermite}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+polynomial.
+}{
+\spadpaste{[hermiteH(i, z) for i in 0..5]}
+}
+\xtc{
+}{
+\spadpaste{hermiteH(100, 1.0)}
+}
+\xtc{
+The expression \axiom{laguerreL(n,z)} evaluates to the \eth{\axiom{n}} Laguerre
+%-% \HDindex{polynomial!Laguerre}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+polynomial.
+}{
+\spadpaste{[laguerreL(i, z) for i in 0..4]}
+}
+\xtc{
+}{
+\spadpaste{laguerreL(4, 1.2)}
+}
+\xtc{
+}{
+\spadpaste{[laguerreL(j, 3, z) for j in 0..4]}
+}
+\xtc{
+}{
+\spadpaste{laguerreL(1, 3, 2.1)}
+}
+\xtc{
+The expression
+%-% \HDindex{polynomial!Legendre}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+\axiom{legendreP(n,z)} evaluates to the \eth{\axiom{n}} Legendre polynomial,
+}{
+\spadpaste{[legendreP(i,z) for i in 0..5]}
+}
+\xtc{
+}{
+\spadpaste{legendreP(3, 3.0*\%i)}
+}
+%
+
+Finally, three number-theoretic polynomial operations may be evaluated.
+%-% \HDindex{number theory}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+The following operations are provided by the package
+\axiomType{NumberTheoreticPolynomialFunctions}.
+%-% \HDexptypeindex{NumberTheoreticPolynomialFunctions}.{ugProblemNumericPage}{8.1.}{Numeric Functions}
+
+\noindent
+\axiomFun{bernoulliB}: \axiom{(NonNegativeInteger, R) -> R} \hbox{}\hfill\newline
+   \axiom{bernoulliB(n,z)} is the \eth{\axiom{n}} Bernoulli polynomial,
+%-% \HDindex{polynomial!Bernoulli}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+   \texht{$B_n (z)$}{\axiom{B[n](z)}}.  These are defined by
+   \texht{\narrowDisplay{\frac{t e^{z t}}{e^t - 1} = \sum_{n=0}^{\infty} B_n (z) \frac{t^n}{n!}.}
+   }{
+\newline
+\centerline{{   \axiom{t*exp(z*t)/(exp t - 1) = sum(B[n](z)*t**n/n! for n - 0..)}}}
+   }%
+
+\noindent
+\axiomFun{eulerE}: \axiom{(NonNegativeInteger, R) -> R} \hbox{}\hfill\newline
+   \axiom{eulerE(n,z)} is the \eth{\axiom{n}} Euler polynomial,
+%-% \HDindex{Euler!polynomial}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+   \texht{$E_n (z)$}{\axiom{E[n](z)}}.  These are defined by
+%-% \HDindex{polynomial!Euler}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+   \texht{\narrowDisplay{\frac{2 e^{z t}}{e^t + 1} = \sum_{n=0}^{\infty} E_n (z) \frac{t^n}{n!}.}%
+   }{
+\newline
+\centerline{{   \axiom{2*exp(z*t)/(exp t + 1) = sum(E[n](z)*t**n/n! for n - 0..)}}}
+   }%
+
+\noindent
+\axiomFun{cyclotomic}: \axiom{(NonNegativeInteger, R) -> R}\hbox{}\hfill\newline
+   \axiom{cyclotomic(n,z)} is the \eth{\axiom{n}} cyclotomic polynomial
+   \texht{$\Phi_n (z)$}{\axiom{phi(n,z)}}.  This is the polynomial whose
+   roots are precisely the primitive \eth{\axiom{n}} roots of unity.
+%-% \HDindex{Euler!totient function}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+   This polynomial has degree given by the Euler totient function
+%-% \HDindex{function!totient}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+   \texht{$\phi(n)$}{\axiom{phi(n)}}.
+
+\xtc{
+The expression \axiom{bernoulliB(n,z)} evaluates to the \eth{\axiom{n}} Bernoulli
+%-% \HDindex{polynomial!Bernouilli}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+polynomial.
+}{
+\spadpaste{bernoulliB(3, z)}
+}
+\xtc{
+}{
+\spadpaste{bernoulliB(3, 0.7 + 0.4 * \%i)}
+}
+\xtc{
+The expression
+%-% \HDindex{polynomial!Euler}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+\axiom{eulerE(n,z)} evaluates to the \eth{\axiom{n}} Euler polynomial.
+}{
+\spadpaste{eulerE(3, z)}
+}
+\xtc{
+}{
+\spadpaste{eulerE(3, 0.7 + 0.4 * \%i)}
+}
+\xtc{
+The expression
+%-% \HDindex{polynomial!cyclotomic}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+\axiom{cyclotomic(n,z)} evaluates to the \eth{\axiom{n}} cyclotomic polynomial.
+%-% \HDindex{cyclotomic polynomial}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+}{
+\spadpaste{cyclotomic(3, z)}
+}
+\xtc{
+}{
+\spadpaste{cyclotomic(3, (-1.0 + 0.0 * \%i)**(2/3))}
+}
+
+Drawing complex functions in \Language{} is presently somewhat
+awkward compared to drawing real functions.
+It is necessary to use the \axiomFun{draw} operations that operate
+on functions rather than expressions.
+
+\psXtc{
+This is the complex exponential function\texht{ (rotated interactively).}{.}
+%-% \HDindex{function!complex exponential}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+When this is displayed in color, the height is the value of the real part of
+the function and the color is the imaginary part.
+Red indicates large negative imaginary values, green indicates imaginary
+values near zero and blue/violet indicates large positive imaginary values.
+}{
+\graphpaste{draw((x,y)+-> real exp complex(x,y), -2..2, -2*\%pi..2*\%pi, colorFunction == (x, y) +->  imag exp complex(x,y), title=="exp(x+\%i*y)", style=="smooth")}
+}{
+\epsffile[0 0 295 295]{../ps/compexp.ps}
+}
+
+\psXtc{
+This is the complex arctangent function.
+%-% \HDindex{function!complex arctangent}{ugProblemNumericPage}{8.1.}{Numeric Functions}
+Again, the height is the real part of the function value but here
+the color indicates the function value's phase.
+The position of the branch cuts are clearly visible and one can
+see that the function is real only for a real argument.
+}{
+\graphpaste{vp := draw((x,y) +-> real atan complex(x,y), -\%pi..\%pi, -\%pi..\%pi, colorFunction==(x,y) +->argument atan complex(x,y), title=="atan(x+\%i*y)", style=="shade"); rotate(vp,-160,-45); vp}
+}{
+\epsffile[0 0 295 295]{../ps/compatan.ps}
+}
+
+\psXtc{
+This is the complex Gamma function.
+}{
+\graphpaste{draw((x,y) +-> max(min(real Gamma complex(x,y),4),-4), -\%pi..\%pi, -\%pi..\%pi, style=="shade", colorFunction == (x,y) +-> argument Gamma complex(x,y), title == "Gamma(x+\%i*y)", var1Steps == 50, var2Steps== 50)}
+}{
+\epsffile[0 0 295 295]{../ps/compgamm.ps}
+}
+
+\psXtc{
+This shows the real Beta function near the origin.
+}{
+\graphpaste{draw(Beta(x,y)/100, x=-1.6..1.7, y = -1.6..1.7, style=="shade", title=="Beta(x,y)", var1Steps==40, var2Steps==40)}
+}{
+\epsffile[0 0 295 295]{../ps/realbeta.ps}
+}
+
+\psXtc{
+This is the Bessel function \texht{$J_\alpha (x)$}{\axiom{J(alpha,x)}}
+for index \texht{$\alpha$}{\axiom{alpha}} in the range \axiom{-6..4} and
+argument \texht{$x$}{\axiom{x}} in the range \axiom{2..14}.
+}{
+\graphpaste{draw((alpha,x) +-> min(max(besselJ(alpha, x+8), -6), 6), -6..4, -6..6, title=="besselJ(alpha,x)", style=="shade", var1Steps==40, var2Steps==40)}
+}{
+\epsffile[0 0 295 295]{../ps/bessel.ps}
+}
+
+\psXtc{
+This is the modified Bessel function
+\texht{$I_\alpha (x)$}{\axiom{I(alpha,x)}}
+evaluated for various real values of the index \texht{$\alpha$}{\axiom{alpha}}
+and fixed argument \texht{$x = 5$}{\axiom{x = 5}}.
+}{
+\graphpaste{draw(besselI(alpha, 5), alpha = -12..12, unit==[5,20])}
+}{
+\epsffile[0 0 295 295]{../ps/modbess.ps}
+}
+
+\psXtc{
+This is similar to the last example
+except the index \texht{$\alpha$}{\axiom{alpha}}
+takes on complex values in a \axiom{6 x 6} rectangle  centered on the origin.
+}{
+\graphpaste{draw((x,y) +-> real besselI(complex(x/20, y/20),5), -60..60, -60..60, colorFunction == (x,y)+-> argument besselI(complex(x/20,y/20),5), title=="besselI(x+i*y,5)", style=="shade")}
+}{
+\epsffile[0 0 295 295]{../ps/modbessc.ps}
+}
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugProblemFactorTitle}{Polynomial Factorization}
+\newcommand{\ugProblemFactorNumber}{8.2.}
+%
+% =====================================================================
+\begin{page}{ugProblemFactorPage}{8.2. Polynomial Factorization}
+% =====================================================================
+\beginscroll
+%
+The \Language{} polynomial factorization
+%-% \HDindex{polynomial!factorization}{ugProblemFactorPage}{8.2.}{Polynomial Factorization}
+facilities are available for all polynomial types and a wide variety of
+coefficient domains.
+%-% \HDindex{factorization}{ugProblemFactorPage}{8.2.}{Polynomial Factorization}
+Here are some examples.
+
+\beginmenu
+    \menudownlink{{8.2.1. Integer and Rational Number Coefficients}}{ugProblemFactorIntRatPage}
+    \menudownlink{{8.2.2. Finite Field Coefficients}}{ugProblemFactorFFPage}
+    \menudownlink{{8.2.3. Simple Algebraic Extension Field Coefficients}}{ugProblemFactorAlgPage}
+    \menudownlink{{8.2.4. Factoring Rational Functions}}{ugProblemFactorRatFunPage}
+\endmenu
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugProblemFactorIntRatTitle}{Integer and Rational Number Coefficients}
+\newcommand{\ugProblemFactorIntRatNumber}{8.2.1.}
+%
+% =====================================================================
+\begin{page}{ugProblemFactorIntRatPage}{8.2.1. Integer and Rational Number Coefficients}
+% =====================================================================
+\beginscroll
+
+\labelSpace{4pc}
+\xtc{
+Polynomials with integer
+%-% \HDindex{polynomial!factorization!integer coefficients}{ugProblemFactorIntRatPage}{8.2.1.}{Integer and Rational Number Coefficients}
+coefficients can be be factored.
+}{
+\spadpaste{v := (4*x**3+2*y**2+1)*(12*x**5-x**3*y+12) \bound{v}}
+}
+\xtc{
+}{
+\spadpaste{factor v \free{v}}
+}
+\xtc{
+Also, \Language{} can factor polynomials with
+%-% \HDindex{polynomial!factorization!rational number coefficients}{ugProblemFactorIntRatPage}{8.2.1.}{Integer and Rational Number Coefficients}
+rational number coefficients.
+}{
+\spadpaste{w := (4*x**3+(2/3)*x**2+1)*(12*x**5-(1/2)*x**3+12) \bound{w}}
+}
+\xtc{
+}{
+\spadpaste{factor w \free{w}}
+}
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugProblemFactorFFTitle}{Finite Field Coefficients}
+\newcommand{\ugProblemFactorFFNumber}{8.2.2.}
+%
+% =====================================================================
+\begin{page}{ugProblemFactorFFPage}{8.2.2. Finite Field Coefficients}
+% =====================================================================
+\beginscroll
+
+Polynomials with coefficients in a finite field
+%-% \HDindex{polynomial!factorization!finite field coefficients}{ugProblemFactorFFPage}{8.2.2.}{Finite Field Coefficients}
+can be also be factored.
+%-% \HDindex{finite field!factoring polynomial with coefficients in}{ugProblemFactorFFPage}{8.2.2.}{Finite Field Coefficients}
+
+\labelSpace{3pc}
+\xtc{
+}{
+\spadpaste{u : POLY(PF(19)) :=3*x**4+2*x**2+15*x+18 \bound{u}}
+}
+\xtc{
+These include the integers mod \axiom{p}, where \axiom{p} is prime, and
+extensions of these fields.
+}{
+\spadpaste{factor u \free{u}}
+}
+\xtc{
+Convert this to have coefficients in the finite
+field with \texht{$19^3$}{\axiom{19**3}} elements.
+See \downlink{``\ugProblemFiniteTitle''}{ugProblemFinitePage} in Section \ugProblemFiniteNumber\ignore{ugProblemFinite} for more information
+about finite fields.
+}{
+\spadpaste{factor(u :: POLY FFX(PF 19,3)) \free{u}}
+}
+%
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugProblemFactorAlgTitle}{Simple Algebraic Extension Field Coefficients}
+\newcommand{\ugProblemFactorAlgNumber}{8.2.3.}
+%
+% =====================================================================
+\begin{page}{ugProblemFactorAlgPage}{8.2.3. Simple Algebraic Extension Field Coefficients}
+% =====================================================================
+\beginscroll
+
+Polynomials with coefficients in simple algebraic extensions
+%-% \HDindex{polynomial!factorization!algebraic extension field coefficients}{ugProblemFactorAlgPage}{8.2.3.}{Simple Algebraic Extension Field Coefficients}
+of the rational numbers can be factored.
+%-% \HDindex{algebraic number}{ugProblemFactorAlgPage}{8.2.3.}{Simple Algebraic Extension Field Coefficients}
+%-% \HDindex{number!algebraic}{ugProblemFactorAlgPage}{8.2.3.}{Simple Algebraic Extension Field Coefficients}
+
+\labelSpace{2pc}
+\xtc{
+Here, \axiom{aa} and \axiom{bb} are symbolic roots of polynomials.
+}{
+\spadpaste{aa := rootOf(aa**2+aa+1) \bound{aa}}
+}
+\xtc{
+}{
+\spadpaste{p:=(x**3+aa**2*x+y)*(aa*x**2+aa*x+aa*y**2)**2 \free{aa}\bound{p}}
+}
+\xtc{
+Note that the second argument to factor can be a list of
+algebraic extensions to factor over.
+}{
+\spadpaste{factor(p,[aa]) \free{p aa}}
+}
+\xtc{
+This factors \axiom{x**2+3} over the integers.
+}{
+\spadpaste{factor(x**2+3)}
+}
+\xtc{
+Factor the same polynomial over the field obtained by adjoining
+\axiom{aa} to the rational numbers.
+}{
+\spadpaste{factor(x**2+3,[aa]) \free{aa}}
+}
+\xtc{
+Factor \axiom{x**6+108} over the same field.
+}{
+\spadpaste{factor(x**6+108,[aa]) \free{aa}}
+}
+\xtc{
+}{
+\spadpaste{bb:=rootOf(bb**3-2) \bound{bb}}
+}
+\xtc{
+}{
+\spadpaste{factor(x**6+108,[bb]) \free{bb}}
+}
+\xtc{
+Factor again over the field obtained by adjoining both \axiom{aa}
+and \axiom{bb} to the rational numbers.
+}{
+\spadpaste{factor(x**6+108,[aa,bb]) \free{aa bb}}
+}
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugProblemFactorRatFunTitle}{Factoring Rational Functions}
+\newcommand{\ugProblemFactorRatFunNumber}{8.2.4.}
+%
+% =====================================================================
+\begin{page}{ugProblemFactorRatFunPage}{8.2.4. Factoring Rational Functions}
+% =====================================================================
+\beginscroll
+
+Since fractions of polynomials form a field, every element (other than zero)
+%-% \HDindex{rational function!factoring}{ugProblemFactorRatFunPage}{8.2.4.}{Factoring Rational Functions}
+divides any other, so there is no useful notion of irreducible factors.
+Thus the \axiomFun{factor} operation is not very useful for fractions
+of polynomials.
+
+\xtc{
+There is, instead, a specific operation \axiomFun{factorFraction}
+that separately factors the numerator and denominator and returns
+a fraction of the factored results.
+}{
+\spadpaste{factorFraction((x**2-4)/(y**2-4))}
+}
+\xtc{
+You can also use \axiomFun{map}. This expression
+applies the \axiomFun{factor} operation
+to the numerator and denominator.
+}{
+\spadpaste{map(factor,(x**2-4)/(y**2-4))}
+}
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugProblemSymRootTitle}{Manipulating Symbolic Roots of a Polynomial}
+\newcommand{\ugProblemSymRootNumber}{8.3.}
+%
+% =====================================================================
+\begin{page}{ugProblemSymRootPage}{8.3. Manipulating Symbolic Roots of a Polynomial}
+% =====================================================================
+\beginscroll
+%
+In this section we show you how to work with one root or all roots
+%-% \HDindex{root!symbolic}{ugProblemSymRootPage}{8.3.}{Manipulating Symbolic Roots of a Polynomial}
+of a polynomial.
+These roots are represented symbolically (as opposed to being
+numeric approximations).
+See \downlink{``\ugxProblemOnePolTitle''}{ugxProblemOnePolPage} in Section \ugxProblemOnePolNumber\ignore{ugxProblemOnePol} and \downlink{``\ugxProblemPolSysTitle''}{ugxProblemPolSysPage} in Section \ugxProblemPolSysNumber\ignore{ugxProblemPolSys} for
+information about solving for the roots of one or more
+polynomials.
+
+\beginmenu
+    \menudownlink{{8.3.1. Using a Single Root of a Polynomial}}{ugxProblemSymRootOnePage}
+    \menudownlink{{8.3.2. Using All Roots of a Polynomial}}{ugxProblemSymRootAllPage}
+\endmenu
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugxProblemSymRootOneTitle}{Using a Single Root of a Polynomial}
+\newcommand{\ugxProblemSymRootOneNumber}{8.3.1.}
+%
+% =====================================================================
+\begin{page}{ugxProblemSymRootOnePage}{8.3.1. Using a Single Root of a Polynomial}
+% =====================================================================
+\beginscroll
+
+Use \axiomFun{rootOf} to get a symbolic root of a polynomial:
+\axiom{rootOf(p, x)} returns a root of \axiom{p(x)}.
+
+\labelSpace{2pc}
+\xtc{
+This creates an algebraic number \axiom{a}.
+%-% \HDindex{algebraic number}{ugxProblemSymRootOnePage}{8.3.1.}{Using a Single Root of a Polynomial}
+%-% \HDindex{number!algebraic}{ugxProblemSymRootOnePage}{8.3.1.}{Using a Single Root of a Polynomial}
+}{
+\spadpaste{a := rootOf(a**4+1,a) \bound{a}}
+}
+\xtc{
+To find the algebraic relation that defines \axiom{a},
+use \axiomFun{definingPolynomial}.
+}{
+\spadpaste{definingPolynomial a \free{a}}
+}
+\xtc{
+You can use \axiom{a} in any further expression,
+including a nested \axiomFun{rootOf}.
+}{
+\spadpaste{b := rootOf(b**2-a-1,b) \free{a}\bound{b}}
+}
+\xtc{
+Higher powers of the roots are automatically reduced during
+calculations.
+}{
+\spadpaste{a + b \free{a b}\bound{c}}
+}
+\xtc{
+}{
+\spadpaste{\% ** 5 \free{c}}
+}
+\xtc{
+The operation \axiomFun{zeroOf} is similar to \axiomFun{rootOf},
+except that it may express the root using radicals in some cases.
+%-% \HDindex{radical}{ugxProblemSymRootOnePage}{8.3.1.}{Using a Single Root of a Polynomial}
+}{
+\spadpaste{rootOf(c**2+c+1,c)}
+}
+\xtc{
+}{
+\spadpaste{zeroOf(d**2+d+1,d)}
+}
+\xtc{
+}{
+\spadpaste{rootOf(e**5-2,e)}
+}
+\xtc{
+}{
+\spadpaste{zeroOf(f**5-2,f)}
+}
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugxProblemSymRootAllTitle}{Using All Roots of a Polynomial}
+\newcommand{\ugxProblemSymRootAllNumber}{8.3.2.}
+%
+% =====================================================================
+\begin{page}{ugxProblemSymRootAllPage}{8.3.2. Using All Roots of a Polynomial}
+% =====================================================================
+\beginscroll
+
+Use \axiomFun{rootsOf} to get all symbolic roots of a polynomial:
+\axiom{rootsOf(p, x)} returns a
+list of all the roots of \axiom{p(x)}.
+If \axiom{p(x)} has a multiple root of order \axiom{n}, then that root
+%-% \HDindex{root!multiple}{ugxProblemSymRootAllPage}{8.3.2.}{Using All Roots of a Polynomial}
+appears \axiom{n} times in the list.
+\texht{\typeout{Make sure these variables are x0 etc}}{}
+
+\xtc{
+Compute all the roots of \axiom{x**4 + 1}.
+}{
+\spadpaste{l := rootsOf(x**4+1,x) \bound{l}}
+}
+\xtc{
+As a side effect, the variables \axiom{\%x0, \%x1} and \axiom{\%x2} are bound
+to the first three roots of \axiom{x**4+1}.
+}{
+\spadpaste{\%x0**5 \free{l}}
+}
+\xtc{
+Although they all satisfy \axiom{x**4 + 1 = 0, \%x0, \%x1,}
+and \axiom{\%x2} are different algebraic numbers.
+To find the algebraic relation that defines each of them,
+use \axiomFun{definingPolynomial}.
+}{
+\spadpaste{definingPolynomial \%x0 \free{l}}
+}
+\xtc{
+}{
+\spadpaste{definingPolynomial \%x1 \free{l}}
+}
+\xtc{
+}{
+\spadpaste{definingPolynomial \%x2 \free{l}}
+}
+\xtc{
+We can check that the sum and product of the roots of \axiom{x**4+1} are
+its trace and norm.
+}{
+\spadpaste{x3 := last l \free{l} \bound{x3}}
+}
+\xtc{
+}{
+\spadpaste{\%x0 + \%x1 + \%x2 + x3 \free{x3}}
+}
+\xtc{
+}{
+\spadpaste{\%x0 * \%x1 * \%x2 * x3 \free{x3}}
+}
+\xtc{
+Corresponding to the pair of operations
+\axiomFun{rootOf}/\axiomFun{zeroOf} in
+\downlink{``\ugxProblemOnePolTitle''}{ugxProblemOnePolPage} in Section \ugxProblemOnePolNumber\ignore{ugxProblemOnePol}, there is
+an operation \axiomFun{zerosOf} that, like \axiomFun{rootsOf},
+computes all the roots
+of a given polynomial, but which expresses some of them in terms of
+radicals.
+}{
+\spadpaste{zerosOf(y**4+1,y) \bound{z}}
+}
+\xtc{
+As you see, only one implicit algebraic number was created
+(\axiom{\%y1}), and its defining equation is this.
+The other three roots are expressed in radicals.
+}{
+\spadpaste{definingPolynomial \%y1 \free{z}}
+}
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugProblemEigenTitle}{Computation of Eigenvalues and Eigenvectors}
+\newcommand{\ugProblemEigenNumber}{8.4.}
+%
+% =====================================================================
+\begin{page}{ugProblemEigenPage}{8.4. Computation of Eigenvalues and Eigenvectors}
+% =====================================================================
+\beginscroll
+%
+In this section we show you
+some of \Language{}'s facilities for computing and
+%-% \HDindex{eigenvalue}{ugProblemEigenPage}{8.4.}{Computation of Eigenvalues and Eigenvectors}
+manipulating eigenvalues and eigenvectors, also called
+%-% \HDindex{eigenvector}{ugProblemEigenPage}{8.4.}{Computation of Eigenvalues and Eigenvectors}
+characteristic values and characteristic vectors,
+%-% \HDindex{characteristic!value}{ugProblemEigenPage}{8.4.}{Computation of Eigenvalues and Eigenvectors}
+respectively.
+%-% \HDindex{characteristic!vector}{ugProblemEigenPage}{8.4.}{Computation of Eigenvalues and Eigenvectors}
+
+\texht{\vskip 4pc}{}
+
+\xtc{
+Let's first create a matrix with integer entries.
+}{
+\spadpaste{m1 := matrix [[1,2,1],[2,1,-2],[1,-2,4]] \bound{m1}}
+}
+\xtc{
+To get a list of the {\it rational} eigenvalues,
+use the operation \axiomFun{eigenvalues}.
+}{
+\spadpaste{leig := eigenvalues(m1) \free{m1} \bound{leig}}
+}
+\xtc{
+Given an explicit eigenvalue, \axiomFun{eigenvector} computes the eigenvectors
+corresponding to it.
+}{
+\spadpaste{eigenvector(first(leig),m1) \free{m1 leig}}
+}
+
+The operation \axiomFun{eigenvectors} returns a list of pairs of values and
+vectors. When an eigenvalue is rational, \Language{} gives you
+the value explicitly; otherwise, its minimal polynomial is given,
+(the polynomial of lowest degree with the eigenvalues as roots),
+together with a parametric representation of the eigenvector using the
+eigenvalue.
+This means that if you ask \Language{} to \axiomFun{solve}
+the minimal polynomial, then you can substitute these roots
+%-% \HDindex{polynomial!minimal}{ugProblemEigenPage}{8.4.}{Computation of Eigenvalues and Eigenvectors}
+into the parametric form of the corresponding eigenvectors.
+%-% \HDindex{minimal polynomial}{ugProblemEigenPage}{8.4.}{Computation of Eigenvalues and Eigenvectors}
+
+\xtc{
+You must be aware that unless an exact eigenvalue has been computed,
+the eigenvector may be badly in error.
+}{
+\spadpaste{eigenvectors(m1) \free{m1}}
+}
+\xtc{
+Another possibility is to use the operation
+\axiomFun{radicalEigenvectors}
+tries to compute explicitly the eigenvectors
+in terms of radicals.
+%-% \HDindex{radical}{ugProblemEigenPage}{8.4.}{Computation of Eigenvalues and Eigenvectors}
+}{
+\spadpaste{radicalEigenvectors(m1) \free{m1}}
+}
+
+Alternatively, \Language{} can compute real or complex approximations to the
+%-% \HDindex{approximation}{ugProblemEigenPage}{8.4.}{Computation of Eigenvalues and Eigenvectors}
+eigenvectors and eigenvalues using the operations \axiomFun{realEigenvectors}
+or \axiomFun{complexEigenvectors}.
+They each take an additional argument \texht{$\epsilon$}{\axiom{epsilon}}
+to specify the ``precision'' required.
+%-% \HDindex{precision}{ugProblemEigenPage}{8.4.}{Computation of Eigenvalues and Eigenvectors}
+In the real case, this means that each approximation will be within
+\texht{$\pm\epsilon$}{plus or minus \axiom{epsilon}} of the actual
+result.
+In the complex case, this means that each approximation will be within
+\texht{$\pm\epsilon$}{plus or minus \axiom{epsilon}} of the actual result
+in each of the real and imaginary parts.
+
+\xtc{
+The precision can be specified as a \axiomType{Float} if the results are
+desired in floating-point notation, or as \axiomType{Fraction Integer} if the
+results are to be expressed using rational (or complex rational) numbers.
+}{
+\spadpaste{realEigenvectors(m1,1/1000) \free{m1}}
+}
+\xtc{
+If an \axiom{n} by \axiom{n} matrix has \axiom{n} distinct eigenvalues (and
+therefore \axiom{n} eigenvectors) the operation \axiomFun{eigenMatrix}
+gives you a matrix of the eigenvectors.
+}{
+\spadpaste{eigenMatrix(m1) \free{m1}}
+}
+\xtc{
+}{
+\spadpaste{m2 := matrix [[-5,-2],[18,7]] \bound{m2}}
+}
+\xtc{
+}{
+\spadpaste{eigenMatrix(m2) \free{m2}}
+}
+%
+%
+\xtc{
+If a symmetric matrix
+%-% \HDindex{matrix!symmetric}{ugProblemEigenPage}{8.4.}{Computation of Eigenvalues and Eigenvectors}
+has a basis of orthonormal eigenvectors, then
+%-% \HDindex{basis!orthonormal}{ugProblemEigenPage}{8.4.}{Computation of Eigenvalues and Eigenvectors}
+\axiomFun{orthonormalBasis} computes a list of these vectors.
+%-% \HDindex{orthonormal basis}{ugProblemEigenPage}{8.4.}{Computation of Eigenvalues and Eigenvectors}
+}{
+\spadpaste{m3 := matrix [[1,2],[2,1]] \bound{m3}}
+}
+\xtc{
+}{
+\spadpaste{orthonormalBasis(m3) \free{m3}}
+}
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugProblemLinPolEqnTitle}{Solution of Linear and Polynomial Equations}
+\newcommand{\ugProblemLinPolEqnNumber}{8.5.}
+%
+% =====================================================================
+\begin{page}{ugProblemLinPolEqnPage}{8.5. Solution of Linear and Polynomial Equations}
+% =====================================================================
+\beginscroll
+%
+In this section we discuss the \Language{} facilities for solving
+systems of linear equations, finding the roots of polynomials and
+%-% \HDindex{linear equation}{ugProblemLinPolEqnPage}{8.5.}{Solution of Linear and Polynomial Equations}
+solving systems of polynomial equations.
+For a discussion of the solution of differential equations, see
+\downlink{``\ugProblemDEQTitle''}{ugProblemDEQPage} in Section \ugProblemDEQNumber\ignore{ugProblemDEQ}.
+
+\beginmenu
+    \menudownlink{{8.5.1. Solution of Systems of Linear Equations}}{ugxProblemLinSysPage}
+    \menudownlink{{8.5.2. Solution of a Single Polynomial Equation}}{ugxProblemOnePolPage}
+    \menudownlink{{8.5.3. Solution of Systems of Polynomial Equations}}{ugxProblemPolSysPage}
+\endmenu
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugxProblemLinSysTitle}{Solution of Systems of Linear Equations}
+\newcommand{\ugxProblemLinSysNumber}{8.5.1.}
+%
+% =====================================================================
+\begin{page}{ugxProblemLinSysPage}{8.5.1. Solution of Systems of Linear Equations}
+% =====================================================================
+\beginscroll
+
+You can use the operation \axiomFun{solve} to solve systems of linear equations.
+%-% \HDindex{equation!linear!solving}{ugxProblemLinSysPage}{8.5.1.}{Solution of Systems of Linear Equations}
+
+The operation \axiomFun{solve} takes two arguments, the list of equations and the
+list of the unknowns to be solved for.
+A system of linear equations need not have a unique solution.
+
+\xtc{
+To solve the linear system:
+\centerline{{\axiom{  x +   y +   z = 8} }}
+\centerline{{\axiom{3*x - 2*y +   z = 0} }}
+\centerline{{\axiom{  x + 2*y + 2*z = 17}}}
+evaluate this expression.
+}{
+\spadpaste{solve([x+y+z=8,3*x-2*y+z=0,x+2*y+2*z=17],[x,y,z])}
+}
+
+Parameters are given as new variables starting with a percent sign and
+\spadSyntax{\%} and
+the variables are expressed in terms of the parameters.
+If the system has no solutions then the empty list is returned.
+
+\xtc{
+When you solve the linear system
+\centerline{{\axiom{  x + 2*y + 3*z = 2} }}
+\centerline{{\axiom{2*x + 3*y + 4*z = 2} }}
+\centerline{{\axiom{3*x + 4*y + 5*z = 2}}}
+with this expression
+you get a solution involving a parameter.
+}{
+\spadpaste{solve([x+2*y+3*z=2,2*x+3*y+4*z=2,3*x+4*y+5*z=2],[x,y,z])}
+}
+
+The system can also be presented as a matrix and a vector.
+The matrix contains the coefficients of the linear equations and
+the vector contains the numbers appearing on the right-hand sides
+of the equations.
+You may input the matrix as a list of rows and the vector as a
+list of its elements.
+
+\xtc{
+To solve the system:
+\centerline{{\axiom{  x +   y +   z = 8} }}
+\centerline{{\axiom{3*x - 2*y +   z = 0} }}
+\centerline{{\axiom{  x + 2*y + 2*z = 17}}}
+in matrix form you would evaluate this expression.
+}{
+\spadpaste{solve([[1,1,1],[3,-2,1],[1,2,2]],[8,0,17])}
+}
+
+The solutions are presented as a \pspadtype{Record} with two
+components: the component \texht{{\it particular}}{{\tt
+particular}} contains a particular solution of the given system or
+the item {\tt "failed"} if there are no solutions, the component
+\texht{{\it basis}}{{\tt basis}} contains a list of vectors that
+are a basis for the space of solutions of the corresponding
+homogeneous system.
+If the system of linear equations does not have a unique solution,
+then the \texht{{\it basis}}{{\tt basis}} component contains
+non-trivial vectors.
+
+\xtc{
+This happens when you solve the linear system
+\centerline{{\axiom{  x + 2*y + 3*z = 2} }}
+\centerline{{\axiom{2*x + 3*y + 4*z = 2} }}
+\centerline{{\axiom{3*x + 4*y + 5*z = 2}}}
+with this command.
+}{
+\spadpaste{solve([[1,2,3],[2,3,4],[3,4,5]],[2,2,2])}
+}
+
+All solutions of this system are obtained by adding the particular
+solution with a linear combination of the \texht{{\it basis}}{{\tt
+basis}} vectors.
+
+When no solution exists then {\tt "failed"} is returned as the
+\texht{{\it particular}}{{\tt particular}} component, as follows:
+
+\xtc{
+}{
+\spadpaste{solve([[1,2,3],[2,3,4],[3,4,5]],[2,3,2])}
+}
+
+When you want to solve a system of homogeneous equations (that is,
+a system where the numbers on the right-hand sides of the
+%-% \HDindex{nullspace}{ugxProblemLinSysPage}{8.5.1.}{Solution of Systems of Linear Equations}
+equations are all zero) in the matrix form you can omit the second
+argument and use the \axiomFun{nullSpace} operation.
+
+\xtc{
+This computes the solutions of the following system of equations:
+\centerline{{\axiom{  x + 2*y + 3*z = 0} }}
+\centerline{{\axiom{2*x + 3*y + 4*z = 0} }}
+\centerline{{\axiom{3*x + 4*y + 5*z = 0}}}
+The result is given as a list of vectors and
+these vectors form a basis for the solution space.
+}{
+\spadpaste{nullSpace([[1,2,3],[2,3,4],[3,4,5]])}
+}
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugxProblemOnePolTitle}{Solution of a Single Polynomial Equation}
+\newcommand{\ugxProblemOnePolNumber}{8.5.2.}
+%
+% =====================================================================
+\begin{page}{ugxProblemOnePolPage}{8.5.2. Solution of a Single Polynomial Equation}
+% =====================================================================
+\beginscroll
+
+\Language{} can solve polynomial equations producing either approximate
+%-% \HDindex{polynomial!root finding}{ugxProblemOnePolPage}{8.5.2.}{Solution of a Single Polynomial Equation}
+or exact solutions.
+%-% \HDindex{equation!polynomial!solving}{ugxProblemOnePolPage}{8.5.2.}{Solution of a Single Polynomial Equation}
+Exact solutions are either members of the ground
+field or can be presented symbolically as roots of irreducible polynomials.
+
+\xtc{
+This returns the one rational root along with an irreducible
+polynomial describing the other solutions.
+}{
+\spadpaste{solve(x**3  = 8,x)}
+}
+\xtc{
+If you want solutions expressed in terms of radicals you would use this
+instead.
+%-% \HDindex{radical}{ugxProblemOnePolPage}{8.5.2.}{Solution of a Single Polynomial Equation}
+}{
+\spadpaste{radicalSolve(x**3  = 8,x)}
+}
+
+The \axiomFun{solve} command always returns a value but
+\axiomFun{radicalSolve} returns only the solutions that it is
+able to express in terms of radicals.
+%-% \HDindex{radical}{ugxProblemOnePolPage}{8.5.2.}{Solution of a Single Polynomial Equation}
+
+If the polynomial equation has rational coefficients
+you can ask for approximations to its real roots by calling
+solve with a second argument that specifies the ``precision''
+%-% \HDindex{precision}{ugxProblemOnePolPage}{8.5.2.}{Solution of a Single Polynomial Equation}
+\texht{$\epsilon$}{\axiom{epsilon}}.
+This means that each approximation will be within
+\texht{$\pm\epsilon$}{plus or minus \axiom{epsilon}} of the actual
+result.
+
+\xtc{
+Notice that the type of second argument controls the type of the result.
+}{
+\spadpaste{solve(x**4 - 10*x**3 + 35*x**2 - 50*x + 25,.0001)}
+}
+\xtc{
+If you give a floating-point precision you get a floating-point result;
+if you give the precision as a rational number you get a rational result.
+}{
+\spadpaste{solve(x**3-2,1/1000)}
+}
+\xtc{
+If you want approximate complex results you should use the
+%-% \HDindex{approximation}{ugxProblemOnePolPage}{8.5.2.}{Solution of a Single Polynomial Equation}
+command \axiomFun{complexSolve} that takes the same precision argument
+\texht{$\epsilon$}{\axiom{epsilon}}.
+}{
+\spadpaste{complexSolve(x**3-2,.0001)}
+}
+\xtc{
+Each approximation will be within
+\texht{$\pm\epsilon$}{plus or minus \axiom{epsilon}} of the actual result
+in each of the real and imaginary parts.
+}{
+\spadpaste{complexSolve(x**2-2*\%i+1,1/100)}
+}
+
+Note that if you omit the \axiomOp{=} from the first argument
+\Language{} generates an equation by equating the first argument to zero.
+Also, when only one variable is present in the equation, you
+do not need to specify the variable to be solved for, that is,
+you can omit the second argument.
+
+\xtc{
+\Language{} can also solve equations involving rational functions.
+Solutions where the denominator vanishes are discarded.
+}{
+\spadpaste{radicalSolve(1/x**3 + 1/x**2 + 1/x = 0,x)}
+}
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugxProblemPolSysTitle}{Solution of Systems of Polynomial Equations}
+\newcommand{\ugxProblemPolSysNumber}{8.5.3.}
+%
+% =====================================================================
+\begin{page}{ugxProblemPolSysPage}{8.5.3. Solution of Systems of Polynomial Equations}
+% =====================================================================
+\beginscroll
+
+Given a system of equations of rational functions with exact coefficients:
+%-% \HDindex{equation!polynomial!solving}{ugxProblemPolSysPage}{8.5.3.}{Solution of Systems of Polynomial Equations}
+\centerline{{\axiom{p1(x1,...,xn)} }}
+\centerline{{\axiom{.}             }}
+\centerline{{\axiom{.}             }}
+\centerline{{\axiom{.}             }}
+\centerline{{\axiom{pm(x1,...,xn)}}}
+\Language{} can find
+numeric or symbolic solutions.
+The system is first split into irreducible components, then for
+each component, a triangular system of equations is found that reduces
+the problem to sequential solution of univariate polynomials resulting
+from substitution of partial solutions from the previous stage.
+\centerline{{\axiom{q1(x1,...,xn)} }}
+\centerline{{\axiom{.} }}
+\centerline{{\axiom{.} }}
+\centerline{{\axiom{.} }}
+\centerline{{\axiom{qm(xn)}}}
+
+Symbolic solutions can be presented using ``implicit'' algebraic numbers
+defined as roots of irreducible polynomials or in terms of radicals.
+\Language{} can also find approximations to the real or complex roots
+of a system of polynomial equations to any user-specified accuracy.
+
+The operation \axiomFun{solve} for systems is used in a way similar
+to \axiomFun{solve} for single equations.
+Instead of a polynomial equation, one has to give a list of
+equations and instead of a single variable to solve for, a list of
+variables.
+For solutions of single equations see \downlink{``\ugxProblemOnePolTitle''}{ugxProblemOnePolPage} in Section \ugxProblemOnePolNumber\ignore{ugxProblemOnePol}.
+
+%
+\xtc{
+Use the operation \axiomFun{solve} if you want implicitly presented
+solutions.
+}{
+\spadpaste{solve([3*x**3 + y + 1,y**2 -4],[x,y])}
+}
+\xtc{
+}{
+\spadpaste{solve([x = y**2-19,y = z**2+x+3,z = 3*x],[x,y,z])}
+}
+\xtc{
+Use \axiomFun{radicalSolve} if you want your solutions expressed
+in terms of radicals.
+}{
+\spadpaste{radicalSolve([3*x**3 + y + 1,y**2 -4],[x,y])}
+}
+
+To get numeric solutions you only need to give the list of
+equations and the precision desired.
+The list of variables would be redundant information since there
+can be no parameters for the numerical solver.
+
+\xtc{
+If the precision is expressed as a floating-point number you get
+results expressed as floats.
+}{
+\spadpaste{solve([x**2*y - 1,x*y**2 - 2],.01)}
+}
+\xtc{
+To get complex numeric solutions, use the operation \axiomFun{complexSolve},
+which takes the same arguments as in the real case.
+}{
+\spadpaste{complexSolve([x**2*y - 1,x*y**2 - 2],1/1000)}
+}
+\xtc{
+It is also possible to solve systems of equations in rational functions
+over the rational numbers.
+Note that \axiom{[x = 0.0, a = 0.0]} is not returned as a solution since
+the denominator vanishes there.
+}{
+\spadpaste{solve([x**2/a = a,a = a*x],.001)}
+}
+\xtc{
+When solving equations with
+denominators, all solutions where the denominator vanishes are
+discarded.
+}{
+\spadpaste{radicalSolve([x**2/a + a + y**3 - 1,a*y + a + 1],[x,y])}
+}
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugProblemLimitsTitle}{Limits}
+\newcommand{\ugProblemLimitsNumber}{8.6.}
+%
+% =====================================================================
+\begin{page}{ugProblemLimitsPage}{8.6. Limits}
+% =====================================================================
+\beginscroll
+%
+To compute a limit, you must specify a functional expression,
+%-% \HDindex{limit}{ugProblemLimitsPage}{8.6.}{Limits}
+a variable, and a limiting value for that variable.
+If you do not specify a direction, \Language{} attempts to
+compute a two-sided limit.
+
+\xtc{
+Issue this to compute the limit
+\texht{$$\lim_{x \rightarrow 1}{{\displaystyle x^2 - 3x +
+2}\over{\displaystyle x^2 - 1}}.$$}{
+of \axiom{(x**2 - 3*x + 2)/(x**2 - 1)} as \axiom{x} approaches \axiom{1}.}
+}{
+\spadpaste{limit((x**2 - 3*x + 2)/(x**2 - 1),x = 1)}
+}
+
+Sometimes the limit when approached from the left is different from the
+limit from the right and, in this case, you may wish to ask for a
+one-sided limit.
+Also,
+if you have a function that is only defined on one side of a particular value,
+%-% \HDindex{limit!one-sided vs. two-sided}{ugProblemLimitsPage}{8.6.}{Limits}
+you can compute a one-sided limit.
+
+\xtc{
+The function \axiom{log(x)} is only defined to the right of zero,
+that is, for \axiom{x > 0}.
+Thus, when computing limits of functions involving \axiom{log(x)},
+you probably want a ``right-hand'' limit.
+}{
+\spadpaste{limit(x * log(x),x = 0,"right")}
+}
+\xtc{
+When you do not specify \axiom{"right"} or \axiom{"left"} as the
+optional fourth argument, \axiomFun{limit} tries to compute a
+two-sided limit.
+Here the limit from the left does not exist, as \Language{}
+indicates when you try to take a two-sided limit.
+}{
+\spadpaste{limit(x * log(x),x = 0)}
+}
+
+A function can be defined on both sides of a particular value, but
+tend to different limits as its variable approaches that value from the
+left and from the right.
+We can construct an example of this as follows:
+Since
+\texht{$\sqrt{y^2}$}{\axiom{sqrt(y**2)}}
+is simply the absolute value of \axiom{y},
+the function
+\texht{$\sqrt{y^2} / y$}{\axiom{sqrt(y**2) / y}}
+is simply the sign (\axiom{+1} or \axiom{-1}) of the nonzero
+real number \axiom{y}.
+Therefore,
+\texht{$\sqrt{y^2} / y = -1$}{\axiom{sqrt(y**2) / y = -1}}
+for \axiom{y < 0} and
+\texht{$\sqrt{y^2} / y = +1$}{\axiom{sqrt(y**2) / y = +1}}
+for \axiom{y > 0}.
+\xtc{
+This is what happens when we take the limit at \axiom{y = 0}.
+The answer returned by \Language{} gives both a
+``left-hand'' and a ``right-hand'' limit.
+}{
+\spadpaste{limit(sqrt(y**2)/y,y = 0)}
+}
+\xtc{
+Here is another example, this time using a more complicated function.
+}{
+\spadpaste{limit(sqrt(1 - cos(t))/t,t = 0)}
+}
+
+You can compute limits at infinity by passing either
+%-% \HDindex{limit!at infinity}{ugProblemLimitsPage}{8.6.}{Limits}
+\texht{$+\infty$ or $-\infty$}{``plus infinity'' or ``minus infinity''}
+as the third argument of \axiomFun{limit}.
+\xtc{
+To do this, use the constants \axiom{\%plusInfinity} and \axiom{\%minusInfinity}.
+}{
+\spadpaste{limit(sqrt(3*x**2 + 1)/(5*x),x = \%plusInfinity)}
+}
+\xtc{
+}{
+\spadpaste{limit(sqrt(3*x**2 + 1)/(5*x),x = \%minusInfinity)}
+}
+
+\xtc{
+You can take limits of functions with parameters.
+%-% \HDindex{limit!of function with parameters}{ugProblemLimitsPage}{8.6.}{Limits}
+As you can see, the limit is expressed in terms of the parameters.
+}{
+\spadpaste{limit(sinh(a*x)/tan(b*x),x = 0)}
+}
+
+When you use \axiomFun{limit}, you are taking the limit of a real
+function of a real variable.
+\xtc{
+When you compute this,
+\Language{} returns \axiom{0} because, as a function of a real variable,
+\axiom{sin(1/z)} is always between \axiom{-1} and \axiom{1}, so \axiom{z * sin(1/z)}
+tends to \axiom{0} as \axiom{z} tends to \axiom{0}.
+}{
+\spadpaste{limit(z * sin(1/z),z = 0)}
+}
+However, as a function of a {\it complex} variable, \axiom{sin(1/z)} is badly
+%-% \HDindex{limit!real vs. complex}{ugProblemLimitsPage}{8.6.}{Limits}
+behaved near \axiom{0} (one says that \axiom{sin(1/z)} has an
+%-% \HDindex{essential singularity}{ugProblemLimitsPage}{8.6.}{Limits}
+{\it essential singularity} at \axiom{z = 0}).
+%-% \HDindex{singularity!essential}{ugProblemLimitsPage}{8.6.}{Limits}
+\xtc{
+When viewed as a function of a complex variable, \axiom{z * sin(1/z)}
+does not approach any limit as \axiom{z} tends to \axiom{0} in the complex plane.
+\Language{} indicates this when we call \axiomFun{complexLimit}.
+}{
+\spadpaste{complexLimit(z * sin(1/z),z = 0)}
+}
+
+%% This is used in chapter 1
+
+You can also take complex limits at infinity, that is, limits of a function of
+\axiom{z} as \axiom{z} approaches infinity on the Riemann sphere.
+Use the symbol \axiom{\%infinity} to denote ``complex infinity.''
+\xtc{
+As above, to compute complex limits rather than real limits, use
+\axiomFun{complexLimit}.
+}{
+\spadpaste{complexLimit((2 + z)/(1 - z),z = \%infinity)}
+}
+\xtc{
+In many cases, a limit of a real function of a real variable
+exists when the corresponding complex limit does not.
+This limit exists.
+}{
+\spadpaste{limit(sin(x)/x,x = \%plusInfinity)}
+}
+\xtc{
+But this limit does not.
+}{
+\spadpaste{complexLimit(sin(x)/x,x = \%infinity)}
+}
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugProblemLaplaceTitle}{Laplace Transforms}
+\newcommand{\ugProblemLaplaceNumber}{8.7.}
+%
+% =====================================================================
+\begin{page}{ugProblemLaplacePage}{8.7. Laplace Transforms}
+% =====================================================================
+\beginscroll
+%
+\Language{} can compute some forward Laplace transforms, mostly
+%-% \HDindex{Laplace transform}{ugProblemLaplacePage}{8.7.}{Laplace Transforms}
+of elementary
+%-% \HDindex{function!elementary}{ugProblemLaplacePage}{8.7.}{Laplace Transforms}
+functions
+%-% \HDindex{transform!Laplace}{ugProblemLaplacePage}{8.7.}{Laplace Transforms}
+not involving logarithms, although some cases of
+special functions are handled.
+\xtc{
+To compute the forward Laplace transform of \axiom{F(t)} with respect to
+\axiom{t} and express the result as \axiom{f(s)}, issue the command
+\axiom{laplace(F(t), t, s)}.
+}{
+\spadpaste{laplace(sin(a*t)*cosh(a*t)-cos(a*t)*sinh(a*t), t, s)}
+}
+\xtc{
+Here are some other non-trivial examples.
+}{
+\spadpaste{laplace((exp(a*t) - exp(b*t))/t, t, s)}
+}
+\xtc{
+}{
+\spadpaste{laplace(2/t * (1 - cos(a*t)), t, s)}
+}
+\xtc{
+}{
+\spadpaste{laplace(exp(-a*t) * sin(b*t) / b**2, t, s)}
+}
+\xtc{
+}{
+\spadpaste{laplace((cos(a*t) - cos(b*t))/t, t, s)}
+}
+\xtc{
+\Language{} also knows about a few special functions.
+}{
+\spadpaste{laplace(exp(a*t+b)*Ei(c*t), t, s)}
+}
+\xtc{
+}{
+\spadpaste{laplace(a*Ci(b*t) + c*Si(d*t), t, s)}
+}
+\xtc{
+When \Language{} does not know about a particular transform,
+it keeps it as a formal transform in the answer.
+}{
+\spadpaste{laplace(sin(a*t) - a*t*cos(a*t) + exp(t**2), t, s)}
+}
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugProblemIntegrationTitle}{Integration}
+\newcommand{\ugProblemIntegrationNumber}{8.8.}
+%
+% =====================================================================
+\begin{page}{ugProblemIntegrationPage}{8.8. Integration}
+% =====================================================================
+\beginscroll
+%
+Integration is the reverse process of differentiation, that is,
+%-% \HDindex{integration}{ugProblemIntegrationPage}{8.8.}{Integration}
+an {\it integral} of a function \axiom{f} with respect to a variable
+\axiom{x} is any function \axiom{g} such that \axiom{D(g,x)}
+is equal to \axiom{f}.
+\xtc{
+The package \axiomType{FunctionSpaceIntegration} provides the top-level
+integration operation, \axiomFunFrom{integrate}{FunctionSpaceIntegration},
+for integrating real-valued elementary functions.
+%-% \HDexptypeindex{FunctionSpaceIntegration}{ugProblemIntegrationPage}{8.8.}{Integration}
+}{
+\spadpaste{integrate(cosh(a*x)*sinh(a*x), x)}
+}
+\xtc{
+Unfortunately, antiderivatives of most functions cannot be expressed in
+terms of elementary functions.
+}{
+\spadpaste{integrate(log(1 + sqrt(a * x + b)) / x, x)}
+}
+Given an elementary function to integrate, \Language{} returns a formal
+integral as above only when it can prove that the integral is not
+elementary and not when it cannot determine the integral.
+In this rare case it prints a message that it cannot
+determine if an elementary integral exists.
+%
+\xtc{
+Similar functions may have antiderivatives
+%-% \HDindex{antiderivative}{ugProblemIntegrationPage}{8.8.}{Integration}
+that look quite different because the form of the antiderivative
+depends on the sign of a constant that appears in the function.
+}{
+\spadpaste{integrate(1/(x**2 - 2),x)}
+}
+\xtc{
+}{
+\spadpaste{integrate(1/(x**2 + 2),x)}
+}
+%
+If the integrand contains parameters, then there may be several possible
+antiderivatives, depending on the signs of expressions of the parameters.
+\xtc{
+In this case \Language{} returns a list of answers that cover all
+the possible cases.
+Here you
+use the answer involving the square root of \axiom{a} when \axiom{a > 0} and
+%-% \HDindex{integration!result as list of real functions}{ugProblemIntegrationPage}{8.8.}{Integration}
+the answer involving the square root of \axiom{-a} when \axiom{a < 0}.
+}{
+\spadpaste{integrate(x**2 / (x**4 - a**2), x)}
+}
+
+If the parameters and the variables of integration can be complex
+numbers rather than real, then the notion of sign is not defined.
+In this case all the possible answers can be expressed as one
+complex function.
+To get that function, rather than a list of real functions, use
+\axiomFunFrom{complexIntegrate}{FunctionSpaceComplexIntegration},
+which is provided by the package
+%-% \HDindex{integration!result as a complex functions}{ugProblemIntegrationPage}{8.8.}{Integration}
+\axiomType{FunctionSpaceComplexIntegration}.
+%-% \HDexptypeindex{FunctionSpaceComplexIntegration}{ugProblemIntegrationPage}{8.8.}{Integration}
+
+\xtc{
+This operation is used for
+integrating complex-valued elementary functions.
+}{
+\spadpaste{complexIntegrate(x**2 / (x**4 - a**2), x)}
+}
+\xtc{
+As with the real case,
+antiderivatives for most complex-valued functions cannot be expressed
+in terms of elementary functions.
+}{
+\spadpaste{complexIntegrate(log(1 + sqrt(a * x + b)) / x, x)}
+}
+
+Sometimes \axiomFun{integrate} can involve symbolic algebraic numbers
+such as those returned by \axiomFunFrom{rootOf}{Expression}.
+To see how to work with these strange generated symbols (such as
+\axiom{\%\%a0}), see \downlink{``\ugxProblemSymRootAllTitle''}{ugxProblemSymRootAllPage} in Section \ugxProblemSymRootAllNumber\ignore{ugxProblemSymRootAll}.
+
+Definite integration is the process of computing the area between
+%-% \HDindex{integration!definite}{ugProblemIntegrationPage}{8.8.}{Integration}
+the \axiom{x}-axis and the curve of a function \axiom{f(x)}.
+The fundamental theorem of calculus states that if \axiom{f} is
+continuous on an interval \axiom{a..b} and if there exists a function \axiom{g}
+that is differentiable on \axiom{a..b} and such that \axiom{D(g, x)}
+is equal to \axiom{f}, then the definite integral of \axiom{f}
+for \axiom{x} in the interval \axiom{a..b} is equal to \axiom{g(b) - g(a)}.
+
+\xtc{
+The package \axiomType{RationalFunctionDefiniteIntegration} provides
+the top-level definite integration operation,
+\axiomFunFrom{integrate}{RationalFunctionDefiniteIntegration},
+for integrating real-valued rational functions.
+}{
+\spadpaste{integrate((x**4 - 3*x**2 + 6)/(x**6-5*x**4+5*x**2+4), x = 1..2)}
+}
+\Language{} checks beforehand that the function you are integrating is
+defined on the interval \axiom{a..b}, and prints an error message if it
+finds that this is not case, as in the following example:
+\begin{verbatim}
+integrate(1/(x**2-2), x = 1..2)
+
+    >> Error detected within library code:
+       Pole in path of integration
+       You are being returned to the top level
+       of the interpreter.
+\end{verbatim}
+When parameters are present in the function, the function may or may not be
+defined on the interval of integration.
+
+\xtc{
+If this is the case, \Language{} issues a warning that a pole might
+lie in the path of integration, and does not compute the integral.
+}{
+\spadpaste{integrate(1/(x**2-a), x = 1..2)}
+}
+
+If you know that you are using values of the parameter for which
+the function has no pole in the interval of integration, use the
+string {\tt "noPole"} as a third argument to
+\axiomFunFrom{integrate}{RationalFunctionDefiniteIntegration}:
+
+%
+\xtc{
+The value here is, of course, incorrect if \axiom{sqrt(a)} is between
+\axiom{1} and \axiom{2.}
+}{
+\spadpaste{integrate(1/(x**2-a), x = 1..2, "noPole")}
+}
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugProblemSeriesTitle}{Working with Power Series}
+\newcommand{\ugProblemSeriesNumber}{8.9.}
+%
+% =====================================================================
+\begin{page}{ugProblemSeriesPage}{8.9. Working with Power Series}
+% =====================================================================
+\beginscroll
+%
+\Language{} has very sophisticated facilities for working with power
+%-% \HDindex{series}{ugProblemSeriesPage}{8.9.}{Working with Power Series}
+series.
+%-% \HDindex{power series}{ugProblemSeriesPage}{8.9.}{Working with Power Series}
+Infinite series are represented by a list of the
+coefficients that have already been determined, together with a
+function for computing the additional coefficients if needed.
+%
+The system command that determines how many terms of a series is displayed
+is \spadcmd{)set streams calculate}.
+For the purposes of this book, we have used this system command to display
+fewer than ten terms.
+%-% \HDsyscmdindex{set streams calculate}{ugProblemSeriesPage}{8.9.}{Working with Power Series}
+Series can be created from expressions, from functions for the
+series coefficients, and from applications of operations on
+existing series.
+The most general function for creating a series is called
+\axiomFun{series}, although you can also use \axiomFun{taylor},
+\axiomFun{laurent} and \axiomFun{puiseux} in situations where you
+know what kind of exponents are involved.
+
+For information about solving differential equations in terms of
+power series, see \downlink{``\ugxProblemDEQSeriesTitle''}{ugxProblemDEQSeriesPage} in Section \ugxProblemDEQSeriesNumber\ignore{ugxProblemDEQSeries}.
+
+\beginmenu
+    \menudownlink{{8.9.1. Creation of Power Series}}{ugxProblemSeriesCreatePage}
+    \menudownlink{{8.9.2. Coefficients of Power Series}}{ugxProblemSeriesCoefficientsPage}
+    \menudownlink{{8.9.3. Power Series Arithmetic}}{ugxProblemSeriesArithmeticPage}
+    \menudownlink{{8.9.4. Functions on Power Series}}{ugxProblemSeriesFunctionsPage}
+    \menudownlink{{8.9.5. Converting to Power Series}}{ugxProblemSeriesConversionsPage}
+    \menudownlink{{8.9.6. Power Series from Formulas}}{ugxProblemSeriesFormulaPage}
+    \menudownlink{{8.9.7. Substituting Numerical Values in Power Series}}{ugxProblemSeriesSubstitutePage}
+    \menudownlink{{8.9.8. Example: Bernoulli Polynomials and Sums of Powers}}{ugxProblemSeriesBernoulliPage}
+\endmenu
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugxProblemSeriesCreateTitle}{Creation of Power Series}
+\newcommand{\ugxProblemSeriesCreateNumber}{8.9.1.}
+%
+% =====================================================================
+\begin{page}{ugxProblemSeriesCreatePage}{8.9.1. Creation of Power Series}
+% =====================================================================
+\beginscroll
+
+\labelSpace{4pc}
+\xtc{
+This is the easiest way to create a power series.
+This tells \Language{} that \axiom{x} is to be treated as a power series,
+%-% \HDindex{series!creating}{ugxProblemSeriesCreatePage}{8.9.1.}{Creation of Power Series}
+so functions of \axiom{x} are again power series.
+}{
+\spadpaste{x := series 'x \bound{x}}
+}
+%
+We didn't say anything about the coefficients of the power
+series, so the coefficients are general expressions over the integers.
+This allows us to introduce denominators, symbolic constants, and other
+variables as needed.
+\xtc{
+Here the coefficients are integers (note that the coefficients
+are the Fibonacci
+%-% \HDindex{Fibonacci numbers}{ugxProblemSeriesCreatePage}{8.9.1.}{Creation of Power Series}
+numbers).
+}{
+\spadpaste{1/(1 - x - x**2) \free{x}}
+}
+\xtc{
+This series has coefficients that are rational numbers.
+}{
+\spadpaste{sin(x) \free{x}}
+}
+\xtc{
+When you enter this expression
+you introduce the symbolic constants \axiom{sin(1)} and \axiom{cos(1).}
+}{
+\spadpaste{sin(1 + x) \free{x}}
+}
+\xtc{
+When you enter the expression
+the variable \axiom{a} appears in the resulting series expansion.
+}{
+\spadpaste{sin(a * x) \free{x}}
+}
+
+\xtc{
+You can also convert an expression into a series expansion.
+This expression creates the series expansion of \axiom{1/log(y)}
+about \axiom{y = 1}.
+For details and more examples, see
+\downlink{``\ugxProblemSeriesConversionsTitle''}{ugxProblemSeriesConversionsPage} in Section \ugxProblemSeriesConversionsNumber\ignore{ugxProblemSeriesConversions}.
+}{
+\spadpaste{series(1/log(y),y = 1)}
+}
+
+You can create power series with more general coefficients.
+You normally accomplish this via a type declaration (see
+\downlink{``\ugTypesDeclareTitle''}{ugTypesDeclarePage} in Section \ugTypesDeclareNumber\ignore{ugTypesDeclare}).
+See \downlink{``\ugxProblemSeriesFunctionsTitle''}{ugxProblemSeriesFunctionsPage} in Section \ugxProblemSeriesFunctionsNumber\ignore{ugxProblemSeriesFunctions} for some warnings about
+working with declared series.
+
+\xtc{
+We declare that \axiom{y} is a one-variable Taylor series
+%-% \HDindex{series!Taylor}{ugxProblemSeriesCreatePage}{8.9.1.}{Creation of Power Series}
+(\axiomType{UTS} is the abbreviation for \axiomType{UnivariateTaylorSeries})
+in the variable \axiom{z} with \axiomType{FLOAT}
+(that is, floating-point) coefficients, centered about \axiom{0.}
+Then, by assignment, we obtain the Taylor expansion of
+\axiom{exp(z)} with floating-point coefficients.
+%-% \HDexptypeindex{UnivariateTaylorSeries}{ugxProblemSeriesCreatePage}{8.9.1.}{Creation of Power Series}
+}{
+\spadpaste{y : UTS(FLOAT,'z,0) := exp(z) \bound{y}}
+}
+
+You can also create a power series by giving an explicit formula
+for its \eth{\axiom{n}} coefficient.
+For details and more examples, see
+\downlink{``\ugxProblemSeriesFormulaTitle''}{ugxProblemSeriesFormulaPage} in Section \ugxProblemSeriesFormulaNumber\ignore{ugxProblemSeriesFormula}.
+
+\xtc{
+To create a series about
+\axiom{w = 0} whose \eth{\axiom{n}} Taylor coefficient
+is \axiom{1/n!}, you can evaluate this expression.
+This is the Taylor expansion of \axiom{exp(w)} at \axiom{w = 0}.
+}{
+\spadpaste{series(1/factorial(n),n,w = 0)}
+}
+%
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugxProblemSeriesCoefficientsTitle}{Coefficients of Power Series}
+\newcommand{\ugxProblemSeriesCoefficientsNumber}{8.9.2.}
+%
+% =====================================================================
+\begin{page}{ugxProblemSeriesCoefficientsPage}{8.9.2. Coefficients of Power Series}
+% =====================================================================
+\beginscroll
+
+You can extract any coefficient from a power series---even one
+that hasn't been computed yet.
+This is possible because in \Language{}, infinite series are
+represented by a list of the coefficients that have already been
+determined, together with a function for computing the additional
+coefficients.
+(This is known as {\it lazy evaluation}.) When you ask for a
+%-% \HDindex{series!lazy evaluation}{ugxProblemSeriesCoefficientsPage}{8.9.2.}{Coefficients of Power Series}
+coefficient that hasn't yet been computed, \Language{} computes
+%-% \HDindex{lazy evaluation}{ugxProblemSeriesCoefficientsPage}{8.9.2.}{Coefficients of Power Series}
+whatever additional coefficients it needs and then stores them in
+the representation of the power series.
+
+\xtc{
+Here's an example of how to extract the coefficients of a power series.
+%-% \HDindex{series!extracting coefficients}{ugxProblemSeriesCoefficientsPage}{8.9.2.}{Coefficients of Power Series}
+}{
+\spadpaste{x := series(x) \bound{x}}
+}
+\xtc{
+}{
+\spadpaste{y := exp(x) * sin(x) \free{x} \bound{y}}
+}
+\xtc{
+This coefficient is readily available.
+}{
+\spadpaste{coefficient(y,6) \free{y}}
+}
+\xtc{
+But let's get the fifteenth coefficient of \axiom{y}.
+}{
+\spadpaste{coefficient(y,15) \free{y} \bound{y15}}
+}
+\xtc{
+If you look at \axiom{y}
+then you see that the coefficients up to order \axiom{15}
+have all been computed.
+}{
+\spadpaste{y \free{y15}}
+}
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugxProblemSeriesArithmeticTitle}{Power Series Arithmetic}
+\newcommand{\ugxProblemSeriesArithmeticNumber}{8.9.3.}
+%
+% =====================================================================
+\begin{page}{ugxProblemSeriesArithmeticPage}{8.9.3. Power Series Arithmetic}
+% =====================================================================
+\beginscroll
+
+You can manipulate power series using the usual arithmetic operations
+%-% \HDindex{series!arithmetic}{ugxProblemSeriesArithmeticPage}{8.9.3.}{Power Series Arithmetic}
+\axiomOpFrom{+}{UnivariatePuiseuxSeries},
+\axiomOpFrom{-}{UnivariatePuiseuxSeries},
+\axiomOpFrom{*}{UnivariatePuiseuxSeries}, and
+\axiomOpFrom{/}{UnivariatePuiseuxSeries}.
+%
+
+\labelSpace{1pc}
+\xtc{
+The results of these operations are also power series.
+}{
+\spadpaste{x := series x \bound{x}}
+}
+\xtc{
+}{
+\spadpaste{(3 + x) / (1 + 7*x)}
+}
+\xtc{
+You can also compute \axiom{f(x) ** g(x)}, where \axiom{f(x)} and \axiom{g(x)}
+are two power series.
+}{
+\spadpaste{base := 1 / (1 - x) \free{x} \bound{base}}
+}
+%
+\xtc{
+}{
+\spadpaste{expon := x * base \free{x base} \bound{expon}}
+}
+%
+\xtc{
+}{
+\spadpaste{base ** expon \free{base expon}}
+}
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugxProblemSeriesFunctionsTitle}{Functions on Power Series}
+\newcommand{\ugxProblemSeriesFunctionsNumber}{8.9.4.}
+%
+% =====================================================================
+\begin{page}{ugxProblemSeriesFunctionsPage}{8.9.4. Functions on Power Series}
+% =====================================================================
+\beginscroll
+
+Once you have created a power series, you can apply transcendental
+functions
+(for example, \axiomFun{exp}, \axiomFun{log}, \axiomFun{sin}, \axiomFun{tan},
+\axiomFun{cosh}, etc.) to it.
+
+\labelSpace{1pc}
+%
+\xtc{
+To demonstrate this, we first create the power series
+expansion of the rational function
+\texht{
+${\displaystyle x^2} \over {\displaystyle 1 - 6x + x^2}$
+}{
+\axiom{x**2/(1 - 6*x + x**2)}
+}
+about \axiom{x = 0}.
+}{
+\spadpaste{x := series 'x \bound{x}}
+}
+%
+\xtc{
+}{
+\spadpaste{rat := x**2 / (1 - 6*x + x**2) \free{x} \bound{rat}}
+}
+%
+%
+\xtc{
+If you want to compute the series expansion of
+\texht{
+$\sin\left({\displaystyle x^2} \over {\displaystyle 1 - 6x + x^2}\right)$
+}{
+\axiom{sin(x**2/(1 - 6*x + x**2))}
+}
+you simply compute the sine of \axiom{rat}.
+}{
+\spadpaste{sin(rat) \free{rat}}
+}
+%
+
+\beginImportant
+\noindent {\bf Warning:}
+the type of the coefficients of a power series may
+affect the kind of computations that you can do with that series.
+This can only happen when you have made a declaration to
+specify a series domain with a certain type of coefficient.
+\endImportant
+
+\xtc{
+If you evaluate
+then you have declared that \axiom{y} is a one variable Taylor series
+%-% \HDindex{series!Taylor}{ugxProblemSeriesFunctionsPage}{8.9.4.}{Functions on Power Series}
+(\axiomType{UTS} is the abbreviation for \axiomType{UnivariateTaylorSeries})
+in the variable \axiom{y} with \axiomType{FRAC INT}
+(that is, fractions of integer) coefficients, centered about \axiom{0}.
+}{
+\spadpaste{y : UTS(FRAC INT,y,0) := y \bound{y}}
+}
+%
+\xtc{
+You can now compute certain power series in \axiom{y},
+{\it provided} that these series have rational coefficients.
+}{
+\spadpaste{exp(y) \free{y}}
+}
+%
+\xtc{
+You can get examples of such series
+by applying transcendental functions to
+series in \axiom{y} that have no constant terms.
+}{
+\spadpaste{tan(y**2) \free{y}}
+}
+%
+\xtc{
+}{
+\spadpaste{cos(y + y**5) \free{y}}
+}
+%
+%
+\xtc{
+Similarly, you can compute the logarithm of a power series with rational
+coefficients if the constant coefficient is \axiom{1.}
+}{
+\spadpaste{log(1 + sin(y)) \free{y}}
+}
+%
+If you wanted to apply, say, the operation \axiomFun{exp} to a power
+series with a nonzero constant coefficient \texht{$a_0$}{\axiom{a0}},
+then the constant coefficient of the result would be
+\texht{$e^{a_0}$}{\axiom{exp(a0)}}, which is {\it not} a rational number.
+Therefore, evaluating \axiom{exp(2 + tan(y))} would generate an error
+message.
+
+If you want to compute the Taylor expansion of \axiom{exp(2 + tan(y))}, you must
+ensure that the coefficient domain has an operation \axiomFun{exp} defined
+for it.
+An example of such a domain is \axiomType{Expression Integer}, the type
+of formal functional expressions over the integers.
+%
+\xtc{
+When working with coefficients of this type,
+}{
+\spadpaste{z : UTS(EXPR INT,z,0) := z \bound{z}}
+}
+\xtc{
+this presents no problems.
+}{
+\spadpaste{exp(2 + tan(z)) \free{z}}
+}
+%
+Another way to create Taylor series whose coefficients are expressions over
+the integers is to use \axiomFun{taylor} which works similarly to
+%-% \HDindex{series!Taylor}{ugxProblemSeriesFunctionsPage}{8.9.4.}{Functions on Power Series}
+\axiomFun{series}.
+%
+\xtc{
+This is equivalent to the previous computation, except that now we
+are using the variable \axiom{w} instead of \axiom{z}.
+}{
+\spadpaste{w := taylor 'w \bound{w}}
+}
+\xtc{
+}{
+\spadpaste{exp(2 + tan(w)) \free{w}}
+}
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugxProblemSeriesConversionsTitle}{Converting to Power Series}
+\newcommand{\ugxProblemSeriesConversionsNumber}{8.9.5.}
+%
+% =====================================================================
+\begin{page}{ugxProblemSeriesConversionsPage}{8.9.5. Converting to Power Series}
+% =====================================================================
+\beginscroll
+
+The \axiomType{ExpressionToUnivariatePowerSeries} package provides
+operations for computing series expansions of functions.
+%-% \HDexptypeindex{ExpressionToUnivariatePowerSeries}{ugxProblemSeriesConversionsPage}{8.9.5.}{Converting to Power Series}
+
+\labelSpace{1pc}
+\xtc{
+Evaluate this
+to compute the Taylor expansion of \axiom{sin x} about
+%-% \HDindex{series!Taylor}{ugxProblemSeriesConversionsPage}{8.9.5.}{Converting to Power Series}
+\axiom{x = 0}.
+The first argument, \axiom{sin(x)}, specifies the function whose series
+expansion is to be computed and the second argument, \axiom{x = 0},
+specifies that the series is to be expanded in power of \axiom{(x - 0)},
+that is, in power of \axiom{x}.
+}{
+\spadpaste{taylor(sin(x),x = 0)}
+}
+\xtc{
+Here is the Taylor expansion of \axiom{sin x} about
+\texht{$x = \frac{\pi}{6}$}{\axiom{x = \%pi/6}}:
+}{
+\spadpaste{taylor(sin(x),x = \%pi/6)}
+}
+
+The function to be expanded into a series may have variables other than
+%-% \HDindex{series!multiple variables}{ugxProblemSeriesConversionsPage}{8.9.5.}{Converting to Power Series}
+the series variable.
+%
+\xtc{
+For example, we may expand \axiom{tan(x*y)} as a Taylor series in
+\axiom{x}
+}{
+\spadpaste{taylor(tan(x*y),x = 0)}
+}
+%
+\xtc{
+or as a Taylor series in \axiom{y}.
+}{
+\spadpaste{taylor(tan(x*y),y = 0)}
+}
+\xtc{
+A more interesting function is
+\texht{${\displaystyle t e^{x t}} \over {\displaystyle e^t - 1}$}{(t *
+\%e**(x*t))/(\%e**t - 1)}.
+When we expand this function as a Taylor series in \axiom{t}
+the \eth{\axiom{n}} order coefficient is the \eth{\axiom{n}} Bernoulli
+%-% \HDindex{Bernoulli!polynomial}{ugxProblemSeriesConversionsPage}{8.9.5.}{Converting to Power Series}
+polynomial
+%-% \HDindex{polynomial!Bernoulli}{ugxProblemSeriesConversionsPage}{8.9.5.}{Converting to Power Series}
+divided by \axiom{n!}.
+}{
+\spadpaste{bern := taylor(t*exp(x*t)/(exp(t) - 1),t = 0) \bound{bern}}
+}
+\xtc{
+Therefore, this and the next expression
+produce the same result.
+}{
+\spadpaste{factorial(6) * coefficient(bern,6) \free{bern}}
+}
+\xtc{
+}{
+\spadpaste{bernoulliB(6,x)}
+}
+
+Technically, a series with terms of negative degree is not considered to
+be a Taylor series, but, rather, a
+%-% \HDindex{series!Laurent}{ugxProblemSeriesConversionsPage}{8.9.5.}{Converting to Power Series}
+{\it Laurent series}.
+%-% \HDindex{Laurent series}{ugxProblemSeriesConversionsPage}{8.9.5.}{Converting to Power Series}
+If you try to compute a Taylor series expansion of
+\texht{$\frac{x}{\log x}$}{x/log(x)}
+at \axiom{x = 1} via \axiom{taylor(x/log(x),x = 1)}
+you get an error message.
+The reason is that the function has a {\it pole} at \axiom{x =
+1}, meaning that
+its series expansion about this point has terms of negative degree.
+A series with finitely many terms of negative degree is called a Laurent
+series.
+\xtc{
+You get the desired series expansion by issuing this.
+}{
+\spadpaste{laurent(x/log(x),x = 1)}
+}
+
+Similarly, a series with terms of fractional degree is neither a Taylor
+series nor a Laurent series.
+Such a series is called a
+%-% \HDindex{series!Puiseux}{ugxProblemSeriesConversionsPage}{8.9.5.}{Converting to Power Series}
+{\it Puiseux series}.
+%-% \HDindex{Puiseux series}{ugxProblemSeriesConversionsPage}{8.9.5.}{Converting to Power Series}
+The expression \axiom{laurent(sqrt(sec(x)),x = 3 * \%pi/2)}
+results in an error message
+because the series expansion about this point has terms of fractional degree.
+\xtc{
+However, this command produces what you want.
+}{
+\spadpaste{puiseux(sqrt(sec(x)),x = 3 * \%pi/2)}
+}
+
+Finally, consider the case of functions that do not have Puiseux
+expansions about certain points.
+An example of this is \texht{$x^x$}{\axiom{x^x}} about \axiom{x = 0}.
+\axiom{puiseux(x**x,x=0)}
+produces an error message because of the
+type of singularity of the function at \axiom{x = 0}.
+\xtc{
+The general function \axiomFun{series} can be used in this case.
+Notice that the series returned is not, strictly speaking, a power series
+because of the \axiom{log(x)} in the expansion.
+}{
+\spadpaste{series(x**x,x=0)}
+}
+
+\beginImportant
+The operation \axiomFun{series} returns the most general type of
+infinite series.
+The user who is not interested in distinguishing
+between various types of infinite series may wish to use this operation
+exclusively.
+\endImportant
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugxProblemSeriesFormulaTitle}{Power Series from Formulas}
+\newcommand{\ugxProblemSeriesFormulaNumber}{8.9.6.}
+%
+% =====================================================================
+\begin{page}{ugxProblemSeriesFormulaPage}{8.9.6. Power Series from Formulas}
+% =====================================================================
+\beginscroll
+
+The \axiomType{GenerateUnivariatePowerSeries} package enables you to
+%-% \HDindex{series!giving formula for coefficients}{ugxProblemSeriesFormulaPage}{8.9.6.}{Power Series from Formulas}
+create power series from explicit formulas for their
+\eth{\axiom{n}} coefficients.
+In what follows, we construct series expansions for certain
+transcendental functions by giving formulas for their
+coefficients.
+You can also compute such series expansions directly simply by
+specifying the function and the point about which the series is to
+be expanded.
+%-% \HDexptypeindex{GenerateUnivariatePowerSeries}{ugxProblemSeriesFormulaPage}{8.9.6.}{Power Series from Formulas}
+See \downlink{``\ugxProblemSeriesConversionsTitle''}{ugxProblemSeriesConversionsPage} in Section \ugxProblemSeriesConversionsNumber\ignore{ugxProblemSeriesConversions} for more information.
+
+Consider the Taylor expansion of \texht{$e^x$}{\axiom{\%e**x}}
+%-% \HDindex{series!Taylor}{ugxProblemSeriesFormulaPage}{8.9.6.}{Power Series from Formulas}
+about \axiom{x = 0}:
+\texht{\narrowDisplay{%
+\begin{array}{ccl}
+e^x &=& \displaystyle 1 + x + \frac{x^2}{2} + \frac{x^3}{6} + \cdots \\ \\
+    &=& \displaystyle\sum_{n=0}^\infty \frac{x^n}{n!}
+\end{array}}%
+}{
+\begin{verbatim}
+%e**x = 1 + x + x**2/2 + x**3/6 + ...
+      = sum from n=0 to n=%infty of x**n / n!
+\end{verbatim}
+}
+The \eth{\axiom{n}} Taylor coefficient is \axiom{1/n!}.
+%
+\xtc{
+This is how you create this series in \Language{}.
+}{
+\spadpaste{series(n +-> 1/factorial(n),x = 0)}
+}
+
+The first argument specifies a formula for the \eth{\axiom{n}}
+coefficient by giving a function that maps \axiom{n} to
+\axiom{1/n!}.
+The second argument specifies that the series is to be expanded in
+powers of \axiom{(x - 0)}, that is, in powers of \axiom{x}.
+Since we did not specify an initial degree, the first term in the
+series was the term of degree 0 (the constant term).
+Note that the formula was given as an anonymous function.
+These are discussed in \downlink{``\ugUserAnonTitle''}{ugUserAnonPage} in Section \ugUserAnonNumber\ignore{ugUserAnon}.
+
+Consider the Taylor expansion of \axiom{log x} about \axiom{x = 1}:
+\texht{\narrowDisplay{%
+\begin{array}{ccl}
+\log(x) &=& \displaystyle (x - 1) - \frac{(x - 1)^2}{2} + \frac{(x - 1)^3}{3} - \cdots \\ \\
+        &=& \displaystyle\sum_{n = 1}^\infty (-1)^{n-1} \frac{(x - 1)^n}{n}
+\end{array}}%
+}{
+\begin{verbatim}
+log x = (x - 1) - \(x - 1)^2/2 + (x - 1)^3/3 - ...
+      = sum from n=1 to n=%infty of (-1)**(n-1) (x - 1)**n / n
+\end{verbatim}
+}
+If you were to evaluate the expression
+\axiom{series(n +-> (-1)**(n-1) / n, x = 1)}
+you would get an error message because \Language{} would try to
+calculate a term of degree \axiom{0} and therefore divide by \axiom{0.}
+
+\xtc{
+Instead, evaluate this.
+The third argument, \axiom{1..}, indicates that only terms of degree
+\axiom{n = 1, ...} are to be computed.
+}{
+\spadpaste{series(n +-> (-1)**(n-1)/n,x = 1,1..)}
+}
+%
+
+Next consider the Taylor expansion of an odd function, say,
+\axiom{sin(x)}:
+\begin{verbatim}
+sin x = x - x**3/3! + x**5/5! - ...
+\end{verbatim}
+Here every other coefficient is zero and we would like to give an
+explicit formula only for the odd Taylor coefficients.
+%
+\xtc{
+This is one way to do it.
+The third argument, \axiom{1..}, specifies that the first term to be computed
+is the term of degree 1.
+The fourth argument, \axiom{2}, specifies that we
+increment by \axiom{2} to find the degrees of subsequent terms, that is, the next
+term
+is of degree \axiom{1 + 2}, the next of degree \axiom{1 + 2 + 2}, etc.
+}{
+\spadpaste{series(n +-> (-1)**((n-1)/2)/factorial(n),x = 0,1..,2)}
+}
+%
+
+\xtc{
+The initial degree and the increment do not have to be integers.
+For example, this expression produces a series expansion of
+\texht{$\sin(x^{\frac{1}{3}})$}{\axiom{sin(x**(1/3))}}.
+}{
+\spadpaste{series(n +-> (-1)**((3*n-1)/2)/factorial(3*n),x = 0,1/3..,2/3)}
+}
+\xtc{
+While the increment must be positive, the initial degree may be negative.
+This yields the Laurent expansion of \axiom{csc(x)} at
+\axiom{x = 0}.
+%
+}{
+\spadpaste{cscx := series(n +-> (-1)**((n-1)/2) * 2 * (2**n-1) * bernoulli(numer(n+1)) / factorial(n+1), x=0, -1..,2) \bound{cscx}}
+}
+\xtc{
+Of course, the reciprocal of this power series is the Taylor expansion
+of \axiom{sin(x)}.
+}{
+\spadpaste{1/cscx \free{cscx}}
+}
+%
+\xtc{
+As a final example,
+here is the Taylor expansion of \axiom{asin(x)} about \axiom{x = 0}.
+}{
+\spadpaste{asinx := series(n +-> binomial(n-1,(n-1)/2)/(n*2**(n-1)),x=0,1..,2) \bound{asinx}}
+}
+\xtc{
+When we compute the \axiom{sin} of this series, we get \axiom{x}
+(in the sense that all higher terms computed so far are zero).
+}{
+\spadpaste{sin(asinx) \free{asinx}}
+}
+
+As we discussed in \downlink{``\ugxProblemSeriesConversionsTitle''}{ugxProblemSeriesConversionsPage} in Section \ugxProblemSeriesConversionsNumber\ignore{ugxProblemSeriesConversions},
+you can also use the operations \axiomFun{taylor}, \axiomFun{laurent} and
+\axiomFun{puiseux} instead of \axiomFun{series} if you know ahead of time
+what kind of exponents a series has.
+You can't go wrong using \axiomFun{series}, though.
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugxProblemSeriesSubstituteTitle}{Substituting Numerical Values in Power Series}
+\newcommand{\ugxProblemSeriesSubstituteNumber}{8.9.7.}
+%
+% =====================================================================
+\begin{page}{ugxProblemSeriesSubstitutePage}{8.9.7. Substituting Numerical Values in Power Series}
+% =====================================================================
+\beginscroll
+
+Use \axiomFunFrom{eval}{UnivariatePowerSeriesCategory}
+%-% \HDindex{approximation}{ugxProblemSeriesSubstitutePage}{8.9.7.}{Substituting Numerical Values in Power Series}
+to substitute a numerical value for a variable in
+%-% \HDindex{series!numerical approximation}{ugxProblemSeriesSubstitutePage}{8.9.7.}{Substituting Numerical Values in Power Series}
+a power series.
+For example, here's a way to obtain numerical approximations of
+\axiom{\%e} from the Taylor series
+expansion of \axiom{exp(x)}.
+
+\labelSpace{1pc}
+\xtc{
+First you create the desired Taylor expansion.
+}{
+\spadpaste{f := taylor(exp(x)) \bound{f}}
+}
+\xtc{
+Then you evaluate the series at the value \axiom{1.0}.
+The result is a sequence of the partial sums.
+}{
+\spadpaste{eval(f,1.0)}
+}
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugxProblemSeriesBernoulliTitle}{Example: Bernoulli Polynomials and Sums of Powers}
+\newcommand{\ugxProblemSeriesBernoulliNumber}{8.9.8.}
+%
+% =====================================================================
+\begin{page}{ugxProblemSeriesBernoulliPage}{8.9.8. Example: Bernoulli Polynomials and Sums of Powers}
+% =====================================================================
+\beginscroll
+
+\Language{} provides operations for computing definite and
+%-% \HDindex{summation!definite}{ugxProblemSeriesBernoulliPage}{8.9.8.}{Example: Bernoulli Polynomials and Sums of Powers}
+indefinite sums.
+%-% \HDindex{summation!indefinite}{ugxProblemSeriesBernoulliPage}{8.9.8.}{Example: Bernoulli Polynomials and Sums of Powers}
+
+\labelSpace{3pc}
+\xtc{
+You can compute the sum of the first
+ten fourth powers by evaluating this.
+This creates a list whose entries are
+\texht{$m^4$}{\axiom{m**4}} as \texht{$m$}{\axiom{m}} ranges from 1
+to 10, and then computes the sum of the entries of that list.
+}{
+\spadpaste{reduce(+,[m**4 for m in 1..10])}
+}
+\xtc{
+You can also compute a formula for the sum of the first
+\texht{$k$}{\axiom{k}} fourth powers, where \texht{$k$}{\axiom{k}} is an
+unspecified positive integer.
+}{
+\spadpaste{sum4 := sum(m**4, m = 1..k) \bound{sum4}}
+}
+\xtc{
+This formula is valid for any positive integer \texht{$k$}{\axiom{k}}.
+For instance, if we replace \texht{$k$}{\axiom{k}} by 10,
+%-% \HDindex{summation!definite}{ugxProblemSeriesBernoulliPage}{8.9.8.}{Example: Bernoulli Polynomials and Sums of Powers}
+we obtain the number we computed earlier.
+}{
+\spadpaste{eval(sum4, k = 10) \free{sum4}}
+}
+
+You can compute a formula for the sum of the first
+\texht{$k$}{\axiom{k}} \eth{\axiom{n}} powers in a similar fashion.
+Just replace the \axiom{4} in the definition of \userfun{sum4} by
+any expression not involving \texht{$k$}{\axiom{k}}.
+\Language{} computes these formulas using Bernoulli polynomials;
+%-% \HDindex{Bernoulli!polynomial}{ugxProblemSeriesBernoulliPage}{8.9.8.}{Example: Bernoulli Polynomials and Sums of Powers}
+we
+%-% \HDindex{polynomial!Bernoulli}{ugxProblemSeriesBernoulliPage}{8.9.8.}{Example: Bernoulli Polynomials and Sums of Powers}
+use the rest of this section to describe this method.
+
+%
+\xtc{
+First consider this function of \axiom{t} and \axiom{x}.
+}{
+\spadpaste{f := t*exp(x*t) / (exp(t) - 1) \bound{f}}
+}
+\noOutputXtc{
+Since the expressions involved get quite large, we tell
+\Language{} to show us only terms of degree up to \axiom{5.}
+}{
+\spadpaste{)set streams calculate 5 \bound{set}}
+}
+%-% \HDsyscmdindex{set streams calculate}{ugxProblemSeriesBernoulliPage}{8.9.8.}{Example: Bernoulli Polynomials and Sums of Powers}
+%
+%
+\xtc{
+If we look at the Taylor expansion of \axiom{f(x, t)} about \axiom{t = 0,}
+we see that the coefficients of the powers of \axiom{t} are polynomials
+in \axiom{x}.
+}{
+\spadpaste{ff := taylor(f,t = 0) \free{f set} \bound{ff}}
+}
+%
+In fact, the \eth{\axiom{n}} coefficient in this series is essentially
+the \eth{\axiom{n}} Bernoulli polynomial:
+the \eth{\axiom{n}} coefficient of the series is
+\texht{${1 \over {n!}} B_n(x)$}{\axiom{1/n! * Bn(x)}}, where
+\texht{$B_n(x)$}{\axiom{Bn(x)}}
+is the \eth{\axiom{n}} Bernoulli polynomial.
+Thus, to obtain the \eth{\axiom{n}} Bernoulli polynomial, we multiply
+the \eth{\axiom{n}} coefficient
+of the series \axiom{ff} by \axiom{n!}.
+%
+\xtc{
+For example, the sixth Bernoulli polynomial is this.
+}{
+\spadpaste{factorial(6) * coefficient(ff,6) \free{ff}}
+}
+%
+\xtc{
+We derive some properties of the function \axiom{f(x,t)}.
+First we compute \axiom{f(x + 1,t) - f(x,t)}.
+}{
+\spadpaste{g := eval(f, x = x + 1) - f \bound{g} \free{f}}
+}
+%
+\xtc{
+If we normalize \axiom{g}, we see that it has a particularly simple form.
+}{
+\spadpaste{normalize(g) \free{g}}
+}
+%
+From this it follows that the \eth{\axiom{n}}
+coefficient in the Taylor expansion of
+\axiom{g(x,t)} at \axiom{t = 0} is
+\texht{${1\over{(n-1)\:!}}\:x^{n-1}$}{\axiom{1/(n-1)! * x**(n-1)}}.
+\xtc{
+If you want to check this, evaluate the next expression.
+}{
+\spadpaste{taylor(g,t = 0) \free{g}}
+}
+%
+However, since \axiom{g(x,t) = f(x+1,t)-f(x,t)}, it follows that the
+\eth{\axiom{n}} coefficient is
+\texht{${1 \over {n!}}\:(B_n(x+1)-B_n(x))$}{\axiom{1/n! * (Bn(x + 1) -
+Bn(x))}}.
+Equating coefficients, we see that
+\texht{${1\over{(n-1)\:!}}\:x^{n-1} = {1\over{n!}}\:(B_n(x + 1) -
+B_n(x))$}{\axiom{1/(n-1)! * x**(n-1) = 1/n! * (Bn(x + 1) - Bn(x))}}
+and, therefore,
+\texht{$x^{n-1} = {1\over{n}}\:(B_n(x + 1) -
+B_n(x))$}{\axiom{x**(n-1) = 1/n * (Bn(x + 1) - Bn(x))}}.
+%
+Let's apply this formula repeatedly, letting \axiom{x} vary between two
+integers \axiom{a} and \axiom{b}, with \axiom{a < b}:
+%
+\begin{verbatim}
+  a**(n-1)       = 1/n * (Bn(a + 1) - Bn(a))
+  (a + 1)**(n-1) = 1/n * (Bn(a + 2) - Bn(a + 1))
+  (a + 2)**(n-1) = 1/n * (Bn(a + 3) - Bn(a + 2))
+     .
+     .
+     .
+  (b - 1)**(n-1) = 1/n * (Bn(b) - Bn(b - 1))
+  b**(n-1)       = 1/n * (Bn(b + 1) - Bn(b))
+\end{verbatim}
+
+When we add these equations we find that
+the sum of the left-hand sides is
+\texht{$\sum_{m=a}^{b} m^{n-1},$}{\axiom{sum(m = a..b,m ** (n-1))},}%
+the sum of the
+\texht{$(n-1)^{\hbox{\small\rm st}}$}{\axiom{(n-1)}-st}
+powers from \axiom{a} to \axiom{b}.
+The sum of the right-hand sides is a ``telescoping series.''
+After cancellation, the sum is simply
+\texht{${1\over{n}}\:(B_n(b + 1) -
+B_n(a))$}{\axiom{1/n * (Bn(b + 1) - Bn(a))}}.
+
+Replacing \axiom{n} by \axiom{n + 1}, we have shown that
+\centerline{{\axiom{sum(m = a..b,m ** n) = 1/(n + 1) * (B<n+1>(b + 1) - B<n+1>(a))}.}}
+
+Let's use this to obtain the formula for the sum of fourth powers.
+\xtc{
+First we obtain the Bernoulli polynomial \texht{$B_5$}{\axiom{B5}}.
+}{
+\spadpaste{B5 := factorial(5) * coefficient(ff,5) \free{ff} \bound{B5}}
+}
+%
+\xtc{
+To find the sum of the first \texht{$k$}{\axiom{k}} 4th powers,
+we multiply \axiom{1/5} by
+\texht{$B_5(k+1) - B_5(1)$}{\axiom{B5(k + 1) - B5(1)}}.
+}{
+\spadpaste{1/5 * (eval(B5, x = k + 1) - eval(B5, x = 1)) \free{B5}}
+}
+%
+\xtc{
+This is the same formula that we obtained via \axiom{sum(m**4, m = 1..k)}.
+}{
+\spadpaste{sum4 \free{sum4}}
+}
+
+At this point you may want to do the same computation, but with an
+exponent other than \axiom{4.}
+For example, you might try to find a formula for the sum of the
+first \texht{$k$}{\axiom{k}} 20th powers.
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugProblemDEQTitle}{Solution of Differential Equations}
+\newcommand{\ugProblemDEQNumber}{8.10.}
+%
+% =====================================================================
+\begin{page}{ugProblemDEQPage}{8.10. Solution of Differential Equations}
+% =====================================================================
+\beginscroll
+%
+In this section we discuss \Language{}'s facilities for
+%-% \HDindex{equation!differential!solving}{ugProblemDEQPage}{8.10.}{Solution of Differential Equations}
+solving
+%-% \HDindex{differential equation}{ugProblemDEQPage}{8.10.}{Solution of Differential Equations}
+differential equations in closed-form and in series.
+
+\Language{} provides facilities for closed-form solution of
+%-% \HDindex{equation!differential!solving in closed-form}{ugProblemDEQPage}{8.10.}{Solution of Differential Equations}
+single differential equations of the following kinds:
+\indent{4}
+\beginitems
+\item[-] linear ordinary differential equations, and
+\item[-] non-linear first order ordinary differential equations
+when integrating factors can be found just by integration.
+\enditems
+\indent{0}
+
+For a discussion of the solution of systems of linear and polynomial
+equations, see
+\downlink{``\ugProblemLinPolEqnTitle''}{ugProblemLinPolEqnPage} in Section \ugProblemLinPolEqnNumber\ignore{ugProblemLinPolEqn}.
+
+\beginmenu
+    \menudownlink{{8.10.1. Closed-Form Solutions of Linear Differential Equations}}{ugxProblemLDEQClosedPage}
+    \menudownlink{{8.10.2. Closed-Form Solutions of Non-Linear Differential Equations}}{ugxProblemNLDEQClosedPage}
+    \menudownlink{{8.10.3. Power Series Solutions of Differential Equations}}{ugxProblemDEQSeriesPage}
+\endmenu
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugxProblemLDEQClosedTitle}{Closed-Form Solutions of Linear Differential Equations}
+\newcommand{\ugxProblemLDEQClosedNumber}{8.10.1.}
+%
+% =====================================================================
+\begin{page}{ugxProblemLDEQClosedPage}{8.10.1. Closed-Form Solutions of Linear Differential Equations}
+% =====================================================================
+\beginscroll
+
+A {\it differential equation} is an equation involving an unknown {\it
+function} and one or more of its derivatives.
+%-% \HDindex{differential equation}{ugxProblemLDEQClosedPage}{8.10.1.}{Closed-Form Solutions of Linear Differential Equations}
+The equation is called {\it ordinary} if derivatives with respect to
+%-% \HDindex{equation!differential}{ugxProblemLDEQClosedPage}{8.10.1.}{Closed-Form Solutions of Linear Differential Equations}
+only one dependent variable appear in the equation (it is called {\it
+partial} otherwise).
+The package \axiomType{ElementaryFunctionODESolver} provides the
+top-level operation \spadfun {solve} for finding closed-form solutions
+of ordinary differential equations.
+%-% \HDexptypeindex{ElementaryFunctionODESolver}{ugxProblemLDEQClosedPage}{8.10.1.}{Closed-Form Solutions of Linear Differential Equations}
+
+To solve a differential equation, you must first create an operator for
+%-% \HDindex{operator}{ugxProblemLDEQClosedPage}{8.10.1.}{Closed-Form Solutions of Linear Differential Equations}
+the unknown function.
+%
+\xtc{
+We let \axiom{y} be the unknown function in terms of \axiom{x}.
+}{
+\spadpaste{y := operator 'y \bound{y}}
+}
+%
+You then type the equation using \axiomFun{D} to create the
+derivatives of the unknown function \axiom{y(x)} where \axiom{x} is any
+symbol you choose (the so-called {\it dependent variable}).
+%
+\xtc{
+This is how you enter
+the equation \axiom{y'' + y' + y = 0}.
+}{
+\spadpaste{deq := D(y x, x, 2) + D(y x, x) + y x = 0\bound{e1}\free{y}}
+}
+%
+The simplest way to invoke the \axiomFun{solve} command is with three
+arguments.
+\begin{items}
+\item the differential equation,
+\item the operator representing the unknown function,
+\item the dependent variable.
+\end{items}
+%
+\xtc{
+So, to solve the above equation, we enter this.
+}{
+\spadpaste{solve(deq, y, x) \free{e1}\free{y}}
+}
+%
+Since linear ordinary differential equations have infinitely many
+solutions, \axiomFun{solve} returns a {\it particular solution}
+\texht{$f_p$}{\axiom{f_p}}
+and a basis
+\texht{$f_1,\dots,f_n$}{\axiom{f1,...,fn}}
+for the solutions of the corresponding homogenuous equation.
+Any expression of the form
+\texht{$f_p + c_1 f_1 + \dots c_n f_n$}{\axiom{fp + c1 f1 + ... + cn fn}}
+where the \texht{$c_i$}{\axiom{ci}} do not involve
+the dependent variable is also a solution.
+This is similar to what you get when you solve systems of linear
+algebraic equations.
+
+A way to select a unique solution is to specify {\it initial
+conditions}: choose a value \axiom{a} for the dependent variable
+and specify the values of the unknown function and its derivatives
+at \axiom{a}.
+If the number of initial conditions is equal to the order of the
+equation, then the solution is unique (if it exists in closed
+form!) and \axiomFun{solve} tries to find it.
+To specify initial conditions to \axiomFun{solve}, use an
+\axiomType{Equation} of the form \axiom{x = a} for the third
+parameter instead of the dependent variable, and add a fourth
+parameter consisting of the list of values \axiom{y(a), y'(a), ...}.
+
+\xtc{
+To find the solution of \axiom{y'' + y = 0} satisfying \axiom{y(0) = y'(0) = 1},
+do this.
+}{
+\spadpaste{deq := D(y x, x, 2) + y x \bound{e2}\free{y}}
+}
+\xtc{
+You can omit the \axiom{= 0} when you enter the equation to be solved.
+}{
+\spadpaste{solve(deq, y, x = 0, [1, 1]) \free{e2}\free{y}}
+}
+%
+
+\Language{} is not limited to linear differential equations with
+constant coefficients.
+It can also find solutions when the coefficients are rational or
+algebraic functions of the dependent variable.
+Furthermore, \Language{} is not limited by the order of the equation.
+%
+\xtc{
+\Language{} can solve the following third order equations with
+polynomial coefficients.
+}{
+\spadpaste{deq := x**3 * D(y x, x, 3) + x**2 * D(y x, x, 2) - 2 * x * D(y x, x) + 2 * y x = 2 * x**4 \bound{e3}\free{y}}
+}
+\xtc{
+}{
+\spadpaste{solve(deq, y, x) \free{e3}\free{y}}
+}
+%
+%
+\xtc{
+Here we are solving a homogeneous equation.
+}{
+\spadpaste{deq := (x**9+x**3) * D(y x, x, 3) + 18 * x**8 * D(y x, x, 2) - 90 * x * D(y x, x) - 30 * (11 * x**6 - 3) * y x \bound{e4}\free{y}}
+}
+\xtc{
+}{
+\spadpaste{solve(deq, y, x) \free{e4}\free{y}}
+}
+%
+On the other hand, and in contrast with the operation
+\axiomFun{integrate}, it can happen that \Language{} finds no solution
+and that some closed-form solution still exists.
+While it is mathematically complicated to describe exactly when the
+solutions are guaranteed to be found, the following statements are
+correct and form good guidelines for linear ordinary differential
+equations:
+\begin{items}
+\item If the coefficients are constants, \Language{} finds a complete basis
+of solutions (i,e, all solutions).
+\item If the coefficients are rational functions in the dependent variable,
+\Language{} at least finds all solutions that do not involve algebraic
+functions.
+\end{items}
+%
+Note that this last statement does not mean that \Language{} does not find
+the solutions that are algebraic functions.
+It means that it is not
+guaranteed that the algebraic function solutions will be found.
+%
+\xtc{
+This is an example where all the algebraic solutions are found.
+}{
+\spadpaste{deq := (x**2 + 1) * D(y x, x, 2) + 3 * x * D(y x, x) + y x = 0 \bound{e5}\free{y}}
+}
+\xtc{
+}{
+\spadpaste{solve(deq, y, x) \free{e5}\free{y}}
+}
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugxProblemNLDEQClosedTitle}{Closed-Form Solutions of Non-Linear Differential Equations}
+\newcommand{\ugxProblemNLDEQClosedNumber}{8.10.2.}
+%
+% =====================================================================
+\begin{page}{ugxProblemNLDEQClosedPage}{8.10.2. Closed-Form Solutions of Non-Linear Differential Equations}
+% =====================================================================
+\beginscroll
+
+This is an example that shows how to solve a non-linear
+first order ordinary differential equation manually when an integrating
+factor can be found just by integration.
+At the end, we show you how to solve it directly.
+
+Let's solve the differential equation \axiom{y' = y / (x + y log y)}.
+%
+\xtc{
+Using the notation
+\axiom{m(x, y) + n(x, y) y' = 0},
+we have \axiom{m = -y} and \axiom{n = x + y log y}.
+}{
+\spadpaste{m := -y \bound{m}}
+}
+\xtc{
+}{
+\spadpaste{n := x + y * log y \bound{n}}
+}
+%
+\xtc{
+We first check for exactness, that is, does \axiom{dm/dy = dn/dx}?
+}{
+\spadpaste{D(m, y) - D(n, x) \free{m n}}
+}
+%
+This is not zero, so the equation is not exact.
+Therefore we must look for
+an integrating factor: a function \axiom{mu(x,y)} such that
+\axiom{d(mu m)/dy = d(mu n)/dx}.
+Normally, we first search for \axiom{mu(x,y)} depending only on
+\axiom{x} or only on \axiom{y}.
+%
+\xtc{
+Let's search for such a \axiom{mu(x)} first.
+}{
+\spadpaste{mu := operator 'mu \bound{mu}}
+}
+\xtc{
+}{
+\spadpaste{a := D(mu(x) * m, y) - D(mu(x) * n, x) \bound{a}\free{m n mu}}
+}
+%
+%
+\xtc{
+If the above is zero for a function
+\axiom{mu} that does {\it not} depend on \axiom{y}, then
+\axiom{mu(x)} is an integrating factor.
+}{
+\spadpaste{solve(a = 0, mu, x) \free{mu a}}
+}
+%
+The solution depends on \axiom{y}, so there is no integrating
+factor that depends on \axiom{x} only.
+%
+\xtc{
+Let's look for one that depends on \axiom{y} only.
+}{
+\spadpaste{b := D(mu(y) * m, y) - D(mu(y) * n, x) \bound{b}\free{mu m}}
+}
+\xtc{
+}{
+\spadpaste{sb := solve(b = 0, mu, y) \free{mu b}\bound{sb}}
+}
+\noindent
+We've found one!
+%
+\xtc{
+The above \axiom{mu(y)} is an integrating factor.
+We must multiply our initial equation
+(that is, \axiom{m} and \axiom{n}) by the integrating factor.
+}{
+\spadpaste{intFactor := sb.basis.1 \bound{intFactor}\free{sb}}
+}
+\xtc{
+}{
+\spadpaste{m := intFactor * m \bound{m1}\free{m intFactor}}
+}
+\xtc{
+}{
+\spadpaste{n := intFactor * n \bound{n1}\free{n intFactor}}
+}
+%
+\xtc{
+Let's check for exactness.
+}{
+\spadpaste{D(m, y) - D(n, x) \free{m1 n1}}
+}
+%
+We must solve the exact equation, that is, find a function
+\axiom{s(x,y)} such that
+\axiom{ds/dx = m}  and \axiom{ds/dy = n}.
+%
+\xtc{
+We start by writing \axiom{s(x, y) = h(y) + integrate(m, x)}
+where \axiom{h(y)} is an unknown function of \axiom{y}.
+This guarantees that \axiom{ds/dx = m}.
+}{
+\spadpaste{h := operator 'h \bound{h}}
+}
+\xtc{
+}{
+\spadpaste{sol := h y + integrate(m, x) \bound{sol}\free{h m1}}
+}
+%
+%
+\xtc{
+All we want is to find \axiom{h(y)} such that
+\axiom{ds/dy = n}.
+}{
+\spadpaste{dsol := D(sol, y) \free{sol}\bound{dsol}}
+}
+\xtc{
+}{
+\spadpaste{nsol := solve(dsol = n, h, y) \free{dsol n1 h}\bound{nsol}}
+}
+%
+\xtc{
+The above particular solution is the \axiom{h(y)} we want, so we just replace
+\axiom{h(y)} by it in the implicit solution.
+}{
+\spadpaste{eval(sol, h y = nsol.particular) \free{sol h nsol}}
+}
+%
+A first integral of the initial equation is obtained by setting
+this result equal to an arbitrary constant.
+
+Now that we've seen how to solve the equation ``by hand,''
+we show you how to do it with the \axiomFun{solve} operation.
+\xtc{
+First define \axiom{y} to be an operator.
+}{
+\spadpaste{y := operator 'y \bound{y}}
+}
+\xtc{
+Next we create the differential equation.
+}{
+\spadpaste{deq := D(y x, x) = y(x) / (x + y(x) * log y x) \bound{deqi}\free{y}}
+}
+\xtc{
+Finally, we solve it.
+}{
+\spadpaste{solve(deq, y, x) \free{deqi y}}
+}
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugxProblemDEQSeriesTitle}{Power Series Solutions of Differential Equations}
+\newcommand{\ugxProblemDEQSeriesNumber}{8.10.3.}
+%
+% =====================================================================
+\begin{page}{ugxProblemDEQSeriesPage}{8.10.3. Power Series Solutions of Differential Equations}
+% =====================================================================
+\beginscroll
+
+The command to solve differential equations in power
+%-% \HDindex{equation!differential!solving in power series}{ugxProblemDEQSeriesPage}{8.10.3.}{Power Series Solutions of Differential Equations}
+series
+%-% \HDindex{power series}{ugxProblemDEQSeriesPage}{8.10.3.}{Power Series Solutions of Differential Equations}
+around
+%-% \HDindex{series!power}{ugxProblemDEQSeriesPage}{8.10.3.}{Power Series Solutions of Differential Equations}
+a particular initial point with specific initial conditions is called
+\axiomFun{seriesSolve}.
+It can take a variety of parameters, so we illustrate
+its use with some examples.
+
+\labelSpace{1pc}
+\noOutputXtc{
+Since the coefficients of some solutions
+are quite large, we reset the default to compute only seven terms.
+}{
+\spadpaste{)set streams calculate 7 \bound{c7}}
+}
+
+You can solve a single nonlinear equation of any order. For example,
+we solve  \axiom{y''' = sin(y'') * exp(y) + cos(x)}
+subject to \axiom{y(0) = 1, y'(0) = 0, y''(0) = 0}.
+
+\xtc{
+We first tell \Language{}
+that the symbol \axiom{'y} denotes a new operator.
+}{
+\spadpaste{y := operator 'y \bound{y}}
+}
+\xtc{
+Enter the differential equation using \axiom{y} like any system
+function.
+}{
+\spadpaste{eq := D(y(x), x, 3) - sin(D(y(x), x, 2))*exp(y(x)) = cos(x)\bound{eq}\free{y}}
+}
+%
+\xtc{
+Solve it around \axiom{x = 0} with the initial conditions
+\axiom{y(0) = 1, y'(0) = y''(0) = 0}.
+}{
+\spadpaste{seriesSolve(eq, y, x = 0, [1, 0, 0])\free{y}\free{eq}\free{c7}}
+}
+
+You can also solve a system of nonlinear first order equations.
+For example, we solve a system that has \axiom{tan(t)} and
+\axiom{sec(t)} as solutions.
+
+\xtc{
+We tell \Language{} that \axiom{x} is also an operator.
+}{
+\spadpaste{x := operator 'x\bound{x}}
+}
+\xtc{
+Enter the two equations forming our system.
+}{
+\spadpaste{eq1 := D(x(t), t) = 1 + x(t)**2\free{x}\free{y}\bound{eq1}}
+}
+%
+\xtc{
+}{
+\spadpaste{eq2 := D(y(t), t) = x(t) * y(t)\free{x}\free{y}\bound{eq2}}
+}
+%
+\xtc{
+Solve the system around \axiom{t = 0} with the initial conditions
+\axiom{x(0) = 0} and \axiom{y(0) = 1}.
+Notice that since we give the unknowns in the
+order \axiom{[x, y]}, the answer is a list of two series in the order
+\axiom{[series for x(t), series for y(t)]}.
+}{
+\spadpaste{seriesSolve([eq2, eq1], [x, y], t = 0, [y(0) = 1, x(0) = 0])\free{x}\free{y}\free{eq1}\free{eq2}\free{c7}}
+}
+\noindent
+The order in which we give the
+equations and the initial conditions has no effect on the order of
+the solution.
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugProblemFiniteTitle}{Finite Fields}
+\newcommand{\ugProblemFiniteNumber}{8.11.}
+%
+% =====================================================================
+\begin{page}{ugProblemFinitePage}{8.11. Finite Fields}
+% =====================================================================
+\beginscroll
+%
+
+A {\it finite field} (also called a {\it Galois field}) is a
+finite algebraic structure where one can add, multiply and divide
+under the same laws (for example, commutativity, associativity or
+distributivity) as apply to the rational, real or complex numbers.
+Unlike those three fields, for any finite field there exists a
+positive prime integer \smath{p}, called the
+\axiomFun{characteristic}, such that \texht{$p \: x = 0$}{\axiom{p *
+x = 0}} for any element \smath{x} in the finite field.
+In fact, the number of elements in a finite field is a power of
+the characteristic and for each prime \smath{p} and positive
+integer \smath{n} there exists exactly one finite field with
+\texht{$p^n$}{\axiom{p**n}} elements, up to
+isomorphism.\footnote{For more information about the algebraic
+structure and properties of finite fields, see, for example, S.
+Lang, {\it Algebra}, Second Edition, New York: Addison-Wesley
+Publishing Company, Inc., 1984, ISBN 0 201 05487 6; or R.
+Lidl, H.
+Niederreiter, {\it Finite Fields}, Encyclopedia of Mathematics and
+Its Applications, Vol.
+20, Cambridge: Cambridge Univ.
+Press, 1983, ISBN 0 521 30240 4.}
+
+When \axiom{n = 1,} the field has \smath{p} elements and is
+called a {\it prime field}, discussed in \texht{the next
+section}{\downlink{``\ugxProblemFinitePrimeTitle''}{ugxProblemFinitePrimePage} in Section \ugxProblemFinitePrimeNumber\ignore{ugxProblemFinitePrime}}.
+There are several ways of implementing extensions of finite
+fields, and \Language{} provides quite a bit of freedom to allow
+you to choose the one that is best for your application.
+Moreover, we provide operations for converting among the different
+representations of extensions and different extensions of a single
+field.
+Finally, note that you usually need to package-call operations
+from finite fields if the operations do not take as an argument an
+object of the field.
+See \downlink{``\ugTypesPkgCallTitle''}{ugTypesPkgCallPage} in Section \ugTypesPkgCallNumber\ignore{ugTypesPkgCall} for more information on
+package-calling.
+
+\beginmenu
+    \menudownlink{{8.11.1. Modular Arithmetic and Prime Fields}}{ugxProblemFinitePrimePage}
+    \menudownlink{{8.11.2. Extensions of Finite Fields}}{ugxProblemFiniteExtensionFinitePage}
+    \menudownlink{{8.11.3. Irreducible Modulus Polynomial Representations}}{ugxProblemFiniteModulusPage}
+    \menudownlink{{8.11.4. Cyclic Group Representations}}{ugxProblemFiniteCyclicPage}
+    \menudownlink{{8.11.5. Normal Basis Representations}}{ugxProblemFiniteNormalPage}
+    \menudownlink{{8.11.6. Conversion Operations for Finite Fields}}{ugxProblemFiniteConversionPage}
+    \menudownlink{{8.11.7. Utility Operations for Finite Fields}}{ugxProblemFiniteUtilityPage}
+\endmenu
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugxProblemFinitePrimeTitle}{Modular Arithmetic and Prime Fields}
+\newcommand{\ugxProblemFinitePrimeNumber}{8.11.1.}
+%
+% =====================================================================
+\begin{page}{ugxProblemFinitePrimePage}{8.11.1. Modular Arithmetic and Prime Fields}
+% =====================================================================
+\beginscroll
+%-% \HDindex{finite field}{ugxProblemFinitePrimePage}{8.11.1.}{Modular Arithmetic and Prime Fields}
+%-% \HDindex{Galois!field}{ugxProblemFinitePrimePage}{8.11.1.}{Modular Arithmetic and Prime Fields}
+%-% \HDindex{field!finite!prime}{ugxProblemFinitePrimePage}{8.11.1.}{Modular Arithmetic and Prime Fields}
+%-% \HDindex{field!prime}{ugxProblemFinitePrimePage}{8.11.1.}{Modular Arithmetic and Prime Fields}
+%-% \HDindex{field!Galois}{ugxProblemFinitePrimePage}{8.11.1.}{Modular Arithmetic and Prime Fields}
+%-% \HDindex{prime field}{ugxProblemFinitePrimePage}{8.11.1.}{Modular Arithmetic and Prime Fields}
+%-% \HDindex{modular arithmetic}{ugxProblemFinitePrimePage}{8.11.1.}{Modular Arithmetic and Prime Fields}
+%-% \HDindex{arithmetic!modular}{ugxProblemFinitePrimePage}{8.11.1.}{Modular Arithmetic and Prime Fields}
+
+Let \smath{n} be a positive integer.
+It is well known that you can get the same result if you perform
+addition, subtraction or multiplication of integers and then take
+the remainder on dividing by \smath{n} as if
+you had first done such remaindering on the operands, performed the
+arithmetic and then (if necessary) done remaindering again.
+This allows us to speak of arithmetic
+{\it modulo} \smath{n} or, more simply
+{\it mod} \smath{n}.
+\xtc{
+In \Language{}, you use \axiomType{IntegerMod} to do such arithmetic.
+}{
+\spadpaste{(a,b) : IntegerMod 12 \bound{abdec}}
+}
+\xtc{
+}{
+\spadpaste{(a, b) := (16, 7) \free{abdec}\bound{a b}}
+}
+\xtc{
+}{
+\spadpaste{[a - b, a * b] \free{a b}}
+}
+\xtc{
+If \smath{n} is not prime, there is only a limited notion of
+reciprocals and division.
+}{
+\spadpaste{a / b \free{a b}}
+}
+\xtc{
+}{
+\spadpaste{recip a \free{a}}
+}
+\xtc{
+Here \axiom{7} and \axiom{12} are relatively prime, so \axiom{7}
+has a multiplicative inverse modulo \axiom{12}.
+}{
+\spadpaste{recip b \free{b}}
+}
+
+If we take \smath{n} to be a prime number \smath{p},
+then taking inverses and, therefore, division are generally defined.
+\xtc{
+Use \axiomType{PrimeField} instead of \axiomType{IntegerMod}
+for \smath{n} prime.
+}{
+\spadpaste{c : PrimeField 11 := 8 \bound{c}}
+}
+\xtc{
+}{
+\spadpaste{inv c \free{c}}
+}
+\xtc{
+You can also use \axiom{1/c} and \axiom{c**(-1)} for the inverse of
+\smath{c}.
+}{
+\spadpaste{9/c \free{c}}
+}
+
+\axiomType{PrimeField} (abbreviation \axiomType{PF}) checks if its
+argument is prime when you try to use an operation from it.
+If you know the argument is prime (particularly if it is large),
+\axiomType{InnerPrimeField} (abbreviation \axiomType{IPF}) assumes
+the argument has already been verified to be prime.
+If you do use a number that is not prime, you will eventually get
+an error message, most likely a division by zero message.
+For computer science applications, the most important finite fields
+are \axiomType{PrimeField 2} and its extensions.
+
+\xtc{
+In the following examples, we work with the finite field with
+\smath{p = 101} elements.
+}{
+\spadpaste{GF101 := PF 101 \bound{GF101} }
+}
+\xtc{
+Like many domains in \Language{}, finite fields provide an operation
+for returning a random element of the domain.
+}{
+\spadpaste{x := random()\$GF101 \bound{x}\free{GF101}}
+}
+\xtc{
+}{
+\spadpaste{y : GF101 := 37 \bound{y}\free{GF101}}
+}
+\xtc{
+}{
+\spadpaste{z := x/y \bound{z}\free{x y}}
+}
+\xtc{
+}{
+\spadpaste{z * y - x \free{x y z}}
+}
+%
+\xtc{
+The element \axiom{2} is a {\it primitive element} of this field,
+%-% \HDindex{primitive element}{ugxProblemFinitePrimePage}{8.11.1.}{Modular Arithmetic and Prime Fields}
+%-% \HDindex{element!primitive}{ugxProblemFinitePrimePage}{8.11.1.}{Modular Arithmetic and Prime Fields}
+}{
+\spadpaste{pe := primitiveElement()\$GF101 \bound{pe}\free{GF101}}
+}
+%
+\xtc{
+in the sense that its powers enumerate all nonzero elements.
+}{
+\spadpaste{[pe**i for i in 0..99] \free{pe}}
+}
+%
+%
+\xtc{
+If every nonzero element is a power of a primitive element, how do you
+determine what the exponent is?
+Use
+%-% \HDindex{discrete logarithm}{ugxProblemFinitePrimePage}{8.11.1.}{Modular Arithmetic and Prime Fields}
+\axiomFun{discreteLog}.
+%-% \HDindex{logarithm!discrete}{ugxProblemFinitePrimePage}{8.11.1.}{Modular Arithmetic and Prime Fields}
+}{
+\spadpaste{ex := discreteLog(y) \bound{ex}\free{y}}
+}
+\xtc{
+}{
+\spadpaste{pe ** ex \free{ex pe}}
+}
+%
+\xtc{
+The \axiomFun{order} of a nonzero element \smath{x} is the
+smallest positive integer \smath{t} such
+\texht{$x^t = 1$}{\axiom{x**t = 1}}.
+}{
+\spadpaste{order y \free{y}}
+}
+\xtc{
+The order of a primitive element is the defining \smath{p-1}.
+}{
+\spadpaste{order pe \free{pe}}
+}
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugxProblemFiniteExtensionFiniteTitle}{Extensions of Finite Fields}
+\newcommand{\ugxProblemFiniteExtensionFiniteNumber}{8.11.2.}
+%
+% =====================================================================
+\begin{page}{ugxProblemFiniteExtensionFinitePage}{8.11.2. Extensions of Finite Fields}
+% =====================================================================
+\beginscroll
+%-% \HDindex{finite field}{ugxProblemFiniteExtensionFinitePage}{8.11.2.}{Extensions of Finite Fields}
+%-% \HDindex{field!finite!extension of}{ugxProblemFiniteExtensionFinitePage}{8.11.2.}{Extensions of Finite Fields}
+
+When you want to work with an extension of a finite field in \Language{},
+you have three choices to make:
+\indent{4}
+\beginitems
+\item[1. ] Do you want to generate an extension of the prime field
+(for example, \axiomType{PrimeField 2}) or an extension of a given field?
+%
+\item[2. ] Do you want to use a representation that is particularly
+efficient for multiplication, exponentiation and addition but
+uses a lot of computer memory (a representation that models the cyclic
+group structure of the multiplicative group of the field extension
+and uses a Zech logarithm table),  one that
+%-% \HDindex{Zech logarithm}{ugxProblemFiniteExtensionFinitePage}{8.11.2.}{Extensions of Finite Fields}
+uses a normal basis for the vector space structure of the field
+extension, or one that performs arithmetic modulo an irreducible
+polynomial?
+The cyclic group representation is only usable up to ``medium''
+(relative to your machine's performance)
+sized fields.
+If the field is large and the normal basis is relatively simple,
+the normal basis representation is more efficient for exponentiation than
+the irreducible polynomial representation.
+%
+\item[3. ] Do you want to provide a polynomial explicitly, a root of which
+``generates'' the extension in one of the three senses in (2),
+or do you wish to have the polynomial generated for you?
+\enditems
+\indent{0}
+
+This illustrates one of the most important features of \Language{}:
+you can choose exactly the right data-type and representation to
+suit your application best.
+
+We first tell you what domain constructors to use for each case
+above, and then give some examples.
+
+\texht{\hangafter=1\hangindent=2pc}{\noindent}
+Constructors that automatically generate extensions of the prime field:
+\newline
+\axiomType{FiniteField} \newline
+\axiomType{FiniteFieldCyclicGroup} \newline
+\axiomType{FiniteFieldNormalBasis}
+
+\texht{\hangafter=1\hangindent=2pc}{\noindent}
+Constructors that generate extensions of an arbitrary field:
+\newline
+\axiomType{FiniteFieldExtension} \newline
+\axiomType{FiniteFieldExtensionByPolynomial} \newline
+\axiomType{FiniteFieldCyclicGroupExtension} \newline
+\axiomType{FiniteFieldCyclicGroupExtensionByPolynomial} \newline
+\axiomType{FiniteFieldNormalBasisExtension} \newline
+\axiomType{FiniteFieldNormalBasisExtensionByPolynomial}
+
+\texht{\hangafter=1\hangindent=2pc}{\noindent}
+Constructors that use a cyclic group representation:
+\newline
+\axiomType{FiniteFieldCyclicGroup} \newline
+\axiomType{FiniteFieldCyclicGroupExtension} \newline
+\axiomType{FiniteFieldCyclicGroupExtensionByPolynomial}
+
+\texht{\hangafter=1\hangindent=2pc}{\noindent}
+Constructors that use a normal basis representation:
+\newline
+\axiomType{FiniteFieldNormalBasis} \newline
+\axiomType{FiniteFieldNormalBasisExtension} \newline
+\axiomType{FiniteFieldNormalBasisExtensionByPolynomial}
+
+\texht{\hangafter=1\hangindent=2pc}{\noindent}
+Constructors that use an irreducible modulus polynomial representation:
+\newline
+\axiomType{FiniteField} \newline
+\axiomType{FiniteFieldExtension} \newline
+\axiomType{FiniteFieldExtensionByPolynomial}
+
+\texht{\hangafter=1\hangindent=2pc}{\noindent}
+Constructors that generate a polynomial for you:
+\newline
+\axiomType{FiniteField} \newline
+\axiomType{FiniteFieldExtension} \newline
+\axiomType{FiniteFieldCyclicGroup} \newline
+\axiomType{FiniteFieldCyclicGroupExtension} \newline
+\axiomType{FiniteFieldNormalBasis} \newline
+\axiomType{FiniteFieldNormalBasisExtension}
+
+\texht{\hangafter=1\hangindent=2pc}{\noindent}
+Constructors for which you provide a polynomial:
+\newline
+\axiomType{FiniteFieldExtensionByPolynomial} \newline
+\axiomType{FiniteFieldCyclicGroupExtensionByPolynomial} \newline
+\axiomType{FiniteFieldNormalBasisExtensionByPolynomial}
+
+These constructors are discussed in the following sections where
+we collect together descriptions of extension fields that have the
+same underlying representation.\footnote{For
+more information on the implementation aspects of finite
+fields, see
+J. Grabmeier, A. Scheerhorn, {\it Finite Fields in AXIOM,}
+Technical Report, IBM Heidelberg Scientific Center, 1992.}
+
+If you don't really care about all this detail, just use
+\axiomType{FiniteField}.
+As your knowledge of your application and its \Language{} implementation
+grows, you can come back and choose an alternative constructor that
+may improve the efficiency of your code.
+Note that the exported operations are almost the same for all constructors
+of finite field extensions and include the operations exported by
+\axiomType{PrimeField}.
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugxProblemFiniteModulusTitle}{Irreducible Modulus Polynomial Representations}
+\newcommand{\ugxProblemFiniteModulusNumber}{8.11.3.}
+%
+% =====================================================================
+\begin{page}{ugxProblemFiniteModulusPage}{8.11.3. Irreducible Modulus Polynomial Representations}
+% =====================================================================
+\beginscroll
+
+All finite field extension constructors discussed in this
+%-% \HDindex{finite field}{ugxProblemFiniteModulusPage}{8.11.3.}{Irreducible Modulus Polynomial Representations}
+section
+%-% \HDindex{field!finite!extension of}{ugxProblemFiniteModulusPage}{8.11.3.}{Irreducible Modulus Polynomial Representations}
+use a representation that performs arithmetic with univariate
+(one-variable) polynomials modulo an irreducible polynomial.
+This polynomial may be given explicitly by you or automatically
+generated.
+The ground field may be the prime field or one you specify.
+See \downlink{``\ugxProblemFiniteExtensionFiniteTitle''}{ugxProblemFiniteExtensionFinitePage} in Section \ugxProblemFiniteExtensionFiniteNumber\ignore{ugxProblemFiniteExtensionFinite} for general
+information about finite field extensions.
+
+For \axiomType{FiniteField} (abbreviation \axiomType{FF}) you provide a
+prime number \smath{p} and an extension degree \smath{n}.
+This degree can be 1.
+%
+\xtc{
+\Language{} uses the prime field \axiomType{PrimeField(p)},
+here \axiomType{PrimeField 2},
+and it chooses an irreducible polynomial of degree \smath{n},
+here 12, over the ground field.
+}{
+\spadpaste{GF4096 := FF(2,12); \bound{GF4096}}
+}
+%
+
+The objects in the generated field extension are polynomials
+of degree at most \smath{n-1} with coefficients in the
+prime field.
+The polynomial indeterminate is automatically chosen by \Language{} and
+is typically something like \axiom{\%A} or \axiom{\%D}.
+These (strange) variables are {\it only} for output display;
+there are several ways to construct elements of this field.
+
+The operation \axiomFun{index} enumerates the elements of the field
+extension and accepts as argument the integers from 1 to
+\smath{p \texht{^}{**} n}.
+%
+\xtc{
+The expression
+\axiom{index(p)} always gives the indeterminate.
+}{
+\spadpaste{a := index(2)\$GF4096 \bound{a}\free{GF4096}}
+}
+%
+%
+\xtc{
+You can build polynomials in \smath{a} and calculate in
+\axiom{GF4096}.
+}{
+\spadpaste{b := a**12 - a**5 + a \bound{b}\free{a}}
+}
+\xtc{
+}{
+\spadpaste{b ** 1000 \free{b}}
+}
+\xtc{
+}{
+\spadpaste{c := a/b \free{a b}\bound{c}}
+}
+%
+\xtc{
+Among the available operations are \axiomFun{norm} and \axiomFun{trace}.
+}{
+\spadpaste{norm c \free{c}}
+}
+\xtc{
+}{
+\spadpaste{trace c \free{c}}
+}
+%
+%
+
+Since any nonzero element is a power of a primitive element, how
+do we discover what the exponent is?
+%
+\xtc{
+The operation \axiomFun{discreteLog} calculates
+%-% \HDindex{discrete logarithm}{ugxProblemFiniteModulusPage}{8.11.3.}{Irreducible Modulus Polynomial Representations}
+the exponent and,
+%-% \HDindex{logarithm!discrete}{ugxProblemFiniteModulusPage}{8.11.3.}{Irreducible Modulus Polynomial Representations}
+if it is called with only one argument, always refers to the primitive
+element returned by \axiomFun{primitiveElement}.
+}{
+\spadpaste{dL := discreteLog a \free{a}\bound{dL}}
+}
+\xtc{
+}{
+\spadpaste{g ** dL \free{dL g}}
+}
+
+\axiomType{FiniteFieldExtension} (abbreviation \axiomType{FFX}) is
+similar to \axiomType{FiniteField} except that the
+ground-field for \axiomType{FiniteFieldExtension} is arbitrary and
+chosen by you.
+%
+\xtc{
+In case you select the prime field as ground field, there is
+essentially no difference between the constructed two finite field
+extensions.
+}{
+\spadpaste{GF16 := FF(2,4); \bound{GF16}}
+}
+\xtc{
+}{
+\spadpaste{GF4096 := FFX(GF16,3); \bound{GF4096x}\free{GF16}}
+}
+\xtc{
+}{
+\spadpaste{r := (random()\$GF4096) ** 20 \free{GF4096x}\bound{r}}
+}
+\xtc{
+}{
+\spadpaste{norm(r) \free{r}}
+}
+%
+
+\axiomType{FiniteFieldExtensionByPolynomial} (abbreviation \axiomType{FFP})
+is similar to \axiomType{FiniteField} and \axiomType{FiniteFieldExtension}
+but is more general.
+%
+\xtc{
+}{
+\spadpaste{GF4 := FF(2,2); \bound{GF4}}
+}
+\xtc{
+}{
+\spadpaste{f := nextIrreduciblePoly(random(6)\$FFPOLY(GF4))\$FFPOLY(GF4) \free{GF4}\bound{f}}
+}
+\xtc{
+For \axiomType{FFP} you choose both the
+ground field and the irreducible polynomial used in the representation.
+The degree of the extension is the degree of the polynomial.
+}{
+\spadpaste{GF4096 := FFP(GF4,f); \bound{GF4096y}\free{f GF4}}
+}
+\xtc{
+}{
+\spadpaste{discreteLog random()\$GF4096 \free{GF4096y}}
+}
+%
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugxProblemFiniteCyclicTitle}{Cyclic Group Representations}
+\newcommand{\ugxProblemFiniteCyclicNumber}{8.11.4.}
+%
+% =====================================================================
+\begin{page}{ugxProblemFiniteCyclicPage}{8.11.4. Cyclic Group Representations}
+% =====================================================================
+\beginscroll
+%-% \HDindex{finite field}{ugxProblemFiniteCyclicPage}{8.11.4.}{Cyclic Group Representations}
+%-% \HDindex{field!finite!extension of}{ugxProblemFiniteCyclicPage}{8.11.4.}{Cyclic Group Representations}
+
+In every finite field there exist elements whose powers are all the
+nonzero elements of the field.
+Such an element is called a {\it primitive element}.
+
+In \axiomType{FiniteFieldCyclicGroup} (abbreviation \axiomType{FFCG})
+%-% \HDindex{group!cyclic}{ugxProblemFiniteCyclicPage}{8.11.4.}{Cyclic Group Representations}
+the nonzero elements are represented by the
+powers of a fixed primitive
+%-% \HDindex{element!primitive}{ugxProblemFiniteCyclicPage}{8.11.4.}{Cyclic Group Representations}
+element
+%-% \HDindex{primitive element}{ugxProblemFiniteCyclicPage}{8.11.4.}{Cyclic Group Representations}
+of the field (that is, a generator of its
+cyclic multiplicative group).
+Multiplication (and hence exponentiation) using this representation is easy.
+To do addition, we consider our primitive element as the root of a primitive
+polynomial (an irreducible polynomial whose
+roots are all primitive).
+See \downlink{``\ugxProblemFiniteUtilityTitle''}{ugxProblemFiniteUtilityPage} in Section \ugxProblemFiniteUtilityNumber\ignore{ugxProblemFiniteUtility} for examples of how to
+compute such a polynomial.
+
+%
+\xtc{
+To use \axiomType{FiniteFieldCyclicGroup} you provide a prime number and an
+extension degree.
+}{
+\spadpaste{GF81 := FFCG(3,4); \bound{GF81}}
+}
+%
+%
+\xtc{
+\Language{} uses the prime field, here \axiomType{PrimeField 3}, as the
+ground field and it chooses a primitive polynomial of degree
+\smath{n}, here 4, over the prime field.
+}{
+\spadpaste{a := primitiveElement()\$GF81 \bound{a}\free{GF81}}
+}
+%
+%
+\xtc{
+You can calculate in \axiom{GF81}.
+}{
+\spadpaste{b  := a**12 - a**5 + a \bound{b}\free{a}}
+}
+%
+\xtc{
+In this representation of finite fields the discrete logarithm
+of an element can be seen directly in its output form.
+}{
+\spadpaste{b \free{b}}
+}
+\xtc{
+}{
+\spadpaste{discreteLog b \free{b}}
+}
+%
+
+\axiomType{FiniteFieldCyclicGroupExtension} (abbreviation
+\axiomType{FFCGX}) is similar to \axiomType{FiniteFieldCyclicGroup}
+except that the ground field for
+\axiomType{FiniteFieldCyclicGroupExtension} is arbitrary and chosen
+by you.
+In case you select the prime field as ground field, there is
+essentially no difference between the constructed two finite field
+extensions.
+%
+\xtc{
+}{
+\spadpaste{GF9 := FF(3,2); \bound{GF9}}
+}
+\xtc{
+}{
+\spadpaste{GF729 := FFCGX(GF9,3); \bound{GF729}\free{GF9}}
+}
+\xtc{
+}{
+\spadpaste{r := (random()\$GF729) ** 20 \free{GF729}\bound{r}}
+}
+\xtc{
+}{
+\spadpaste{trace(r) \free{r}}
+}
+%
+
+\axiomType{FiniteFieldCyclicGroupExtensionByPolynomial}
+(abbreviation \axiomType{FFCGP})
+is similar to \axiomType{FiniteFieldCyclicGroup} and
+\axiomType{FiniteFieldCyclicGroupExtension}
+but is more general.
+For \axiomType{FiniteFieldCyclicGroupExtensionByPolynomial} you choose both the
+ground field and the irreducible polynomial used in the representation.
+The degree of the extension is the degree of the polynomial.
+%
+\xtc{
+}{
+\spadpaste{GF3  := PrimeField 3; \bound{GF3}}
+}
+\xtc{
+We use a utility operation to generate an irreducible primitive
+polynomial (see \downlink{``\ugxProblemFiniteUtilityTitle''}{ugxProblemFiniteUtilityPage} in Section \ugxProblemFiniteUtilityNumber\ignore{ugxProblemFiniteUtility}).
+The polynomial has one variable that is ``anonymous'': it displays
+as a question mark.
+}{
+\spadpaste{f := createPrimitivePoly(4)\$FFPOLY(GF3) \bound{f}\free{GF3}}
+}
+\xtc{
+}{
+\spadpaste{GF81 := FFCGP(GF3,f); \bound{GF81x}\free{f GF3}}
+}
+\xtc{
+Let's look at a random element from this field.
+}{
+\spadpaste{random()\$GF81 \free{GF81x}}
+}
+%
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugxProblemFiniteNormalTitle}{Normal Basis Representations}
+\newcommand{\ugxProblemFiniteNormalNumber}{8.11.5.}
+%
+% =====================================================================
+\begin{page}{ugxProblemFiniteNormalPage}{8.11.5. Normal Basis Representations}
+% =====================================================================
+\beginscroll
+%-% \HDindex{finite field}{ugxProblemFiniteNormalPage}{8.11.5.}{Normal Basis Representations}
+%-% \HDindex{field!finite!extension of}{ugxProblemFiniteNormalPage}{8.11.5.}{Normal Basis Representations}
+%-% \HDindex{basis!normal}{ugxProblemFiniteNormalPage}{8.11.5.}{Normal Basis Representations}
+%-% \HDindex{normal basis}{ugxProblemFiniteNormalPage}{8.11.5.}{Normal Basis Representations}
+
+Let \smath{K} be a finite extension of degree \smath{n} of the
+finite field \smath{F} and let \smath{F} have \smath{q}
+elements.
+An element \smath{x} of \smath{K} is said to be
+{\it normal} over \smath{F} if the elements
+\centerline{{\axiom{1, x**q, x**(q**2), ..., x**(q**(n-1))}}}
+form a basis of \smath{K} as a vector space over \smath{F}.
+Such a basis is called a {\it normal basis}.\footnote{This
+agrees with the general definition of a normal basis because the
+\smath{n} distinct powers of the automorphism
+\texht{$x \mapsto x^q$}{\axiom{x +-> x**q}}
+constitute the Galois group of \smath{K/F}.}
+
+If \smath{x} is normal over \smath{F}, its minimal
+%-% \HDindex{polynomial!minimal}{ugxProblemFiniteNormalPage}{8.11.5.}{Normal Basis Representations}
+polynomial is also said to be {\it normal} over \smath{F}.
+%-% \HDindex{minimal polynomial}{ugxProblemFiniteNormalPage}{8.11.5.}{Normal Basis Representations}
+There exist normal bases for all finite extensions of arbitrary
+finite fields.
+
+In \axiomType{FiniteFieldNormalBasis} (abbreviation
+\axiomType{FFNB}), the elements of the finite field are represented
+by coordinate vectors with respect to a normal basis.
+
+\xtc{
+You provide a prime \smath{p} and an extension degree
+\smath{n}.
+}{
+\spadpaste{K := FFNB(3,8) \bound{K}}
+}
+%
+\Language{} uses the prime field \axiomType{PrimeField(p)},
+here \axiomType{PrimeField 3},
+and it chooses a normal polynomial of degree
+\smath{n}, here 8, over the ground field.
+The remainder class of the indeterminate is used
+as the normal element.
+The polynomial indeterminate is automatically chosen by \Language{} and
+is typically something like \axiom{\%A} or \axiom{\%D}.
+These (strange) variables are only for output display;
+there are several ways to construct elements of this field.
+The output of the basis elements is something like
+\texht{$\%A^{q^i}.$}{
+\begin{verbatim}
+   i
+  q
+%A  .
+\end{verbatim}
+}
+%
+\xtc{
+}{
+\spadpaste{a := normalElement()\$K \bound{a}\free{K}}
+}
+%
+%
+\xtc{
+You can calculate in \smath{K} using \smath{a}.
+}{
+\spadpaste{b  := a**12 - a**5 + a \bound{b}\free{a}}
+}
+
+\axiomType{FiniteFieldNormalBasisExtension} (abbreviation
+\axiomType{FFNBX}) is
+similar to \axiomType{FiniteFieldNormalBasis} except that the
+groundfield for \axiomType{FiniteFieldNormalBasisExtension} is arbitrary and
+chosen by you.
+In case you select the prime field as ground field, there is
+essentially no difference between the constructed two finite field
+extensions.
+\xtc{
+}{
+\spadpaste{GF9 := FFNB(3,2); \bound{GF9}}
+}
+\xtc{
+}{
+\spadpaste{GF729 := FFNBX(GF9,3); \bound{GF729}\free{GF9}}
+}
+\xtc{
+}{
+\spadpaste{r := random()\$GF729 \bound{r}\free{GF729}}
+}
+\xtc{
+}{
+\spadpaste{r + r**3 + r**9 + r**27 \free{r}}
+}
+
+\axiomType{FiniteFieldNormalBasisExtensionByPolynomial}
+(abbreviation \axiomType{FFNBP}) is similar to
+\axiomType{FiniteFieldNormalBasis} and
+\axiomType{FiniteFieldNormalBasisExtension} but is more general.
+For \axiomType{FiniteFieldNormalBasisExtensionByPolynomial} you
+choose both the ground field and the irreducible polynomial used
+in the representation.
+The degree of the extension is the degree of the polynomial.
+
+%
+\xtc{
+}{
+\spadpaste{GF3 := PrimeField 3; \bound{GF3}}
+}
+\xtc{
+We use a utility operation to generate an irreducible normal
+polynomial (see \downlink{``\ugxProblemFiniteUtilityTitle''}{ugxProblemFiniteUtilityPage} in Section \ugxProblemFiniteUtilityNumber\ignore{ugxProblemFiniteUtility}).
+The polynomial has one variable that is ``anonymous'': it displays
+as a question mark.
+}{
+\spadpaste{f := createNormalPoly(4)\$FFPOLY(GF3) \free{GF3}\bound{f}}
+}
+\xtc{
+}{
+\spadpaste{GF81 := FFNBP(GF3,f); \bound{GF81}\free{f GF3}}
+}
+\xtc{
+Let's look at a random element from this field.
+}{
+\spadpaste{r := random()\$GF81 \free{GF81}\bound{r1}}
+}
+\xtc{
+}{
+\spadpaste{r * r**3 * r**9 * r**27 \free{r1}}
+}
+\xtc{
+}{
+\spadpaste{norm r \free{r1}}
+}
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugxProblemFiniteConversionTitle}{Conversion Operations for Finite Fields}
+\newcommand{\ugxProblemFiniteConversionNumber}{8.11.6.}
+%
+% =====================================================================
+\begin{page}{ugxProblemFiniteConversionPage}{8.11.6. Conversion Operations for Finite Fields}
+% =====================================================================
+\beginscroll
+%-% \HDindex{field!finite!conversions}{ugxProblemFiniteConversionPage}{8.11.6.}{Conversion Operations for Finite Fields}
+
+\labelSpace{5pc}
+%
+\xtc{
+Let \texht{$K$}{\axiom{K}} be a finite field.
+}{
+\spadpaste{K := PrimeField 3 \bound{K}}
+}
+%
+An extension field \texht{$K_m$}{\axiom{Km}} of degree
+\smath{m} over \texht{$K$}{\axiom{K}} is a subfield of an
+extension field \texht{$K_n$}{\axiom{Kn}} of degree \smath{n}
+over \texht{$K$}{\axiom{K}} if and only if \smath{m} divides
+\smath{n}.
+\texht{
+\centerline{{\begin{tabular}{ccc}}}
+\centerline{{$K_n$ }}
+\centerline{{$|$ }}
+\centerline{{$K_m$ & $\Longleftrightarrow$ & $m | n$ }}
+\centerline{{$|$ }}
+\centerline{{K}}
+\centerline{{\end{tabular}}}
+}{
+\begin{verbatim}
+Kn
+|
+Km   <==>  m | n
+|
+K
+\end{verbatim}
+}
+\axiomType{FiniteFieldHomomorphisms} provides conversion operations
+between different extensions of one
+fixed finite ground field and between different representations of
+these finite fields.
+\xtc{
+Let's choose \smath{m} and \smath{n},
+}{
+\spadpaste{(m,n) := (4,8) \bound{m n}}
+}
+\xtc{
+build the field extensions,
+}{
+\spadpaste{Km := FiniteFieldExtension(K,m) \bound{Km}\free{K m}}
+}
+\xtc{
+and pick two random elements from the smaller field.
+}{
+\spadpaste{Kn := FiniteFieldExtension(K,n) \bound{Kn}\free{K n}}
+}
+\xtc{
+}{
+\spadpaste{a1 := random()\$Km \bound{a1}\free{Km}}
+}
+\xtc{
+}{
+\spadpaste{b1 := random()\$Km \bound{b1}\free{Km}}
+}
+%
+\xtc{
+Since \smath{m} divides \smath{n},
+\texht{$K_m$}{\axiom{Km}} is a subfield of \texht{$K_n$}{\axiom{Kn}}.
+}{
+\spadpaste{a2 := a1 :: Kn \bound{a2}\free{a1 Kn}}
+}
+\xtc{
+Therefore we can convert the elements of \texht{$K_m$}{\axiom{Km}}
+into elements of \texht{$K_n$}{\axiom{Kn}}.
+}{
+\spadpaste{b2 := b1 :: Kn \bound{b2}\free{b1 Kn}}
+}
+%
+%
+\xtc{
+To check this, let's do some arithmetic.
+}{
+\spadpaste{a1+b1 - ((a2+b2) :: Km) \free{a1 a2 b1 b2 Km Kn}}
+}
+\xtc{
+}{
+\spadpaste{a1*b1 - ((a2*b2) :: Km) \free{a1 a2 b1 b2 Km Kn}}
+}
+%
+There are also conversions available for the
+situation, when \texht{$K_m$}{\axiom{Km}} and \texht{$K_n$}{\axiom{Kn}}
+are represented in different ways (see
+\downlink{``\ugxProblemFiniteExtensionFiniteTitle''}{ugxProblemFiniteExtensionFinitePage} in Section \ugxProblemFiniteExtensionFiniteNumber\ignore{ugxProblemFiniteExtensionFinite}).
+For example let's choose \texht{$K_m$}{\axiom{Km}} where the
+representation is 0 plus the cyclic multiplicative group and
+\texht{$K_n$}{\axiom{Kn}} with a normal basis representation.
+\xtc{
+}{
+\spadpaste{Km := FFCGX(K,m) \bound{Km2}\free{K m}}
+}
+\xtc{
+}{
+\spadpaste{Kn := FFNBX(K,n) \bound{Kn2}\free{K n}}
+}
+\xtc{
+}{
+\spadpaste{(a1,b1) := (random()\$Km,random()\$Km) \bound{a12 b12}\free{Km2}}
+}
+\xtc{
+}{
+\spadpaste{a2 := a1 :: Kn \bound{a22}\free{a12 Kn2}}
+}
+\xtc{
+}{
+\spadpaste{b2 := b1 :: Kn \bound{b22}\free{b12 Kn2}}
+}
+%
+\xtc{
+Check the arithmetic again.
+}{
+\spadpaste{a1+b1 - ((a2+b2) :: Km) \free{a12 a22 b12 b22 Km2}}
+}
+\xtc{
+}{
+\spadpaste{a1*b1 - ((a2*b2) :: Km) \free{a12 a22 b12 b22 Km2}}
+}
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugxProblemFiniteUtilityTitle}{Utility Operations for Finite Fields}
+\newcommand{\ugxProblemFiniteUtilityNumber}{8.11.7.}
+%
+% =====================================================================
+\begin{page}{ugxProblemFiniteUtilityPage}{8.11.7. Utility Operations for Finite Fields}
+% =====================================================================
+\beginscroll
+
+\axiomType{FiniteFieldPolynomialPackage} (abbreviation
+\axiomType{FFPOLY})
+provides operations for generating, counting and testing polynomials
+over finite fields. Let's start with a couple of definitions:
+\indent{4}
+\beginitems
+\item[-] A polynomial is {\it primitive} if its roots are primitive
+%-% \HDindex{polynomial!primitive}{ugxProblemFiniteUtilityPage}{8.11.7.}{Utility Operations for Finite Fields}
+elements in an extension of the coefficient field of degree equal
+to the degree of the polynomial.
+\item[-] A polynomial is {\it normal} over its coefficient field
+%-% \HDindex{polynomial!normal}{ugxProblemFiniteUtilityPage}{8.11.7.}{Utility Operations for Finite Fields}
+if its roots are linearly independent
+elements in an extension of the coefficient field of degree equal
+to the degree of the polynomial.
+\enditems
+\indent{0}
+In what follows, many of the generated polynomials have one
+``anonymous'' variable.
+This indeterminate is displayed as a question mark (\spadSyntax{?}).
+
+\xtc{
+To fix ideas, let's use the field with five elements for the first
+few examples.
+}{
+\spadpaste{GF5 := PF 5; \bound{GF5}}
+}
+%
+%
+\xtc{
+You can generate irreducible polynomials of any (positive) degree
+%-% \HDindex{polynomial!irreducible}{ugxProblemFiniteUtilityPage}{8.11.7.}{Utility Operations for Finite Fields}
+(within the storage capabilities of the computer and your ability
+to wait) by using
+\axiomFunFrom{createIrreduciblePoly}{FiniteFieldPolynomialPackage}.
+}{
+\spadpaste{f := createIrreduciblePoly(8)\$FFPOLY(GF5) \bound{f}\free{GF5}}
+}
+%
+\xtc{
+Does this polynomial have other important properties? Use
+\axiomFun{primitive?} to test whether it is a primitive polynomial.
+}{
+\spadpaste{primitive?(f)\$FFPOLY(GF5) \free{f}}
+}
+\xtc{
+Use \axiomFun{normal?} to test whether it is a normal polynomial.
+}{
+\spadpaste{normal?(f)\$FFPOLY(GF5) \free{f}}
+}
+\noindent
+Note that this is actually a trivial case,
+because a normal polynomial of degree \smath{n}
+must have a nonzero term of degree \smath{n-1}.
+We will refer back to this later.
+
+\xtc{
+To get a primitive polynomial of degree 8 just issue this.
+}{
+\spadpaste{p := createPrimitivePoly(8)\$FFPOLY(GF5) \bound{p}\free{GF5}}
+}
+\xtc{
+}{
+\spadpaste{primitive?(p)\$FFPOLY(GF5) \free{p}}
+}
+\xtc{
+This polynomial is not normal,
+}{
+\spadpaste{normal?(p)\$FFPOLY(GF5) \free{p}}
+}
+\xtc{
+but if you want a normal one simply write this.
+}{
+\spadpaste{n := createNormalPoly(8)\$FFPOLY(GF5) \bound{n} \free{GF5}}
+}
+\xtc{
+This polynomial is not primitive!
+}{
+\spadpaste{primitive?(n)\$FFPOLY(GF5) \free{n}}
+}
+This could have been seen directly, as
+the constant term is 1 here, which is not a primitive
+element up to the factor (\axiom{-1}) raised to the degree of the
+polynomial.\footnote{Cf. Lidl, R. \& Niederreiter,
+H., {\it Finite Fields,} Encycl. of Math. 20, (Addison-Wesley, 1983),
+p.90, Th. 3.18.}
+
+What about
+polynomials that are both primitive and normal?
+The existence of such a polynomial is by no means obvious.
+\footnote{The existence of such polynomials is proved in
+Lenstra, H. W. \& Schoof, R. J., {\it Primitive
+Normal Bases for Finite Fields,} Math. Comp. 48, 1987, pp. 217-231.}
+%
+\xtc{
+If you really need one use either
+\axiomFunFrom{createPrimitiveNormalPoly}{FiniteFieldPolynomialPackage} or
+\axiomFunFrom{createNormalPrimitivePoly}{FiniteFieldPolynomialPackage}.
+}{
+\spadpaste{createPrimitiveNormalPoly(8)\$FFPOLY(GF5) \free{GF5}}
+}
+%
+
+If you want to obtain additional polynomials of the various types above
+as given by the {\bf create...} operations above, you can use the {\bf
+next...} operations.
+For instance,
+\axiomFunFrom{nextIrreduciblePoly}{FiniteFieldPolynomialPackage} yields
+the next monic irreducible polynomial with the same degree as the input
+polynomial.
+By ``next'' we mean ``next in a natural order using the terms and
+coefficients.''
+This will become more clear in the following examples.
+
+\xtc{
+This is the field with five elements.
+}{
+\spadpaste{GF5 := PF 5; \bound{GF5}}
+}
+%
+\xtc{
+Our first example irreducible polynomial, say
+of degree 3, must be ``greater'' than this.
+}{
+\spadpaste{h := monomial(1,8)\$SUP(GF5) \bound{h}\free{GF5}}
+}
+\xtc{
+You can generate it by doing this.
+}{
+\spadpaste{nh := nextIrreduciblePoly(h)\$FFPOLY(GF5) \bound{nh}\free{h}}
+}
+%
+\xtc{
+Notice that this polynomial is not the same as the one
+\axiomFunFrom{createIrreduciblePoly}{FiniteFieldPolynomialPackage}.
+}{
+\spadpaste{createIrreduciblePoly(3)\$FFPOLY(GF5) \free{GF5}}
+}
+\xtc{
+You can step through all irreducible polynomials of degree 8 over
+the field with 5 elements by repeatedly issuing this.
+}{
+\spadpaste{nh := nextIrreduciblePoly(nh)\$FFPOLY(GF5) \free{nh}}
+}
+\xtc{
+You could also ask for the total number of these.
+}{
+\spadpaste{numberOfIrreduciblePoly(5)\$FFPOLY(GF5) \free{GF5}}
+}
+
+We hope that ``natural order'' on polynomials is now clear:
+first we compare the number of monomials of
+two polynomials (``more'' is ``greater'');
+then, if necessary, the degrees of these monomials (lexicographically),
+and lastly their coefficients (also
+lexicographically, and using the operation \axiomFun{lookup} if
+our field is not a prime field).
+Also note that we make both polynomials monic before looking at the
+coefficients:
+multiplying either polynomial  by a nonzero constant
+produces the same result.
+
+%
+\xtc{
+The package
+\axiomType{FiniteFieldPolynomialPackage} also provides similar
+operations for primitive and normal polynomials. With
+the exception of the number of primitive normal polynomials;
+we're not aware of any known formula for this.
+}{
+\spadpaste{numberOfPrimitivePoly(3)\$FFPOLY(GF5) \free{GF5}}
+}
+%
+%
+\xtc{
+Take these,
+}{
+\spadpaste{m := monomial(1,1)\$SUP(GF5) \bound{m}\free{GF5}}
+}
+\xtc{
+}{
+\spadpaste{f := m**3 + 4*m**2 + m + 2 \bound{fx}\free{m}}
+}
+%
+%
+\xtc{
+and then we have:
+}{
+\spadpaste{f1 := nextPrimitivePoly(f)\$FFPOLY(GF5) \free{fx}\bound{f1}}
+}
+\xtc{
+What happened?
+}{
+\spadpaste{nextPrimitivePoly(f1)\$FFPOLY(GF5) \free{f1}}
+}
+%
+Well, for the ordering used in
+\axiomFunFrom{nextPrimitivePoly}{FiniteFieldPolynomialPackage} we
+use as first criterion a comparison of the constant terms of the
+polynomials.
+Analogously, in
+\axiomFunFrom{nextNormalPoly}{FiniteFieldPolynomialPackage} we first
+compare the monomials of degree 1 less than the degree of the
+polynomials (which is nonzero, by an earlier remark).
+%
+\xtc{
+}{
+\spadpaste{f := m**3 + m**2 + 4*m + 1 \bound{fy} \free{m}}
+}
+\xtc{
+}{
+\spadpaste{f1 := nextNormalPoly(f)\$FFPOLY(GF5) \free{fy}\bound{f1y}}
+}
+\xtc{
+}{
+\spadpaste{nextNormalPoly(f1)\$FFPOLY(GF5) \free{f1y}}
+}
+%
+\noindent
+We don't have to restrict ourselves to prime fields.
+%
+\xtc{
+Let's consider, say, a field with 16 elements.
+}{
+\spadpaste{GF16 := FFX(FFX(PF 2,2),2); \bound{GF16} }
+}
+%
+%
+\xtc{
+We can apply any of the operations described above.
+}{
+\spadpaste{createIrreduciblePoly(5)\$FFPOLY(GF16) \free{GF16}}
+}
+
+\xtc{
+\Language{} also provides operations
+for producing random polynomials of a given degree
+}{
+\spadpaste{random(5)\$FFPOLY(GF16) \free{GF16}}
+}
+\xtc{
+or with degree between two given bounds.
+}{
+\spadpaste{random(3,9)\$FFPOLY(GF16) \free{GF16}}
+}
+
+\axiomType{\axiom{FiniteFieldPolynomialPackage2}} (abbreviation
+\axiomType{FFPOLY2})
+exports an operation \axiomFun{rootOfIrreduciblePoly}
+for finding one root of an irreducible polynomial \axiom{f}
+%-% \HDindex{polynomial!root of}{ugxProblemFiniteUtilityPage}{8.11.7.}{Utility Operations for Finite Fields}
+in an extension field of the coefficient field.
+The degree of the extension has to be a multiple of the degree of \axiom{f}.
+It is not checked whether \axiom{f} actually is irreducible.
+
+%
+\xtc{
+To illustrate this operation, we fix a ground field \axiom{GF}
+}{
+\spadpaste{GF2 := PrimeField 2; \bound{GF2}}
+}
+%
+%
+\xtc{
+and then an extension field.
+}{
+\spadpaste{F := FFX(GF2,12) \bound{F}\free{GF2}}
+}
+%
+%
+\xtc{
+We construct an irreducible polynomial over \axiom{GF2}.
+}{
+\spadpaste{f := createIrreduciblePoly(6)\$FFPOLY(GF2) \bound{fz}\free{GF2}}
+}
+%
+%
+\xtc{
+We compute a root of \axiom{f}.
+}{
+\spadpaste{root := rootOfIrreduciblePoly(f)\$FFPOLY2(F,GF2) \free{F GF2 fz}\bound{root}}
+}
+%
+%and check the result
+%\spadcommand{eval(f, monomial(1,1)\$SUP(F) = root) \free{fz F root}}
+
+%*********************************************************************
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugProblemIdealTitle}{Primary Decomposition of Ideals}
+\newcommand{\ugProblemIdealNumber}{8.12.}
+%
+% =====================================================================
+\begin{page}{ugProblemIdealPage}{8.12. Primary Decomposition of Ideals}
+% =====================================================================
+\beginscroll
+%*********************************************************************
+%
+\Language{} provides a facility for the primary decomposition
+%-% \HDindex{ideal!primary decomposition}{ugProblemIdealPage}{8.12.}{Primary Decomposition of Ideals}
+of
+%-% \HDindex{primary decomposition of ideal}{ugProblemIdealPage}{8.12.}{Primary Decomposition of Ideals}
+polynomial ideals over fields of characteristic zero.
+The algorithm
+%is discussed in \cite{gtz:gbpdpi} and
+works in essentially two steps:
+\indent{4}
+\beginitems
+\item[1. ] the problem is solved for 0-dimensional ideals by ``generic''
+projection on the last coordinate
+\item[2. ] a ``reduction process'' uses localization and ideal quotients
+to reduce the general case to the 0-dimensional one.
+\enditems
+\indent{0}
+The \Language{} constructor \axiomType{PolynomialIdeals}
+represents ideals with coefficients in any field and
+supports the basic ideal operations,
+including intersection, sum and quotient.
+\axiomType{IdealDecompositionPackage} contains the specific
+operations for the primary decomposition and the computation of the
+radical of an ideal with polynomial
+coefficients in a field of characteristic 0 with
+an effective algorithm for factoring polynomials.
+
+The following examples illustrate the capabilities of this facility.
+%
+\xtc{
+First consider the ideal generated by
+\texht{$x^2 + y^2 - 1$}{\axiom{x**2 + y**2 - 1}}
+(which defines a circle in the \axiom{(x,y)}-plane) and the ideal
+generated by \texht{$x^2 - y^2$}{\axiom{x**2 - y**2}} (corresponding to the
+straight lines \axiom{x = y} and \axiom{x = -y}.
+}{
+\spadpaste{(n,m) : List DMP([x,y],FRAC INT) \bound{nm}}
+}
+\xtc{
+}{
+\spadpaste{m := [x**2+y**2-1] \free{nm} \bound{m}}
+}
+\xtc{
+}{
+\spadpaste{n := [x**2-y**2] \free{nm} \bound{n}}
+}
+%
+%
+\xtc{
+We find the equations defining the intersection of the two loci.
+This correspond to the sum of the associated ideals.
+}{
+\spadpaste{id := ideal m  + ideal n \free{n m} \bound{id}}
+}
+%
+%
+\xtc{
+We can check if the locus contains only a finite number of points,
+that is, if the ideal is zero-dimensional.
+}{
+\spadpaste{zeroDim? id \free{id}}
+}
+\xtc{
+}{
+\spadpaste{zeroDim?(ideal m) \free{m}}
+}
+\xtc{
+}{
+\spadpaste{dimension ideal m \free{m}}
+}
+\xtc{
+We can find polynomial relations among the generators
+(\axiom{f} and \axiom{g} are the parametric equations of the knot).
+}{
+\spadpaste{(f,g):DMP([x,y],FRAC INT) \bound{fg}}
+}
+\xtc{
+}{
+\spadpaste{f := x**2-1 \free{fg} \bound{f}}
+}
+\xtc{
+}{
+\spadpaste{g := x*(x**2-1) \free{fg} \bound{g}}
+}
+\xtc{
+}{
+\spadpaste{relationsIdeal [f,g] \free{f g}}
+}
+
+\xtc{
+We can compute the primary decomposition of an ideal.
+}{
+\spadpaste{l: List DMP([x,y,z],FRAC INT) \bound{ll}}
+}
+\xtc{
+}{
+\spadpaste{l:=[x**2+2*y**2,x*z**2-y*z,z**2-4] \free{ll} \bound{l}}
+}
+\xtc{
+}{
+\spadpaste{ld:=primaryDecomp ideal l \free{l} \bound{ld}}
+}
+\xtc{
+We can intersect back.
+}{
+\spadpaste{reduce(intersect,ld) \free{ld}}
+}
+
+\xtc{
+We can compute the radical of every primary component.
+}{
+\spadpaste{reduce(intersect,[radical ld.i for i in 1..2]) \free{ld}}
+}
+\xtc{
+Their intersection is equal to the radical of the ideal of \axiom{l}.
+}{
+\spadpaste{radical ideal l \free{l}}
+}
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugProblemGaloisTitle}{Computation of Galois Groups}
+\newcommand{\ugProblemGaloisNumber}{8.13.}
+%
+% =====================================================================
+\begin{page}{ugProblemGaloisPage}{8.13. Computation of Galois Groups}
+% =====================================================================
+\beginscroll
+%
+As a sample use of \Language{}'s algebraic number facilities,
+%-% \HDindex{group!Galois}{ugProblemGaloisPage}{8.13.}{Computation of Galois Groups}
+we compute
+%-% \HDindex{Galois!group}{ugProblemGaloisPage}{8.13.}{Computation of Galois Groups}
+the Galois group of the polynomial
+\texht{$p(x) = x^5 - 5 x + 12$}{\axiom{p(x) = x**5 - 5*x + 12}}.
+%
+\xtc{
+}{
+\spadpaste{p := x**5 - 5*x + 12 \bound{p}}
+}
+%
+We would like to construct a polynomial \smath{f(x)}
+such that the splitting
+%-% \HDindex{field!splitting}{ugProblemGaloisPage}{8.13.}{Computation of Galois Groups}
+field
+%-% \HDindex{splitting field}{ugProblemGaloisPage}{8.13.}{Computation of Galois Groups}
+of \smath{p(x)} is generated by one root of \smath{f(x)}.
+First we construct a polynomial \smath{r = r(x)} such that one
+root of \smath{r(x)} generates the field generated by two roots of
+the polynomial \smath{p(x)}.
+(As it will turn out, the field generated by two roots of
+\smath{p(x)} is, in fact, the splitting field of
+\smath{p(x)}.)
+
+From the proof of the primitive element theorem we know that
+if \smath{a} and \smath{b} are
+algebraic numbers, then the field
+\texht{${\bf Q}(a,b)$}{\axiom{Q(a,b)}} is equal to
+\texht{${\bf Q}(a+kb)$}{\axiom{Q(a + k*b)}} for an
+appropriately chosen integer \smath{k}.
+In our case, we construct the minimal polynomial of
+\texht{$a_i - a_j$}{\axiom{a[i] - a[j]}}, where
+\texht{$a_i$}{\axiom{a[i]}} and
+\texht{$a_j$}{\axiom{a[j]}} are two roots of \smath{p(x)}.
+We construct this polynomial using \axiomFun{resultant}.
+The main result we  need is the following:
+If \smath{f(x)} is a polynomial with roots
+\texht{$a_i \ldots a_m$}{\axiom{a[1]...a[m]}} and
+\smath{g(x)} is a polynomial
+with roots
+\texht{$b_i \ldots b_n$}{\axiom{b[1]...b[n]}}, then the polynomial
+\axiom{h(x) = resultant(f(y), g(x-y), y)}
+is a polynomial of degree \smath{m*n} with
+roots
+\texht{$a_i + b_j, i = 1 \ldots m, j = 1 \ldots n$}
+{\axiom{a[i] + b[j], 1 <= i <= m, 1 <= j <= n}}.
+
+\xtc{
+For \smath{f(x)} we use the polynomial \smath{p(x)}.
+For \smath{g(x)} we use the polynomial \smath{-p(-x)}.
+Thus, the polynomial we  first construct is
+\axiom{resultant(p(y), -p(y-x), y)}.
+}{
+\spadpaste{q := resultant(eval(p,x,y),-eval(p,x,y-x),y) \free{p} \bound{q}}
+}
+%
+The roots of \smath{q(x)} are
+\texht{$a_i - a_j, i \leq 1, j \leq 5$}
+{\axiom{a[i] - a[j], 1 <= i,j <= 5}}.
+Of course, there are five pairs \smath{(i,j)} with \smath{i = j},
+so \axiom{0} is a 5-fold root of \smath{q(x)}.
+%
+\xtc{
+Let's get rid of this factor.
+}{
+\spadpaste{q1 := exquo(q, x**5) \free{q} \bound{q1}}
+}
+\xtc{
+Factor the polynomial \axiom{q1}.
+}{
+\spadpaste{factoredQ := factor q1 \free{q1} \bound{factoredQ}}
+}
+%
+We see that \axiom{q1} has two irreducible factors, each of degree \axiom{10}.
+(The fact that the polynomial \axiom{q1} has two factors of
+degree \axiom{10} is enough to show
+that the Galois group of \smath{p(x)} is the dihedral group of
+order \axiom{10}.\footnote{See McKay, Soicher,  Computing Galois Groups
+over the Rationals, Journal of Number Theory 20, 273-281 (1983).
+We do not assume the results of this paper, however, and we continue with
+the computation.}
+Note that the type of \axiom{factoredQ} is \axiomType{FR POLY INT}, that is,
+\axiomType{Factored Polynomial Integer}.
+%-% \HDexptypeindex{Factored}{ugProblemGaloisPage}{8.13.}{Computation of Galois Groups}
+This is a special data type for recording factorizations of polynomials with
+integer coefficients (see \downlink{`Factored'}{FactoredXmpPage}\ignore{Factored}).
+%
+\xtc{
+We can access the individual factors using the operation
+\axiomFunFrom{nthFactor}{Factored}.
+}{
+\spadpaste{r := nthFactor(factoredQ,1) \free{factoredQ} \bound{r}}
+}
+%
+
+Consider the polynomial \smath{r = r(x)}.
+This is the minimal polynomial of the difference of two roots of
+\smath{p(x)}.
+Thus, the splitting field of \smath{p(x)} contains a subfield of
+degree \axiom{10}.
+We show that this subfield is, in fact, the splitting field of
+\smath{p(x)} by showing that \smath{p(x)} factors completely
+over this field.
+%
+\xtc{
+First we create a symbolic root of the polynomial \smath{r(x)}.
+(We replaced \axiom{x} by \axiom{b} in the
+polynomial \axiom{r} so that our symbolic root would be
+printed as \axiom{b}.)
+}{
+\spadpaste{beta:AN := rootOf(eval(r,x,b)) \free{r} \bound{beta}}
+}
+\xtc{
+We next tell \Language{} to view \smath{p(x)} as a univariate polynomial
+in \axiom{x}
+with algebraic number coefficients.
+This is accomplished with this type declaration.
+}{
+\spadpaste{p := p::UP(x,INT)::UP(x,AN) \free{p} \bound{declareP}}
+}
+%
+%
+\xtc{
+Factor \smath{p(x)} over the field
+\texht{${\bf Q}(\beta)$}{\axiom{Q(beta)}}.
+(This computation will take some time!)
+}{
+\spadpaste{algFactors := factor(p,[beta]) \free{declareP beta} \bound{algFactors}}
+}
+%
+When factoring over number fields, it is important to specify the
+field over which the polynomial is to be factored, as polynomials
+have different factorizations over different fields.
+When you use the operation \axiomFun{factor}, the field over which
+the polynomial is factored is the field generated by
+\indent{4}
+\beginitems
+\item[1. ] the algebraic numbers that appear
+in the coefficients of the polynomial, and
+\item[2. ] the algebraic numbers that
+appear in a list passed as an optional second argument of the operation.
+\enditems
+\indent{0}
+In our case, the coefficients of \axiom{p}
+are all rational integers and only \axiom{beta}
+appears in the list, so the field is simply
+\texht{${\bf Q}(\beta)$}{\axiom{Q(beta)}}.
+%
+\xtc{
+It was necessary to give the list \axiom{[beta]}
+as a second argument of the operation
+because otherwise the polynomial would have been factored over the field
+generated by its coefficients, namely the rational numbers.
+}{
+\spadpaste{factor(p) \free{declareP}}
+}
+%
+We have shown that the splitting field of \smath{p(x)} has degree
+\axiom{10}.
+Since the symmetric group of degree 5 has only one transitive subgroup
+of order \axiom{10}, we know that the Galois group of \smath{p(x)} must be
+this group, the dihedral group
+%-% \HDindex{group!dihedral}{ugProblemGaloisPage}{8.13.}{Computation of Galois Groups}
+of order \axiom{10}.
+Rather than stop here, we explicitly compute the action of the Galois
+group on the roots of \smath{p(x)}.
+
+First we assign the roots of \smath{p(x)} as the values of five
+%-% \HDindex{root}{ugProblemGaloisPage}{8.13.}{Computation of Galois Groups}
+variables.
+\xtc{
+We can obtain an individual root by negating the constant coefficient of
+one of the factors of \smath{p(x)}.
+}{
+\spadpaste{factor1 := nthFactor(algFactors,1) \free{algFactors} \bound{factor1}}
+}
+\xtc{
+}{
+\spadpaste{root1 := -coefficient(factor1,0) \free{factor1} \bound{root1}}
+}
+%
+%
+\xtc{
+We can obtain a list of all the roots in this way.
+}{
+\spadpaste{roots := [-coefficient(nthFactor(algFactors,i),0) for i in 1..5] \free{algFactors} \bound{roots}}
+}
+
+The expression
+\begin{verbatim}
+- coefficient(nthFactor(algFactors, i), 0)}
+\end{verbatim}
+is the \eth{\axiom{i }} root
+of \smath{p(x)} and the elements of \axiom{roots} are the \eth{\axiom{i }}
+roots of \smath{p(x)} as \axiom{i} ranges from \axiom{1} to \axiom{5}.
+
+\xtc{
+Assign the roots as the values of the variables \axiom{a1,...,a5}.
+}{
+\spadpaste{(a1,a2,a3,a4,a5) := (roots.1,roots.2,roots.3,roots.4,roots.5) \free{roots} \bound{ais}}
+}
+%
+
+Next we express the roots of \smath{r(x)} as polynomials in
+\axiom{beta}.
+We could obtain these roots by calling the operation \axiomFun{factor}:
+\axiom{factor(r, [beta])} factors \axiom{r(x)} over
+\texht{${\bf Q}(\beta)$}{\axiom{Q(beta)}}.
+However, this is a lengthy computation and we can obtain the roots of
+\smath{r(x)} as differences of the roots \axiom{a1,...,a5} of
+\smath{p(x)}.
+Only ten of these differences are roots of \smath{r(x)} and the
+other ten are roots
+of the other irreducible factor of \axiom{q1}.
+We can determine if a given value is a root of \smath{r(x)} by evaluating
+\smath{r(x)} at that particular value.
+(Of course, the order in which factors are returned by the
+operation \axiomFun{factor}
+is unimportant and may change with different implementations of the operation.
+Therefore, we cannot predict in advance which differences are roots of
+\smath{r(x)} and which are not.)
+%
+\xtc{
+Let's look at four examples (two are roots of \smath{r(x)} and
+two are not).
+}{
+\spadpaste{eval(r,x,a1 - a2) \free{ais}}
+}
+\xtc{
+}{
+\spadpaste{eval(r,x,a1 - a3) \free{ais}}
+}
+\xtc{
+}{
+\spadpaste{eval(r,x,a1 - a4) \free{ais}}
+}
+\xtc{
+}{
+\spadpaste{eval(r,x,a1 - a5) \free{ais}}
+}
+%
+
+Take one of the differences that was a root of \smath{r(x)}
+and assign it to the variable \axiom{bb}.
+\xtc{
+For example, if \axiom{eval(r,x,a1 - a4)} returned \axiom{0}, you would
+enter this.
+}{
+\spadpaste{bb := a1 - a4 \free{ais} \bound{bb}}
+}
+Of course, if the difference is, in fact, equal to the root \axiom{beta},
+you should choose another root of \smath{r(x)}.
+
+Automorphisms of the splitting field are given by mapping a generator of
+the field, namely \axiom{beta}, to other roots of its minimal polynomial.
+Let's see what happens when \axiom{beta} is mapped to \axiom{bb}.
+%
+\xtc{
+We compute the images of the roots \axiom{a1,...,a5}
+under this automorphism:
+}{
+\spadpaste{aa1 := subst(a1,beta = bb) \free{beta bb ais} \bound{aa1}}
+}
+\xtc{
+}{
+\spadpaste{aa2 := subst(a2,beta = bb) \free{beta bb ais} \bound{aa2}}
+}
+\xtc{
+}{
+\spadpaste{aa3 := subst(a3,beta = bb) \free{beta bb ais} \bound{aa3}}
+}
+\xtc{
+}{
+\spadpaste{aa4 := subst(a4,beta = bb) \free{beta bb ais} \bound{aa4}}
+}
+\xtc{
+}{
+\spadpaste{aa5 := subst(a5,beta = bb) \free{beta bb ais} \bound{aa5}}
+}
+%
+Of course, the values \axiom{aa1,...,aa5} are simply a permutation of the values
+\axiom{a1,...,a5}.
+%
+\xtc{
+Let's find the value of \axiom{aa1} (execute as many of the following five commands
+as necessary).
+}{
+\spadpaste{(aa1 = a1) :: Boolean \free{aa1}}
+}
+\xtc{
+}{
+\spadpaste{(aa1 = a2) :: Boolean \free{aa1}}
+}
+\xtc{
+}{
+\spadpaste{(aa1 = a3) :: Boolean \free{aa1}}
+}
+\xtc{
+}{
+\spadpaste{(aa1 = a4) :: Boolean \free{aa1}}
+}
+\xtc{
+}{
+\spadpaste{(aa1 = a5) :: Boolean \free{aa1}}
+}
+%
+Proceeding in this fashion, you can find the values of
+\axiom{aa2,...aa5}.\footnote{Here you should use the
+\Clef{} line editor.
+See \downlink{``\ugAvailCLEFTitle''}{ugAvailCLEFPage} in Section \ugAvailCLEFNumber\ignore{ugAvailCLEF} for more information about \Clef{}.}
+You have represented the automorphism \axiom{beta -> bb}
+as a permutation of the roots \axiom{a1,...,a5}.
+If you wish, you can repeat this computation for all the roots of
+\smath{r(x)} and represent the Galois group of
+\smath{p(x)} as a subgroup of the symmetric group on five letters.
+
+Here are two other problems that you may attack in a similar fashion:
+\indent{4}
+\beginitems
+\item[1. ] Show that the Galois group of
+\texht{$p(x) = x^4 + 2 x^3 - 2 x^2 - 3 x + 1$}{\axiom{
+p(x) = x**4 + 2*x**3 - 2*x**2 - 3*x + 1}}
+is the dihedral group of order eight.
+%-% \HDindex{group!dihedral}{ugProblemGaloisPage}{8.13.}{Computation of Galois Groups}
+(The splitting field of this polynomial is the Hilbert class field
+%-% \HDindex{Hilbert class field}{ugProblemGaloisPage}{8.13.}{Computation of Galois Groups}
+of
+%-% \HDindex{field!Hilbert class}{ugProblemGaloisPage}{8.13.}{Computation of Galois Groups}
+the quadratic field
+\texht{${\bf Q}(\sqrt{145})$}{\axiom{Q(sqrt(145))}}.)
+\item[2. ] Show that the Galois group of
+\texht{$p(x) = x^6 + 108$}{\axiom{p(x) = x**6 + 108}}
+has order 6 and is
+isomorphic to \texht{$S_3,$}{} the symmetric group on three letters.
+%-% \HDindex{group!symmetric}{ugProblemGaloisPage}{8.13.}{Computation of Galois Groups}
+(The splitting field of this polynomial is the splitting field of
+\texht{$x^3 - 2$}{\axiom{x**3 - 2}}.)
+\enditems
+\indent{0}
+
+\endscroll
+\autobuttons
+\end{page}
+%
+%
+\newcommand{\ugProblemGeneticTitle}{Non-Associative Algebras and Modelling Genetic Laws}
+\newcommand{\ugProblemGeneticNumber}{8.14.}
+%
+% =====================================================================
+\begin{page}{ugProblemGeneticPage}{8.14. Non-Associative Algebras and Modelling Genetic Laws}
+% =====================================================================
+\beginscroll
+
+Many algebraic structures of mathematics and \Language{}
+have a multiplication operation \axiomOp{*} that satisfies
+the associativity law
+%-% \HDindex{associativity law}{ugProblemGeneticPage}{8.14.}{Non-Associative Algebras and Modelling Genetic Laws}
+\texht{$a*(b*c) = (a*b)*c$}{\spad{a*(b*c) = (a*b)*c}}
+for all \smath{a}, \smath{b} and \smath{c}.
+The octonions (see \downlink{`Octonion'}{OctonionXmpPage}\ignore{Octonion}) are a well known exception.
+There are many other interesting non-associative structures, such as the
+class of
+%-% \HDindex{Lie algebra}{ugProblemGeneticPage}{8.14.}{Non-Associative Algebras and Modelling Genetic Laws}
+Lie algebras.\footnote{Two \Language{} implementations of Lie algebras
+are \spadtype{LieSquareMatrix} and \spadtype{FreeNilpotentLie}.}
+Lie algebras can be used, for example, to analyse Lie symmetry algebras of
+%-% \HDindex{symmetry}{ugProblemGeneticPage}{8.14.}{Non-Associative Algebras and Modelling Genetic Laws}
+partial differential
+%-% \HDindex{differential equation!partial}{ugProblemGeneticPage}{8.14.}{Non-Associative Algebras and Modelling Genetic Laws}
+equations.
+%-% \HDindex{partial differential equation}{ugProblemGeneticPage}{8.14.}{Non-Associative Algebras and Modelling Genetic Laws}
+In this section we show a different application of non-associative algebras,
+%-% \HDindex{non-associative algebra}{ugProblemGeneticPage}{8.14.}{Non-Associative Algebras and Modelling Genetic Laws}
+the modelling of genetic laws.
+%-% \HDindex{algebra!non-associative}{ugProblemGeneticPage}{8.14.}{Non-Associative Algebras and Modelling Genetic Laws}
+
+The \Language{} library contains several constructors for
+creating non-assoc\-i\-a\-tive structures,
+ranging from the categories \spadtype{Monad},
+\spadtype{NonAssociativeRng}, and
+\spadtype{FramedNonAssociativeAlgebra}, to the domains
+\spadtype{AlgebraGivenByStructuralConstants} and
+\spadtype{GenericNonAssociativeAlgebra}.
+Furthermore, the package \spadtype{AlgebraPackage} provides
+operations for analysing the structure of such algebras.\footnote{%
+The interested reader can learn more about these aspects of the
+\Language{} library from the paper
+``Computations in Algebras of Finite Rank,''
+by Johannes Grabmeier and Robert Wisbauer,
+Technical Report, IBM Heidelberg Scientific Center, 1992.}
+
+Mendel's genetic laws are often written in a form like
+\texht{
+\narrowDisplay{Aa \times Aa = {1\over 4}AA + {1\over 2}Aa + {1\over 4}aa.}}{
+\spad{Aa * Aa = (1/4)*AA + (1/2)*Aa + (1/4)*aa.}
+}
+The implementation of general algebras in \Language{} allows us to
+%-% \HDindex{Mendel's genetic laws}{ugProblemGeneticPage}{8.14.}{Non-Associative Algebras and Modelling Genetic Laws}
+use this as the definition for multiplication in an algebra.
+%-% \HDindex{genetics}{ugProblemGeneticPage}{8.14.}{Non-Associative Algebras and Modelling Genetic Laws}
+Hence, it is possible to study
+questions of genetic inheritance using \Language{}.
+To demonstrate this more precisely,  we discuss one example from a
+monograph of \texht{A. W\"orz-Busekros}{A. Woerz-Busekros},
+where you can also find a general setting of this theory.\footnote{%
+\texht{W\"{o}rz-Busekros}{Woerz-Busekros}, A.,
+{\it Algebras in Genetics},
+Springer Lectures Notes in Biomathematics 36, Berlin e.a. (1980).
+In particular, see example 1.3.}
+
+We assume that there is an infinitely large random mating population.
+Random mating of two gametes \texht{$a_i$}{\spad{ai}}  and
+\texht{$a_j$}{\spad{aj}} gives zygotes
+%-% \HDindex{zygote}{ugProblemGeneticPage}{8.14.}{Non-Associative Algebras and Modelling Genetic Laws}
+\texht{$a_ia_j$}{\spad{ai aj}}, which produce new gametes.
+%-% \HDindex{gamete}{ugProblemGeneticPage}{8.14.}{Non-Associative Algebras and Modelling Genetic Laws}
+In classical Mendelian segregation we have
+\texht{$a_ia_j = {1 \over 2}a_i+{1 \over 2}a_j$}{\spad{ai aj = (1/2)*ai+(1/2)*aj}}.
+In general, we have
+\texht{\narrowDisplay{a_ia_j = \sum_{k=1}^n \gamma_{i,j}^k\ a_k.}}%
+{\spad{ai aj = gammaij1 a1 + gammaij2 a2 + ... + gammaijn an}}
+The segregation rates \texht{$\gamma_{i,j}$}{\spad{gammaij}}
+are  the structural constants of an
+\smath{n}-dimensional algebra.
+This is provided in \Language{} by
+the constructor \spadtype{AlgebraGivenByStructuralConstants}
+(abbreviation \spadtype{ALGSC}).
+
+Consider two coupled autosomal loci with alleles
+\smath{A},\smath{a}, \smath{B}, and \smath{b}, building four
+different gametes
+\texht{$a_1 =  AB, a_2 =  Ab, a_3 =  aB,$ and $a_4 =  ab$}%
+{\spad{a1 :=  AB, a2 :=  Ab, a3 :=  aB,} and \spad{a4 :=  ab}}.
+The zygotes \texht{$a_ia_j$}{\spad{ai aj}}
+produce gametes \texht{$a_i$}{\spad{ai}} and
+\texht{$a_j$}{\spad{aj}}
+with classical Mendelian segregation.
+Zygote \texht{$a_1a_4$}{a1 a4} undergoes transition
+to \texht{$a_2a_3$}{\spad{a2 a3}}
+and vice versa with probability
+\texht{$0 \le \theta \le {1\over 2}$}{0 <= theta <= 1/2}.
+
+\xtc{
+Define a list
+\texht{$[(\gamma_{i,j}^k) 1 \le k \le 4]$}{\spad{[(gammaijk) 1 <= k <= 4]}}
+of four four-by-four matrices giving the segregation
+rates.
+We use the value \smath{1/10} for \smath{\theta}.
+}{
+\spadpaste{segregationRates : List SquareMatrix(4,FRAC INT) := [matrix [ [1, 1/2, 1/2, 9/20], [1/2, 0, 1/20, 0], [1/2, 1/20, 0, 0], [9/20, 0, 0, 0] ], matrix [ [0, 1/2, 0, 1/20], [1/2, 1, 9/20, 1/2], [0, 9/20, 0, 0], [1/20, 1/2, 0, 0] ], matrix [ [0, 0, 1/2, 1/20], [0, 0, 9/20, 0], [1/2, 9/20, 1, 1/2], [1/20, 0, 1/2, 0] ], matrix [ [0, 0, 0, 9/20], [0, 0, 1/20, 1/2], [0, 1/20, 0, 1/2], [9/20, 1/2, 1/2, 1] ] ] \bound{segregationRates}}
+}
+\xtc{
+Choose the appropriate symbols for the basis of gametes,
+}{
+\spadpaste{gametes := ['AB,'Ab,'aB,'ab]  \bound{gametes}}
+}
+\xtc{
+Define the algebra.
+}{
+\spadpaste{A := ALGSC(FRAC INT, 4, gametes, segregationRates);\bound{A}\free{gametes, segregationRates}}
+}
+
+\xtc{
+What are the probabilities for zygote
+\texht{$a_1a_4$}{a1 a4} to produce the different gametes?
+}{
+\spadpaste{a := basis()\$A;  a.1*a.4}
+}
+
+Elements in this algebra whose coefficients sum to one play a
+distinguished role.
+They represent a population with the distribution of gametes
+reflected by the coefficients with respect to the basis of
+gametes.
+
+Random mating of different populations \smath{x} and \smath{y} is described by
+their product \smath{x*y}.
+
+\xtc{
+This product is commutative only
+if the gametes are not sex-dependent, as in our example.
+}{
+\spadpaste{commutative?()\$A \free{A}}
+}
+\xtc{
+In general, it is not associative.
+}{
+\spadpaste{associative?()\$A \free{A}}
+}
+
+Random mating within a population \smath{x} is described by
+\smath{x*x.}
+The next generation is \smath{(x*x)*(x*x).}
+
+\xtc{
+Use decimal numbers to compare the distributions more easily.
+}{
+\spadpaste{x : ALGSC(DECIMAL, 4, gametes, segregationRates) :=  convert [3/10, 1/5, 1/10, 2/5]\bound{x}\free{gametes, segregationRates}}
+}
+\xtc{
+To compute directly the gametic distribution in the fifth
+generation, we use \spadfun{plenaryPower}.
+}{
+\spadpaste{plenaryPower(x,5) \free{x}}
+}
+
+We now ask two questions:
+Does this distribution converge to an equilibrium state?
+What are the distributions that are stable?
+
+\xtc{
+This is an invariant of the algebra and it is used to answer
+the first question.
+The new indeterminates describe a symbolic distribution.
+}{
+\spadpaste{q := leftRankPolynomial()\$GCNAALG(FRAC INT, 4, gametes, segregationRates) :: UP(Y, POLY FRAC INT)\bound{q}\free{gametes, segregationRates}}
+}
+\xtc{
+Because the coefficient \texht{${9 \over 20}$}{\axiom{9/20}} has absolute
+value less than 1, all distributions do converge,
+by a theorem of this theory.
+}{
+\spadpaste{factor(q :: POLY FRAC INT) \free{q}}
+}
+\xtc{
+The second question is answered by searching for idempotents in the
+algebra.
+}{
+\spadpaste{cI := conditionsForIdempotents()\$GCNAALG(FRAC INT, 4, gametes, segregationRates)\bound{cI} \free{gametes, segregationRates}}
+}
+\xtc{
+Solve these equations and look at the first solution.
+}{
+\spadpaste{gbs:= groebnerFactorize cI; gbs.1\free{cI}\bound{gbs}}
+}
+
+Further analysis using the package \spadtype{PolynomialIdeals}
+shows that there is a two-dimensional variety of equilibrium states and all
+other solutions are contained in it.
+
+\xtc{
+Choose one equilibrium state by setting two indeterminates to
+concrete values.
+}{
+\spadpaste{sol := solve concat(gbs.1,[\%x1-1/10,\%x2-1/10])\bound{sol} \free{gbs}}
+}
+\xtc{
+}{
+\spadpaste{e : A := represents reverse (map(rhs, sol.1) :: List FRAC INT)\bound{e} \free{A, sol}}
+}
+\xtc{
+Verify the result.
+}{
+\spadpaste{e*e-e \free{e}}
+}
+\endscroll
+\autobuttons
+\end{page}
+%