From ab8cc85adde879fb963c94d15675783f2cf4b183 Mon Sep 17 00:00:00 2001 From: dos-reis Date: Tue, 14 Aug 2007 05:14:52 +0000 Subject: Initial population. --- src/hyper/pages/ug08.ht | 4986 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 4986 insertions(+) create mode 100644 src/hyper/pages/ug08.ht (limited to 'src/hyper/pages/ug08.ht') diff --git a/src/hyper/pages/ug08.ht b/src/hyper/pages/ug08.ht new file mode 100644 index 00000000..f14e0289 --- /dev/null +++ b/src/hyper/pages/ug08.ht @@ -0,0 +1,4986 @@ +% Copyright The Numerical Algorithms Group Limited 1992-94. All rights reserved. +% !! DO NOT MODIFY THIS FILE BY HAND !! Created by ht.awk. +\texht{\setcounter{chapter}{7}}{} % Chapter 8 + +% +\newcommand{\ugProblemTitle}{Advanced Problem Solving} +\newcommand{\ugProblemNumber}{8.} +% +% ===================================================================== +\begin{page}{ugProblemPage}{8. Advanced Problem Solving} +% ===================================================================== +\beginscroll + +In this chapter we describe techniques useful in solving advanced problems +with \Language{}. + +\beginmenu + \menudownlink{{8.1. Numeric Functions}}{ugProblemNumericPage} + \menudownlink{{8.2. Polynomial Factorization}}{ugProblemFactorPage} + \menudownlink{{8.3. Manipulating Symbolic Roots of a Polynomial}}{ugProblemSymRootPage} + \menudownlink{{8.4. Computation of Eigenvalues and Eigenvectors}}{ugProblemEigenPage} + \menudownlink{{8.5. Solution of Linear and Polynomial Equations}}{ugProblemLinPolEqnPage} + \menudownlink{{8.6. Limits}}{ugProblemLimitsPage} + \menudownlink{{8.7. Laplace Transforms}}{ugProblemLaplacePage} + \menudownlink{{8.8. Integration}}{ugProblemIntegrationPage} + \menudownlink{{8.9. Working with Power Series}}{ugProblemSeriesPage} + \menudownlink{{8.10. Solution of Differential Equations}}{ugProblemDEQPage} + \menudownlink{{8.11. Finite Fields}}{ugProblemFinitePage} + \menudownlink{{8.12. Primary Decomposition of Ideals}}{ugProblemIdealPage} + \menudownlink{{8.13. Computation of Galois Groups}}{ugProblemGaloisPage} + \menudownlink{{8.14. Non-Associative Algebras and Modelling Genetic Laws}}{ugProblemGeneticPage} +\endmenu +\endscroll +\autobuttons +\end{page} +% +% +\newcommand{\ugProblemNumericTitle}{Numeric Functions} +\newcommand{\ugProblemNumericNumber}{8.1.} +% +% ===================================================================== +\begin{page}{ugProblemNumericPage}{8.1. Numeric Functions} +% ===================================================================== +\beginscroll +% +\Language{} provides two basic floating-point types: \axiomType{Float} and +\axiomType{DoubleFloat}. This section describes how to use numerical +%-% \HDindex{function!numeric}{ugProblemNumericPage}{8.1.}{Numeric Functions} +operations defined on these types and the related complex types. +%-% \HDindex{numeric operations}{ugProblemNumericPage}{8.1.}{Numeric Functions} +% +As we mentioned in \downlink{``\ugIntroTitle''}{ugIntroPage} in Chapter \ugIntroNumber\ignore{ugIntro}, the \axiomType{Float} type is a software +implementation of floating-point numbers in which the exponent and the +%-% \HDindex{floating-point number}{ugProblemNumericPage}{8.1.}{Numeric Functions} +significand may have any number of digits. +%-% \HDindex{number!floating-point}{ugProblemNumericPage}{8.1.}{Numeric Functions} +See \downlink{`Float'}{FloatXmpPage}\ignore{Float} for detailed information about this domain. +The \axiomType{DoubleFloat} (see \downlink{`DoubleFloat'}{DoubleFloatXmpPage}\ignore{DoubleFloat}) is usually a hardware +implementation of floating point numbers, corresponding to machine double +precision. +The types \axiomType{Complex Float} and \axiomType{Complex DoubleFloat} are +%-% \HDindex{floating-point number!complex}{ugProblemNumericPage}{8.1.}{Numeric Functions} +the corresponding software implementations of complex floating-point numbers. +%-% \HDindex{complex!floating-point number}{ugProblemNumericPage}{8.1.}{Numeric Functions} +In this section the term {\it floating-point type} means any of these +%-% \HDindex{number!complex floating-point}{ugProblemNumericPage}{8.1.}{Numeric Functions} +four types. +% +The floating-point types implement the basic elementary functions. +These include (where \axiomSyntax{\$} means +\axiomType{DoubleFloat}, +\axiomType{Float}, +\axiomType{Complex DoubleFloat}, or +\axiomType{Complex Float}): + +\noindent +\axiomFun{exp}, \axiomFun{log}: \axiom{\$ -> \$} \newline +\axiomFun{sin}, \axiomFun{cos}, \axiomFun{tan}, \axiomFun{cot}, \axiomFun{sec}, \axiomFun{csc}: \axiom{\$ -> \$} \newline +\axiomFun{sin}, \axiomFun{cos}, \axiomFun{tan}, \axiomFun{cot}, \axiomFun{sec}, \axiomFun{csc}: \axiom{\$ -> \$} \newline +\axiomFun{asin}, \axiomFun{acos}, \axiomFun{atan}, \axiomFun{acot}, \axiomFun{asec}, \axiomFun{acsc}: \axiom{\$ -> \$} \newline +\axiomFun{sinh}, \axiomFun{cosh}, \axiomFun{tanh}, \axiomFun{coth}, \axiomFun{sech}, \axiomFun{csch}: \axiom{\$ -> \$} \newline +\axiomFun{asinh}, \axiomFun{acosh}, \axiomFun{atanh}, \axiomFun{acoth}, \axiomFun{asech}, \axiomFun{acsch}: \axiom{\$ -> \$} \newline +\axiomFun{pi}: \axiom{() -> \$} \newline +\axiomFun{sqrt}: \axiom{\$ -> \$} \newline +\axiomFun{nthRoot}: \axiom{(\$, Integer) -> \$} \newline +\axiomFunFrom{**}{Float}: \axiom{(\$, Fraction Integer) -> \$} \newline +\axiomFunFrom{**}{Float}: \axiom{(\$,\$) -> \$} \newline + +The handling of roots depends on whether the floating-point type +%-% \HDindex{root!numeric approximation}{ugProblemNumericPage}{8.1.}{Numeric Functions} +is real or complex: for the real floating-point types, +\axiomType{DoubleFloat} and \axiomType{Float}, if a real root exists +the one with the same sign as the radicand is returned; for the +complex floating-point types, the principal value is returned. +%-% \HDindex{principal value}{ugProblemNumericPage}{8.1.}{Numeric Functions} +Also, for real floating-point types the inverse functions +produce errors if the results are not real. +This includes cases such as \axiom{asin(1.2)}, \axiom{log(-3.2)}, +\axiom{sqrt(-1.1)}. +% +\xtc{ +The default floating-point type is \axiomType{Float} so to evaluate +functions using \axiomType{Float} or \axiomType{Complex Float}, just use +normal decimal notation. +}{ +\spadpaste{exp(3.1)} +} +\xtc{ +}{ +\spadpaste{exp(3.1 + 4.5 * \%i)} +} +\xtc{ +To evaluate functions using \axiomType{DoubleFloat} +or \axiomType{Complex DoubleFloat}, +a declaration or conversion is required. +}{ +\spadpaste{r: DFLOAT := 3.1; t: DFLOAT := 4.5; exp(r + t*\%i)} +} +\xtc{ +}{ +\spadpaste{exp(3.1::DFLOAT + 4.5::DFLOAT * \%i)} +} +% +A number of special functions are provided by the package +\axiomType{DoubleFloatSpecialFunctions} for the machine-precision +%-% \HDindex{special functions}{ugProblemNumericPage}{8.1.}{Numeric Functions} +floating-point types. +%-% \HDexptypeindex{DoubleFloatSpecialFunctions}{ugProblemNumericPage}{8.1.}{Numeric Functions} +The special functions provided are listed below, where \axiom{F} stands for +the types \axiomType{DoubleFloat} and \axiomType{Complex DoubleFloat}. +The real versions of the functions yield an error if the result is not real. +%-% \HDindex{function!special}{ugProblemNumericPage}{8.1.}{Numeric Functions} + +\noindent +\axiomFun{Gamma}: \axiom{F -> F}\hfill\newline +\axiom{Gamma(z)} is the Euler gamma function, +%-% \HDindex{function!Gamma}{ugProblemNumericPage}{8.1.}{Numeric Functions} + \texht{$\Gamma(z)$}{\axiom{Gamma(z)}}, + defined by +%-% \HDindex{Euler!gamma function}{ugProblemNumericPage}{8.1.}{Numeric Functions} + \texht{\narrowDisplay{\Gamma(z) = \int_{0}^{\infty} t^{z-1} e^{-t} dt.}% + }{ +\newline +\centerline{{\axiom{Gamma(z) = integrate(t**(z-1)*exp(-t), t=0..\%infinity).}}} + } + +\noindent +\axiomFun{Beta}: \axiom{F -> F}\hfill\newline + \axiom{Beta(u, v)} is the Euler Beta function, +%-% \HDindex{function!Euler Beta}{ugProblemNumericPage}{8.1.}{Numeric Functions} + \texht{$B(u,v)$}{\axiom{B(u,v)}}, defined by +%-% \HDindex{Euler!Beta function}{ugProblemNumericPage}{8.1.}{Numeric Functions} + \texht{\narrowDisplay{B(u,v) = \int_{0}^{1} t^{u-1} (1-t)^{v-1} dt.}% + }{ +\newline +\centerline{{ \axiom{Beta(u,v) = integrate(t**(u-1)*(1-t)**(b-1), t=0..1).}}} + } + This is related to \texht{$\Gamma(z)$}{\axiom{Gamma(z)}} by + \texht{\narrowDisplay{B(u,v) = \frac{\Gamma(u) \Gamma(v)}{\Gamma(u + v)}.}% + }{ +\newline +\centerline{{ \axiom{Beta(u,v) = Gamma(u)*Gamma(v)/Gamma(u + v).}}} + }% + +\noindent +\axiomFun{logGamma}: \axiom{F -> F}\hfill\newline + \axiom{logGamma(z)} is the natural logarithm of +\texht{$\Gamma(z)$}{\axiom{Gamma(z)}}. + This can often be computed even if \texht{$\Gamma(z)$}{\axiom{Gamma(z)}} +cannot. +% + +\noindent +\axiomFun{digamma}: \axiom{F -> F}\hfill\newline + \axiom{digamma(z)}, also called \axiom{psi(z)}, +%-% \HDindex{psi @ $\psi$}{ugProblemNumericPage}{8.1.}{Numeric Functions} +is the function \texht{$\psi(z)$}{\axiom{psi(z)}}, +%-% \HDindex{function!digamma}{ugProblemNumericPage}{8.1.}{Numeric Functions} + defined by \texht{\narrowDisplay{\psi(z) = \Gamma'(z)/\Gamma(z).}% + }{ +\newline +\centerline{{ \axiom{psi(z) = Gamma'(z)/Gamma(z).}}} + }% + +\noindent +\axiomFun{polygamma}: \axiom{(NonNegativeInteger, F) -> F}\hfill\newline + \axiom{polygamma(n, z)} is the \eth{\axiom{n}} derivative of +%-% \HDindex{function!polygamma}{ugProblemNumericPage}{8.1.}{Numeric Functions} + \texht{$\psi(z)$}{\axiom{digamma(z)}}\texht{, written $\psi^{(n)}(z)$}{}. + +\noindent +\axiomFun{besselJ}: \axiom{(F,F) -> F}\hfill\newline + \axiom{besselJ(v,z)} is the Bessel function of the first kind, +%-% \HDindex{function!Bessel}{ugProblemNumericPage}{8.1.}{Numeric Functions} + \texht{$J_\nu (z)$}{\axiom{J(v,z)}}. + This function satisfies the differential equation + \texht{\narrowDisplay{z^2 w''(z) + z w'(z) + (z^2-\nu^2)w(z) = 0.}% + }{ +\newline +\centerline{{ \axiom{z**2*w''(z) + z*w'(z) + (z**2-v**2)*w(z) = 0.}}} + }% + +\noindent +\axiomFun{besselY}: \axiom{(F,F) -> F}\hfill\newline + \axiom{besselY(v,z)} is the Bessel function of the second kind, +%-% \HDindex{function!Bessel}{ugProblemNumericPage}{8.1.}{Numeric Functions} + \texht{$Y_\nu (z)$}{\axiom{Y(v,z)}}. + This function satisfies the same differential equation as + \axiomFun{besselJ}. + The implementation simply uses the relation + \texht{\narrowDisplay{Y_\nu (z) = \frac{J_\nu (z) \cos(\nu \pi) - J_{-\nu} (z)}{\sin(\nu \pi)}.}% + }{ +\newline +\centerline{{ \axiom{Y(v,z) = (J(v,z)*cos(v*\%pi) - J(-v,z))/sin(v*\%pi).}}} + }% + +\noindent +\axiomFun{besselI}: \axiom{(F,F) -> F}\hfill\newline + \axiom{besselI(v,z)} is the modified Bessel function of the first kind, +%-% \HDindex{function!Bessel}{ugProblemNumericPage}{8.1.}{Numeric Functions} + \texht{$I_\nu (z)$}{\axiom{I(v,z)}}. + This function satisfies the differential equation + \texht{\narrowDisplay{z^2 w''(z) + z w'(z) - (z^2+\nu^2)w(z) = 0.}% + }{ +\newline +\centerline{{ \axiom{z**2 * w''(z) + z * w'(z) - (z**2+v**2)*w(z) = 0.}}} + }% + +\noindent +\axiomFun{besselK}: \axiom{(F,F) -> F}\hfill\newline + \axiom{besselK(v,z)} is the modified Bessel function of the second kind, +%-% \HDindex{function!Bessel}{ugProblemNumericPage}{8.1.}{Numeric Functions} + \texht{$K_\nu (z)$}{\axiom{K(v,z)}}. + This function satisfies the same differential equation as \axiomFun{besselI}. +%-% \HDindex{Bessel function}{ugProblemNumericPage}{8.1.}{Numeric Functions} + The implementation simply uses the relation + \texht{\narrowDisplay{K_\nu (z) = \pi \frac{I_{-\nu} (z) - I_{\nu} (z)}{2 \sin(\nu \pi)}.}% + }{ +\newline +\centerline{{ \axiom{K(v,z) = \%pi*(I(v,z) - I(-v,z))/(2*sin(v*\%pi)).}}} + } + +\noindent +\axiomFun{airyAi}: \axiom{F -> F}\hfill\newline + \axiom{airyAi(z)} is the Airy function \texht{$Ai(z)$}{\axiom{Ai(z)}}. +%-% \HDindex{function!Airy Ai}{ugProblemNumericPage}{8.1.}{Numeric Functions} + This function satisfies the differential equation + \texht{$w''(z) - z w(z) = 0.$}{\axiom{w''(z) - z * w(z) = 0.}} + The implementation simply uses the relation + \texht{\narrowDisplay{Ai(-z) = \frac{1}{3}\sqrt{z} ( J_{-1/3} (\frac{2}{3}z^{3/2}) + J_{1/3} (\frac{2}{3}z^{3/2}) ).}% + }{ +\newline +\centerline{{ \axiom{Ai(-z) = 1/3 * sqrt(z) * (J(-1/3, 2/3*z**(3/2)) + J(1/3, 2/3*z**(3/2)) ).}}} + }% + +\noindent +\axiomFun{airyBi}: \axiom{F -> F}\hfill\newline + \axiom{airyBi(z)} is the Airy function \texht{$Bi(z)$}{\axiom{Bi(z)}}. +%-% \HDindex{function!Airy Bi}{ugProblemNumericPage}{8.1.}{Numeric Functions} + This function satisfies the same differential equation as \axiomFun{airyAi}. +%-% \HDindex{Airy function}{ugProblemNumericPage}{8.1.}{Numeric Functions} + The implementation simply uses the relation + \texht{\narrowDisplay{Bi(-z) = \frac{1}{3}\sqrt{3 z} ( J_{-1/3} (\frac{2}{3}z^{3/2}) - J_{1/3} (\frac{2}{3}z^{3/2}) ).}% + }{ +\newline +\centerline{{ \axiom{Bi(-z) = 1/3 *sqrt(3*z) * (J(-1/3, 2/3*z**(3/2)) - J(1/3, 2/3*z**(3/2)) ).}}} + } + +\noindent +\axiomFun{hypergeometric0F1}: \axiom{(F,F) -> F}\hfill\newline + \axiom{hypergeometric0F1(c,z)} is the hypergeometric function +%-% \HDindex{function!hypergeometric}{ugProblemNumericPage}{8.1.}{Numeric Functions} + \texht{${}_0 F_1 ( ; c; z)$}{\axiom{0F1(; c; z)}}.% + +\xtc{ +The above special functions are defined only for small floating-point types. +If you give \axiomType{Float} arguments, they are converted to +\axiomType{DoubleFloat} by \Language{}. +}{ +\spadpaste{Gamma(0.5)**2} +} +\xtc{ +}{ +\spadpaste{a := 2.1; b := 1.1; besselI(a + \%i*b, b*a + 1)} +} +% +A number of additional operations may be used to compute numerical values. +These are special polynomial functions that can be evaluated for values in +any commutative ring \axiom{R}, and in particular for values in any +floating-point type. +The following operations are provided by the package +\axiomType{OrthogonalPolynomialFunctions}: +%-% \HDexptypeindex{OrthogonalPolynomialFunctions}{ugProblemNumericPage}{8.1.}{Numeric Functions} + +\noindent +\axiomFun{chebyshevT}: \axiom{(NonNegativeInteger, R) -> R}\hbox{}\hfill\newline + \axiom{chebyshevT(n,z)} is the \eth{\axiom{n}} Chebyshev polynomial of the first + kind, \texht{$T_n (z)$}{\axiom{T[n](z)}}. These are defined by + \texht{\narrowDisplay{\frac{1-t z}{1-2 t z+t^2} = \sum_{n=0}^{\infty} T_n (z) t^n.}% + }{ +\newline +\centerline{{ \axiom{(1-t*z)/(1-2*t*z+t**2) = sum(T[n](z) *t**n, n = 0..).}}} + }% + +\noindent +\axiomFun{chebyshevU}: \axiom{(NonNegativeInteger, R) -> R}\hbox{}\hfill\newline + \axiom{chebyshevU(n,z)} is the \eth{\axiom{n}} Chebyshev polynomial of the second + kind, \texht{$U_n (z)$}{\axiom{U[n](z)}}. These are defined by + \texht{\narrowDisplay{\frac{1}{1-2 t z+t^2} = \sum_{n=0}^{\infty} U_n (z) t^n.}% + }{ +\newline +\centerline{{ \axiom{1/(1-2*t*z+t**2) = sum(U[n](z) *t**n, n = 0..).}}} + }% + +\noindent +\axiomFun{hermiteH}: \axiom{(NonNegativeInteger, R) -> R}\hbox{}\hfill\newline + \axiom{hermiteH(n,z)} is the \eth{\axiom{n}} Hermite polynomial, + \texht{$H_n (z)$}{\axiom{H[n](z)}}. + These are defined by + \texht{\narrowDisplay{e^{2 t z - t^2} = \sum_{n=0}^{\infty} H_n (z) \frac{t^n}{n!}.}% + }{ +\newline +\centerline{{ \axiom{exp(2*t*z-t**2) = sum(H[n](z)*t**n/n!, n = 0..).}}} + }% + +\noindent +\axiomFun{laguerreL}: \axiom{(NonNegativeInteger, R) -> R}\hbox{}\hfill\newline + \axiom{laguerreL(n,z)} is the \eth{\axiom{n}} Laguerre polynomial, + \texht{$L_n (z)$}{\axiom{L[n](z)}}. + These are defined by + \texht{\narrowDisplay{\frac{e^{-\frac{t z}{1-t}}}{1-t} = \sum_{n=0}^{\infty} L_n (z) \frac{t^n}{n!}.}% + }{ +\newline +\centerline{{ \axiom{exp(-t*z/(1-t))/(1-t) = sum(L[n](z)*t**n/n!, n = 0..).}}} + }% + +\noindent +\axiomFun{laguerreL}: \axiom{(NonNegativeInteger, NonNegativeInteger, R) -> R}\hbox{}\hfill\newline + \axiom{laguerreL(m,n,z)} is the associated Laguerre polynomial, + \texht{$L^m_n (z)$}{\axiom{L[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(b + 1) - B(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} +% -- cgit v1.2.3