% 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}
%