8 Advanced Problem Solving
In this chapter we describe techniques useful in solving advanced problems with Axiom.
8.1 Numeric Functions
Axiom provides two basic floating-point types: Float and DoubleFloat. This section describes how to use numerical function:numeric operations defined on these types and the related complex types. numeric operations
As we mentioned in Chapter ugIntro , the Float type is a software implementation of floating-point numbers in which the exponent and the floating-point number significand may have any number of digits. number:floating-point See FloatXmpPage for detailed information about this domain. The DoubleFloat (see DoubleFloatXmpPage ) is usually a hardware implementation of floating point numbers, corresponding to machine double precision. The types Complex Float and Complex DoubleFloat are floating-point number:complex the corresponding software implementations of complex floating-point numbers. complex:floating-point number In this section the term floating-point type means any of these number:complex floating-point four types.
The floating-point types implement the basic elementary functions. These include (where $ means DoubleFloat, Float, Complex DoubleFloat, or Complex Float):
exp, log: $ -> $
sin, cos, tan, cot, sec, csc: $ -> $
asin, acos, atan, acot, asec, acsc: $ -> $
sinh, cosh, tanh, coth, sech, csch: $ -> $
asinh, acosh, atanh, acoth, asech, acsch: $ -> $
pi: () -> $
sqrt: $ -> $
nthRoot: ($, Integer) -> $
****Float: ($, Fraction Integer) -> $
****Float: ($,$) -> $
The handling of roots depends on whether the floating-point type root:numeric approximation is real or complex: for the real floating-point types, DoubleFloat and 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. principal value Also, for real floating-point types the inverse functions produce errors if the results are not real. This includes cases such as , , .
The default floating-point type is Float so to evaluate functions using Float or Complex Float, just use normal decimal notation.
To evaluate functions using DoubleFloat or Complex DoubleFloat, a declaration or conversion is required.
A number of special functions are provided by the package DoubleFloatSpecialFunctions for the machine-precision special functions floating-point types. DoubleFloatSpecialFunctions The special functions provided are listed below, where stands for the types DoubleFloat and Complex DoubleFloat. The real versions of the functions yield an error if the result is not real. function:special
Gamma:
is the Euler gamma function,
function:Gamma
,
defined by
Euler:gamma function
Beta:
is the Euler Beta function,
function:Euler Beta
, defined by
Euler:Beta function
This is related to by
logGamma:
is the natural logarithm of
.
This can often be computed even if
cannot.
digamma:
, also called ,
psi @
is the function
function:digamma
defined by
polygamma:
is the -th derivative of
function:polygamma
, written .
besselJ:
is the Bessel function of the first kind,
function:Bessel
.
This function satisfies the differential equation
besselY:
is the Bessel function of the second kind,
function:Bessel
.
This function satisfies the same differential equation as
besselJ.
The implementation simply uses the relation
besselI:
is the modified Bessel function of the first kind,
function:Bessel
.
This function satisfies the differential equation
besselK:
is the modified Bessel function of the second kind,
function:Bessel
.
This function satisfies the same differential equation as besselI.
Bessel function
The implementation simply uses the relation
airyAi:
is the Airy function .
function:Airy Ai
This function satisfies the differential equation
The implementation simply uses the relation
airyBi:
is the Airy function .
function:Airy Bi
This function satisfies the same differential equation as airyAi.
Airy function
The implementation simply uses the relation
hypergeometric0F1:
is the hypergeometric function
function:hypergeometric
.
The above special functions are defined only for small floating-point types. If you give Float arguments, they are converted to DoubleFloat by Axiom.
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 , and in particular for values in any floating-point type. The following operations are provided by the package OrthogonalPolynomialFunctions: OrthogonalPolynomialFunctions
chebyshevT:
is the -th Chebyshev polynomial of the first
kind, . These are defined by
chebyshevU:
is the -th Chebyshev polynomial of the second
kind, . These are defined by
hermiteH:
is the -th Hermite polynomial,
.
These are defined by
laguerreL:
is the -th Laguerre polynomial,
.
These are defined by
laguerreL:
is the associated Laguerre polynomial
.
This is the -th derivative of .
legendreP:
is the -th Legendre polynomial,
. These are defined by
These operations require non-negative integers for the indices, but otherwise the argument can be given as desired.
The expression evaluates to the -th Chebyshev polynomial:Chebyshev:of the first kind polynomial of the first kind.
The expression evaluates to the -th Chebyshev polynomial:Chebyshev:of the second kind polynomial of the second kind.
The expression evaluates to the -th Hermite polynomial:Hermite polynomial.
The expression evaluates to the -th Laguerre polynomial:Laguerre polynomial.
The expression polynomial:Legendre evaluates to the -th Legendre polynomial,
Finally, three number-theoretic polynomial operations may be evaluated. number theory The following operations are provided by the package NumberTheoreticPolynomialFunctions. NumberTheoreticPolynomialFunctions.
bernoulliB:
is the -th Bernoulli polynomial,
polynomial:Bernoulli
. These are defined by
eulerE:
is the -th Euler polynomial,
Euler:polynomial
. These are defined by
polynomial:Euler
cyclotomic:
is the -th cyclotomic polynomial
. This is the polynomial whose
roots are precisely the primitive -th roots of unity.
Euler:totient function
This polynomial has degree given by the Euler totient function
function:totient
.
The expression evaluates to the -th Bernoulli polynomial:Bernouilli polynomial.
The expression polynomial:Euler evaluates to the -th Euler polynomial.
The expression polynomial:cyclotomic evaluates to the -th cyclotomic polynomial. cyclotomic polynomial
Drawing complex functions in Axiom is presently somewhat awkward compared to drawing real functions. It is necessary to use the draw operations that operate on functions rather than expressions.
This is the complex exponential function (rotated interactively). function:complex exponential 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.
This is the complex arctangent function. function:complex arctangent 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.
This is the complex Gamma function.
This shows the real Beta function near the origin.
This is the Bessel function for index in the range and argument in the range .
This is the modified Bessel function evaluated for various real values of the index and fixed argument .
This is similar to the last example except the index takes on complex values in a rectangle centered on the origin.
8.2 Polynomial Factorization
The Axiom polynomial factorization polynomial:factorization facilities are available for all polynomial types and a wide variety of coefficient domains. factorization Here are some examples.
8.2.1 Integer and Rational Number Coefficients
Polynomials with integer polynomial:factorization:integer coefficients coefficients can be be factored.
Also, Axiom can factor polynomials with polynomial:factorization:rational number coefficients rational number coefficients.
8.2.2 Finite Field Coefficients
Polynomials with coefficients in a finite field polynomial:factorization:finite field coefficients can be also be factored. finite field:factoring polynomial with coefficients in
These include the integers mod , where is prime, and extensions of these fields.
Convert this to have coefficients in the finite field with elements. See ugProblemFinite for more information about finite fields.
8.2.3 Simple Algebraic Extension Field Coefficients
Polynomials with coefficients in simple algebraic extensions polynomial:factorization:algebraic extension field coefficients of the rational numbers can be factored. algebraic number number:algebraic
Here, and are symbolic roots of polynomials.
Note that the second argument to factor can be a list of algebraic extensions to factor over.
This factors over the integers.
Factor the same polynomial over the field obtained by adjoining to the rational numbers.
Factor over the same field.
Factor again over the field obtained by adjoining both and to the rational numbers.
8.2.4 Factoring Rational Functions
Since fractions of polynomials form a field, every element (other than zero) rational function:factoring divides any other, so there is no useful notion of irreducible factors. Thus the factor operation is not very useful for fractions of polynomials.
There is, instead, a specific operation factorFraction that separately factors the numerator and denominator and returns a fraction of the factored results.
You can also use map. This expression applies the factor operation to the numerator and denominator.
8.3 Manipulating Symbolic Roots of a Polynomial
In this section we show you how to work with one root or all roots root:symbolic of a polynomial. These roots are represented symbolically (as opposed to being numeric approximations). See ugxProblemOnePol and ugxProblemPolSys for information about solving for the roots of one or more polynomials.
8.3.1 Using a Single Root of a Polynomial
Use rootOf to get a symbolic root of a polynomial: returns a root of .
This creates an algebraic number . algebraic number number:algebraic
To find the algebraic relation that defines , use definingPolynomial.
You can use in any further expression, including a nested rootOf.
Higher powers of the roots are automatically reduced during calculations.
The operation zeroOf is similar to rootOf, except that it may express the root using radicals in some cases. radical
8.3.2 Using All Roots of a Polynomial
Use rootsOf to get all symbolic roots of a polynomial: returns a list of all the roots of . If has a multiple root of order , then that root root:multiple appears times in the list. Make sure these variables are x0 etc.
Compute all the roots of .
As a side effect, the variables and are bound to the first three roots of .
Although they all satisfy and are different algebraic numbers. To find the algebraic relation that defines each of them, use definingPolynomial.
We can check that the sum and product of the roots of are its trace and norm.
Corresponding to the pair of operations rootOf/ zeroOf in ugxProblemOnePol , there is an operation zerosOf that, like rootsOf, computes all the roots of a given polynomial, but which expresses some of them in terms of radicals.
As you see, only one implicit algebraic number was created ( ), and its defining equation is this. The other three roots are expressed in radicals.
8.4 Computation of Eigenvalues and Eigenvectors
In this section we show you some of Axiom's facilities for computing and eigenvalue manipulating eigenvalues and eigenvectors, also called eigenvector characteristic values and characteristic vectors, characteristic:value respectively. characteristic:vector
Let's first create a matrix with integer entries.
To get a list of the rational eigenvalues, use the operation eigenvalues.
Given an explicit eigenvalue, eigenvector computes the eigenvectors corresponding to it.
The operation eigenvectors returns a list of pairs of values and vectors. When an eigenvalue is rational, Axiom 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 Axiom to solve the minimal polynomial, then you can substitute these roots polynomial:minimal into the parametric form of the corresponding eigenvectors. minimal polynomial
You must be aware that unless an exact eigenvalue has been computed, the eigenvector may be badly in error.
Another possibility is to use the operation radicalEigenvectors tries to compute explicitly the eigenvectors in terms of radicals. radical
Alternatively, Axiom can compute real or complex approximations to the approximation eigenvectors and eigenvalues using the operations realEigenvectors or complexEigenvectors. They each take an additional argument to specify the ``precision'' required. precision In the real case, this means that each approximation will be within of the actual result. In the complex case, this means that each approximation will be within of the actual result in each of the real and imaginary parts.
The precision can be specified as a Float if the results are desired in floating-point notation, or as Fraction Integer if the results are to be expressed using rational (or complex rational) numbers.
If an by matrix has distinct eigenvalues (and therefore eigenvectors) the operation eigenMatrix gives you a matrix of the eigenvectors.
If a symmetric matrix matrix:symmetric has a basis of orthonormal eigenvectors, then basis:orthonormal orthonormalBasis computes a list of these vectors. orthonormal basis
8.5 Solution of Linear and Polynomial Equations
In this section we discuss the Axiom facilities for solving systems of linear equations, finding the roots of polynomials and linear equation solving systems of polynomial equations. For a discussion of the solution of differential equations, see ugProblemDEQ .
8.5.1 Solution of Systems of Linear Equations
You can use the operation solve to solve systems of linear equations. equation:linear:solving
The operation 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.
To solve the linear system: evaluate this expression.
Parameters are given as new variables starting with a percent sign and % and the variables are expressed in terms of the parameters. If the system has no solutions then the empty list is returned.
When you solve the linear system with this expression you get a solution involving a parameter.
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.
To solve the system: in matrix form you would evaluate this expression.
The solutions are presented as a Record with two components: the component particular contains a particular solution of the given system or the item "failed" if there are no solutions, the component 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 basis component contains non-trivial vectors.
This happens when you solve the linear system with this command.
All solutions of this system are obtained by adding the particular solution with a linear combination of the basis vectors.
When no solution exists then "failed" is returned as the particular component, as follows:
When you want to solve a system of homogeneous equations (that is, a system where the numbers on the right-hand sides of the nullspace equations are all zero) in the matrix form you can omit the second argument and use the nullSpace operation.
This computes the solutions of the following system of equations: The result is given as a list of vectors and these vectors form a basis for the solution space.
8.5.2 Solution of a Single Polynomial Equation
Axiom can solve polynomial equations producing either approximate polynomial:root finding or exact solutions. equation:polynomial:solving Exact solutions are either members of the ground field or can be presented symbolically as roots of irreducible polynomials.
This returns the one rational root along with an irreducible polynomial describing the other solutions.
If you want solutions expressed in terms of radicals you would use this instead. radical
The solve command always returns a value but radicalSolve returns only the solutions that it is able to express in terms of radicals. radical
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'' precision . This means that each approximation will be within of the actual result.
Notice that the type of second argument controls the type of the result.
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.
If you want approximate complex results you should use the approximation command complexSolve that takes the same precision argument .
Each approximation will be within of the actual result in each of the real and imaginary parts.
Note that if you omit the = from the first argument Axiom 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.
Axiom can also solve equations involving rational functions. Solutions where the denominator vanishes are discarded.
8.5.3 Solution of Systems of Polynomial Equations
Given a system of equations of rational functions with exact coefficients: equation:polynomial:solving
Axiom 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.
Symbolic solutions can be presented using ``implicit'' algebraic numbers defined as roots of irreducible polynomials or in terms of radicals. Axiom can also find approximations to the real or complex roots of a system of polynomial equations to any user-specified accuracy.
The operation solve for systems is used in a way similar to 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 ugxProblemOnePol .
Use the operation solve if you want implicitly presented solutions.
Use radicalSolve if you want your solutions expressed in terms of radicals.
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.
If the precision is expressed as a floating-point number you get results expressed as floats.
To get complex numeric solutions, use the operation complexSolve, which takes the same arguments as in the real case.
It is also possible to solve systems of equations in rational functions over the rational numbers. Note that is not returned as a solution since the denominator vanishes there.
When solving equations with denominators, all solutions where the denominator vanishes are discarded.
8.6 Limits
To compute a limit, you must specify a functional expression, limit a variable, and a limiting value for that variable. If you do not specify a direction, Axiom attempts to compute a two-sided limit.
Issue this to compute the limit
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, limit:one-sided vs. two-sided you can compute a one-sided limit.
The function is only defined to the right of zero, that is, for . Thus, when computing limits of functions involving , you probably want a ``right-hand'' limit.
When you do not specify `` '' or `` '' as the optional fourth argument, limit tries to compute a two-sided limit. Here the limit from the left does not exist, as Axiom indicates when you try to take a two-sided limit.
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 is simply the absolute value of , the function is simply the sign ( or ) of the nonzero real number . Therefore, for and for .
This is what happens when we take the limit at . The answer returned by Axiom gives both a ``left-hand'' and a ``right-hand'' limit.
Here is another example, this time using a more complicated function.
You can compute limits at infinity by passing either limit:at infinity or as the third argument of limit.
To do this, use the constants and .
You can take limits of functions with parameters. limit:of function with parameters As you can see, the limit is expressed in terms of the parameters.
When you use limit, you are taking the limit of a real function of a real variable.
When you compute this, Axiom returns because, as a function of a real variable, is always between and , so tends to as tends to .
However, as a function of a complex variable, is badly limit:real vs. complex behaved near (one says that has an essential singularity essential singularity at ). singularity:essential
When viewed as a function of a complex variable, does not approach any limit as tends to in the complex plane. Axiom indicates this when we call complexLimit.
Here is another example. As approaches along the real axis, tends to .
However, if is allowed to approach along any path in the complex plane, the limiting value of depends on the path taken because the function has an essential singularity at . This is reflected in the error message returned by the function.
You can also take complex limits at infinity, that is, limits of a function of as approaches infinity on the Riemann sphere. Use the symbol to denote ``complex infinity.''
As above, to compute complex limits rather than real limits, use complexLimit.
In many cases, a limit of a real function of a real variable exists when the corresponding complex limit does not. This limit exists.
But this limit does not.
8.7 Laplace Transforms
Axiom can compute some forward Laplace transforms, mostly Laplace transform of elementary function:elementary functions transform:Laplace not involving logarithms, although some cases of special functions are handled.
To compute the forward Laplace transform of with respect to and express the result as , issue the command .
Here are some other non-trivial examples.
Axiom also knows about a few special functions.
When Axiom does not know about a particular transform, it keeps it as a formal transform in the answer.
8.8 Integration
Integration is the reverse process of differentiation, that is, integration an integral of a function with respect to a variable is any function such that is equal to .
The package FunctionSpaceIntegration provides the top-level integration operation, integrateintegrateFunctionSpaceIntegration, for integrating real-valued elementary functions. FunctionSpaceIntegration
Unfortunately, antiderivatives of most functions cannot be expressed in terms of elementary functions.
Given an elementary function to integrate, Axiom 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.
Similar functions may have antiderivatives antiderivative that look quite different because the form of the antiderivative depends on the sign of a constant that appears in the function.
If the integrand contains parameters, then there may be several possible antiderivatives, depending on the signs of expressions of the parameters.
In this case Axiom returns a list of answers that cover all the possible cases. Here you use the answer involving the square root of when and integration:result as list of real functions the answer involving the square root of when .
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 complexIntegratecomplexIntegrateFunctionSpaceComplexIntegration, which is provided by the package integration:result as a complex functions FunctionSpaceComplexIntegration. FunctionSpaceComplexIntegration
This operation is used for integrating complex-valued elementary functions.
As with the real case, antiderivatives for most complex-valued functions cannot be expressed in terms of elementary functions.
Sometimes integrate can involve symbolic algebraic numbers such as those returned by rootOfrootOfExpression. To see how to work with these strange generated symbols (such as ), see ugxProblemSymRootAll .
Definite integration is the process of computing the area between integration:definite the -axis and the curve of a function . The fundamental theorem of calculus states that if is continuous on an interval and if there exists a function that is differentiable on and such that is equal to , then the definite integral of for in the interval is equal to .
The package RationalFunctionDefiniteIntegration provides the top-level definite integration operation, integrateintegrateRationalFunctionDefiniteIntegration, for integrating real-valued rational functions.
Axiom checks beforehand that the function you are integrating is defined on the interval , and prints an error message if it finds that this is not case, as in the following example:
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.
When parameters are present in the function, the function may or may not be defined on the interval of integration.
If this is the case, Axiom issues a warning that a pole might lie in the path of integration, and does not compute the integral.
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 ``noPole'' as a third argument to integrateintegrateRationalFunctionDefiniteIntegration:
The value here is, of course, incorrect if is between and
8.9 Working with Power Series
Axiom has very sophisticated facilities for working with power series series. 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 )set streams calculate. For the purposes of this book, we have used this system command to display fewer than ten terms. set streams calculate 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 series, although you can also use taylor, laurent and puiseux in situations where you know what kind of exponents are involved.
For information about solving differential equations in terms of power series, see ugxProblemDEQSeries .
8.9.1 Creation of Power Series
This is the easiest way to create a power series. This tells Axiom that is to be treated as a power series, series:creating so functions of are again power series.
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.
Here the coefficients are integers (note that the coefficients are the Fibonacci Fibonacci numbers numbers).
This series has coefficients that are rational numbers.
When you enter this expression you introduce the symbolic constants and
When you enter the expression the variable appears in the resulting series expansion.
You can also convert an expression into a series expansion. This expression creates the series expansion of about . For details and more examples, see ugxProblemSeriesConversions .
You can create power series with more general coefficients. You normally accomplish this via a type declaration (see ugTypesDeclare ). See ugxProblemSeriesFunctions for some warnings about working with declared series.
We declare that is a one-variable Taylor series series:Taylor (UTS is the abbreviation for UnivariateTaylorSeries) in the variable with FLOAT (that is, floating-point) coefficients, centered about Then, by assignment, we obtain the Taylor expansion of with floating-point coefficients. UnivariateTaylorSeries
You can also create a power series by giving an explicit formula for its -th coefficient. For details and more examples, see ugxProblemSeriesFormula .
To create a series about whose -th Taylor coefficient is , you can evaluate this expression. This is the Taylor expansion of at .
8.9.2 Coefficients of Power Series
You can extract any coefficient from a power series---even one that hasn't been computed yet. This is possible because in Axiom, 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 lazy evaluation.) When you ask for a series:lazy evaluation coefficient that hasn't yet been computed, Axiom computes lazy evaluation whatever additional coefficients it needs and then stores them in the representation of the power series.
Here's an example of how to extract the coefficients of a power series. series:extracting coefficients
This coefficient is readily available.
But let's get the fifteenth coefficient of .
If you look at then you see that the coefficients up to order have all been computed.
8.9.3 Power Series Arithmetic
You can manipulate power series using the usual arithmetic operations series:arithmetic , , , and (from UnivariatePuiseuxSeries)
The results of these operations are also power series.
You can also compute , where and are two power series.
8.9.4 Functions on Power Series
Once you have created a power series, you can apply transcendental functions (for example, exp, log, sin, tan, cosh, etc.) to it.
To demonstrate this, we first create the power series expansion of the rational function
about .
If you want to compute the series expansion of
you simply compute the sine of .
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.
If you evaluate then you have declared that is a one variable Taylor series series:Taylor (UTS is the abbreviation for UnivariateTaylorSeries) in the variable with FRAC INT (that is, fractions of integer) coefficients, centered about .
You can now compute certain power series in , provided that these series have rational coefficients.
You can get examples of such series by applying transcendental functions to series in that have no constant terms.
Similarly, you can compute the logarithm of a power series with rational coefficients if the constant coefficient is
If you wanted to apply, say, the operation exp to a power series with a nonzero constant coefficient , then the constant coefficient of the result would be , which is not a rational number. Therefore, evaluating would generate an error message.
If you want to compute the Taylor expansion of , you must ensure that the coefficient domain has an operation exp defined for it. An example of such a domain is Expression Integer, the type of formal functional expressions over the integers.
When working with coefficients of this type,
this presents no problems.
Another way to create Taylor series whose coefficients are expressions over the integers is to use taylor which works similarly to series:Taylor series.
This is equivalent to the previous computation, except that now we are using the variable instead of .
8.9.5 Converting to Power Series
The ExpressionToUnivariatePowerSeries package provides operations for computing series expansions of functions. ExpressionToUnivariatePowerSeries
Evaluate this to compute the Taylor expansion of about series:Taylor . The first argument, , specifies the function whose series expansion is to be computed and the second argument, , specifies that the series is to be expanded in power of , that is, in power of .
Here is the Taylor expansion of about :
The function to be expanded into a series may have variables other than series:multiple variables the series variable.
For example, we may expand as a Taylor series in
or as a Taylor series in .
A more interesting function is When we expand this function as a Taylor series in the -th order coefficient is the -th Bernoulli Bernoulli:polynomial polynomial polynomial:Bernoulli divided by .
Therefore, this and the next expression produce the same result.
Technically, a series with terms of negative degree is not considered to be a Taylor series, but, rather, a series:Laurent Laurent series. Laurent series If you try to compute a Taylor series expansion of at via you get an error message. The reason is that the function has a pole at , 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.
You get the desired series expansion by issuing this.
Similarly, a series with terms of fractional degree is neither a Taylor series nor a Laurent series. Such a series is called a series:Puiseux Puiseux series. Puiseux series The expression results in an error message because the series expansion about this point has terms of fractional degree.
However, this command produces what you want.
Finally, consider the case of functions that do not have Puiseux expansions about certain points. An example of this is about . produces an error message because of the type of singularity of the function at .
The general function series can be used in this case. Notice that the series returned is not, strictly speaking, a power series because of the in the expansion.
The operation 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.
8.9.6 Power Series from Formulas
The GenerateUnivariatePowerSeries package enables you to series:giving formula for coefficients create power series from explicit formulas for their -th 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. GenerateUnivariatePowerSeries See ugxProblemSeriesConversions for more information.
Consider the Taylor expansion of series:Taylor about :
The -th Taylor coefficient is .
This is how you create this series in Axiom.
The first argument specifies a formula for the -th coefficient by giving a function that maps to . The second argument specifies that the series is to be expanded in powers of , that is, in powers of . 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 ugUserAnon .
Consider the Taylor expansion of about :
If you were to evaluate the expression you would get an error message because Axiom would try to calculate a term of degree and therefore divide by
Instead, evaluate this. The third argument, , indicates that only terms of degree are to be computed.
Next consider the Taylor expansion of an odd function, say, :
Here every other coefficient is zero and we would like to give an explicit formula only for the odd Taylor coefficients.
This is one way to do it. The third argument, , specifies that the first term to be computed is the term of degree 1. The fourth argument, , specifies that we increment by to find the degrees of subsequent terms, that is, the next term is of degree , the next of degree , etc.
The initial degree and the increment do not have to be integers. For example, this expression produces a series expansion of .
While the increment must be positive, the initial degree may be negative. This yields the Laurent expansion of at . (bernoulli(numer(n+1)) is necessary because bernoulli takes integer arguments.)
Of course, the reciprocal of this power series is the Taylor expansion of .
As a final example,here is the Taylor expansion of about .
When we compute the of this series, we get (in the sense that all higher terms computed so far are zero).
Axiom isn't sufficiently ``symbolic'' in the sense we might wish. It is an open problem to decide that ``x'' is the only surviving term. Two attacks on the problem might be:
(1) Notice that all of the higher terms are identically zero but Axiom can't decide that from the information it knows. Presumably we could attack this problem by looking at the sin function as a taylor series around x=0 and seeing the term cancellation occur. This uses a term-difference mechanism.
(2) Notice that there is no way to decide that the stream for asinx is actually the definition of asin(x). But we could recognize that the stream for asin(x) has a generator term and so will a taylor series expansion of sin(x). From these two generators it may be possible in certain cases to decide that the application of one generator to the other will yield only ``x''. This trick involves finding the correct inverse for the stream functions. If we can find an inverse for the ``remaining tail'' of the stream we could conclude cancellation and thus turn an infinite stream into a finite object.
In general this is the zero-equivalence problem and is undecidable.
As we discussed in ugxProblemSeriesConversions , you can also use the operations taylor, laurent and puiseux instead of series if you know ahead of time what kind of exponents a series has. You can't go wrong using series, though.
8.9.7 Substituting Numerical Values in Power Series
Use evalevalUnivariatePowerSeriesCategory approximation to substitute a numerical value for a variable in series:numerical approximation a power series. For example, here's a way to obtain numerical approximations of from the Taylor series expansion of .
First you create the desired Taylor expansion.
Then you evaluate the series at the value . The result is a sequence of the partial sums.
8.9.8 Example: Bernoulli Polynomials and Sums of Powers
Axiom provides operations for computing definite and summation:definite indefinite sums. summation:indefinite
You can compute the sum of the first ten fourth powers by evaluating this. This creates a list whose entries are as ranges from 1 to 10, and then computes the sum of the entries of that list.
You can also compute a formula for the sum of the first fourth powers, where is an unspecified positive integer.
This formula is valid for any positive integer . For instance, if we replace by 10, summation:definite we obtain the number we computed earlier.
You can compute a formula for the sum of the first -th powers in a similar fashion. Just replace the in the definition of sum4 by any expression not involving . Axiom computes these formulas using Bernoulli polynomials; Bernoulli:polynomial we polynomial:Bernoulli use the rest of this section to describe this method.
First consider this function of and .
Since the expressions involved get quite large, we tell Axiom to show us only terms of degree up to
If we look at the Taylor expansion of about we see that the coefficients of the powers of are polynomials in .
Type: UnivariateTaylorSeries(Expression Integer,t,0)
In fact, the -th coefficient in this series is essentially the -th Bernoulli polynomial: the -th coefficient of the series is , where is the -th Bernoulli polynomial. Thus, to obtain the -th Bernoulli polynomial, we multiply the -th coefficient of the series by .
For example, the sixth Bernoulli polynomial is this.
We derive some properties of the function . First we compute .
If we normalize , we see that it has a particularly simple form.
From this it follows that the -th coefficient in the Taylor expansion of at is .
If you want to check this, evaluate the next expression.
However, since it follows that the -th coefficient is Equating coefficients, we see that and, therefore,
Let's apply this formula repeatedly, letting vary between two integers and , with :
When we add these equations we find that the sum of the left-hand sides is the sum of the powers from to . The sum of the right-hand sides is a ``telescoping series.'' After cancellation, the sum is simply
Replacing by , we have shown that
Let's use this to obtain the formula for the sum of fourth powers.
First we obtain the Bernoulli polynomial .
To find the sum of the first 4th powers, we multiply by .
This is the same formula that we obtained via .
At this point you may want to do the same computation, but with an exponent other than For example, you might try to find a formula for the sum of the first 20th powers.
8.10 Solution of Differential Equations
In this section we discuss Axiom's facilities for equation:differential:solving solving differential equation differential equations in closed-form and in series.
Axiom provides facilities for closed-form solution of equation:differential:solving in closed-form single differential equations of the following kinds:
- linear ordinary differential equations, and
- non-linear first order ordinary differential equations when integrating factors can be found just by integration.
For a discussion of the solution of systems of linear and polynomial equations, see ugProblemLinPolEqn .
8.10.1 Closed-Form Solutions of Linear Differential Equations
A differential equation is an equation involving an unknown function and one or more of its derivatives. differential equation The equation is called ordinary if derivatives with respect to equation:differential only one dependent variable appear in the equation (it is called partial otherwise). The package ElementaryFunctionODESolver provides the top-level operation solve for finding closed-form solutions of ordinary differential equations. ElementaryFunctionODESolver
To solve a differential equation, you must first create an operator for operator the unknown function.
We let be the unknown function in terms of .
You then type the equation using D to create the derivatives of the unknown function where is any symbol you choose (the so-called dependent variable).
This is how you enter the equation .
The simplest way to invoke the solve command is with three arguments.
- the differential equation,
- the operator representing the unknown function,
- the dependent variable.
So, to solve the above equation, we enter this.
Since linear ordinary differential equations have infinitely many solutions, solve returns a particular solution and a basis for the solutions of the corresponding homogenuous equation. Any expression of the form where the 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 initial conditions: choose a value for the dependent variable and specify the values of the unknown function and its derivatives at . 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 solve tries to find it. To specify initial conditions to solve, use an Equation of the form for the third parameter instead of the dependent variable, and add a fourth parameter consisting of the list of values .
To find the solution of satisfying , do this.
You can omit the when you enter the equation to be solved.
Axiom 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, Axiom is not limited by the order of the equation.
Axiom can solve the following third order equations with polynomial coefficients.
Here we are solving a homogeneous equation.
On the other hand, and in contrast with the operation integrate, it can happen that Axiom 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:
- If the coefficients are constants, Axiom finds a complete basis of solutions (i,e, all solutions).
- If the coefficients are rational functions in the dependent variable, Axiom at least finds all solutions that do not involve algebraic functions.
Note that this last statement does not mean that Axiom does not find the solutions that are algebraic functions. It means that it is not guaranteed that the algebraic function solutions will be found.
This is an example where all the algebraic solutions are found.
8.10.2 Closed-Form Solutions of Non-Linear Differential Equations
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 .
Using the notation , we have and .
We first check for exactness, that is, does ?
This is not zero, so the equation is not exact. Therefore we must look for an integrating factor: a function such that . Normally, we first search for depending only on or only on .
Let's search for such a first.
If the above is zero for a function that does not depend on , then is an integrating factor.
The solution depends on , so there is no integrating factor that depends on only.
Let's look for one that depends on only.
We've found one!
The above is an integrating factor. We must multiply our initial equation (that is, and ) by the integrating factor.
Let's check for exactness.
We must solve the exact equation, that is, find a function such that and .
We start by writing where is an unknown function of . This guarantees that .
All we want is to find such that .
The above particular solution is the we want, so we just replace by it in the implicit solution.
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 solve operation.
First define to be an operator.
Next we create the differential equation.
Finally, we solve it.
8.10.3 Power Series Solutions of Differential Equations
The command to solve differential equations in power equation:differential:solving in power series series power series around series:power a particular initial point with specific initial conditions is called seriesSolve. It can take a variety of parameters, so we illustrate its use with some examples.
Since the coefficients of some solutions are quite large, we reset the default to compute only seven terms.
You can solve a single nonlinear equation of any order. For example, we solve subject to
We first tell Axiom that the symbol denotes a new operator.
Enter the differential equation using like any system function.
Solve it around with the initial conditions .
You can also solve a system of nonlinear first order equations. For example, we solve a system that has and as solutions.
We tell Axiom that is also an operator.
Enter the two equations forming our system.
Solve the system around with the initial conditions and . Notice that since we give the unknowns in the order , the answer is a list of two series in the order
Compiling function %BZ with type List UnivariateTaylorSeries(
Expression Integer,t,0) -> UnivariateTaylorSeries(Expression
Integer,t,0)
Compiling function %CA with type List UnivariateTaylorSeries(
Expression Integer,t,0) -> UnivariateTaylorSeries(Expression
Integer,t,0)
The order in which we give the equations and the initial conditions has no effect on the order of the solution.
8.11 Finite Fields
A finite field (also called a 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 , called the characteristic, such that
for any element in the finite field. In fact, the
number of elements in a finite field is a power of the characteristic
and for each prime and positive integer there exists exactly
one finite field with elements, up to isomorphism. For
more information about the algebraic structure and properties of
finite fields, see, for example,
S. Lang, Algebra, Second
Edition, New York: Addison-Wesley Publishing Company, Inc., 1984, ISBN
0 201 05487 6;
or
R. Lidl, H. Niederreiter, Finite Fields,
Encyclopedia of Mathematics and Its Applications, Vol. 20, Cambridge:
Cambridge Univ. Press, 1983, ISBN 0 521 30240 4.
When the field has elements and is called a prime field, discussed in the next section. There are several ways of implementing extensions of finite fields, and Axiom 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 ugTypesPkgCall for more information on package-calling.
8.11.1 Modular Arithmetic and Prime Fields
finite field Galois:field field:finite:prime field:prime field:Galois prime field modular arithmetic arithmetic:modular
Let 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 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 modulo or, more simply mod .
In Axiom, you use IntegerMod to do such arithmetic.
If is not prime, there is only a limited notion of reciprocals and division.
There are 12 exposed and 13 unexposed library operations named /
having 2 argument(s) but none was determined to be applicable.
Use HyperDoc Browse, or issue
)display op /
to learn more about the available operations. Perhaps
package-calling the operation or using coercions on the arguments
will allow you to apply the operation.
Cannot find a definition or applicable library operation named /
with argument type(s)
IntegerMod 12
IntegerMod 12
Perhaps you should use "@" to indicate the required return type,
or "$" to specify which version of the function you need.
Here and are relatively prime, so has a multiplicative inverse modulo .
If we take to be a prime number , then taking inverses and, therefore, division are generally defined.
Use PrimeField instead of IntegerMod for prime.
You can also use and for the inverse of .
PrimeField (abbreviation 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), InnerPrimeField (abbreviation 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 PrimeField 2 and its extensions.
In the following examples, we work with the finite field with elements.
Like many domains in Axiom, finite fields provide an operation for returning a random element of the domain.
The element is a primitive element of this field, primitive element element:primitive
in the sense that its powers enumerate all nonzero elements.
If every nonzero element is a power of a primitive element, how do you determine what the exponent is? Use discrete logarithm discreteLog. logarithm:discrete
The order of a nonzero element is the smallest positive integer such .
The order of a primitive element is the defining .
8.11.2 Extensions of Finite Fields
finite field field:finite:extension of
When you want to work with an extension of a finite field in Axiom, you have three choices to make:
- Do you want to generate an extension of the prime field (for example, PrimeField 2) or an extension of a given field?
- 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 Zech logarithm 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.
- 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?
This illustrates one of the most important features of Axiom: 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.
Constructors that automatically generate extensions of the prime field:
FiniteField
FiniteFieldCyclicGroup
FiniteFieldNormalBasis
Constructors that generate extensions of an arbitrary field:
FiniteFieldExtension
FiniteFieldExtensionByPolynomial
FiniteFieldCyclicGroupExtension
FiniteFieldCyclicGroupExtensionByPolynomial
FiniteFieldNormalBasisExtension
FiniteFieldNormalBasisExtensionByPolynomial
Constructors that use a cyclic group representation:
FiniteFieldCyclicGroup
FiniteFieldCyclicGroupExtension
FiniteFieldCyclicGroupExtensionByPolynomial
Constructors that use a normal basis representation:
FiniteFieldNormalBasis
FiniteFieldNormalBasisExtension
FiniteFieldNormalBasisExtensionByPolynomial
Constructors that use an irreducible modulus polynomial representation:
FiniteField
FiniteFieldExtension
FiniteFieldExtensionByPolynomial
Constructors that generate a polynomial for you:
FiniteField
FiniteFieldExtension
FiniteFieldCyclicGroup
FiniteFieldCyclicGroupExtension
FiniteFieldNormalBasis
FiniteFieldNormalBasisExtension
Constructors for which you provide a polynomial:
FiniteFieldExtensionByPolynomial
FiniteFieldCyclicGroupExtensionByPolynomial
FiniteFieldNormalBasisExtensionByPolynomial
These constructors are discussed in the following sections where we collect together descriptions of extension fields that have the same underlying representation.For more information on the implementation aspects of finite fields, see J. Grabmeier, A. Scheerhorn, Finite Fields in AXIOM, Technical Report, IBM Heidelberg Scientific Center, 1992.
If you don't really care about all this detail, just use FiniteField. As your knowledge of your application and its Axiom 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 PrimeField.
8.11.3 Irreducible Modulus Polynomial Representations
All finite field extension constructors discussed in this finite field section field:finite:extension of 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 ugxProblemFiniteExtensionFinite for general information about finite field extensions.
For FiniteField (abbreviation FF) you provide a prime number and an extension degree . This degree can be 1.
Axiom uses the prime field PrimeField(p), here PrimeField 2, and it chooses an irreducible polynomial of degree , here 12, over the ground field.
The objects in the generated field extension are polynomials of degree at most with coefficients in the prime field. The polynomial indeterminate is automatically chosen by Axiom and is typically something like or . These (strange) variables are only for output display; there are several ways to construct elements of this field.
The operation index enumerates the elements of the field extension and accepts as argument the integers from 1 to .
The expression always gives the indeterminate.
You can build polynomials in and calculate in .
Among the available operations are norm and trace.
Since any nonzero element is a power of a primitive element, how do we discover what the exponent is?
The operation discreteLog calculates discrete logarithm the exponent and, logarithm:discrete if it is called with only one argument, always refers to the primitive element returned by primitiveElement.
FiniteFieldExtension (abbreviation FFX) is similar to FiniteField except that the ground-field for FiniteFieldExtension 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.
FiniteFieldExtensionByPolynomial (abbreviation FFP) is similar to FiniteField and FiniteFieldExtension but is more general.
For 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.
8.11.4 Cyclic Group Representations
finite field field:finite:extension of
In every finite field there exist elements whose powers are all the nonzero elements of the field. Such an element is called a primitive element.
In FiniteFieldCyclicGroup (abbreviation FFCG) group:cyclic the nonzero elements are represented by the powers of a fixed primitive element:primitive element primitive element 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 ugxProblemFiniteUtility for examples of how to compute such a polynomial.
To use FiniteFieldCyclicGroup you provide a prime number and an extension degree.
Axiom uses the prime field, here PrimeField 3, as the ground field and it chooses a primitive polynomial of degree , here 4, over the prime field.
You can calculate in .
In this representation of finite fields the discrete logarithm of an element can be seen directly in its output form.
FiniteFieldCyclicGroupExtension (abbreviation FFCGX) is similar to FiniteFieldCyclicGroup except that the ground field for 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.
FiniteFieldCyclicGroupExtensionByPolynomial (abbreviation FFCGP) is similar to FiniteFieldCyclicGroup and FiniteFieldCyclicGroupExtension but is more general. For 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.
We use a utility operation to generate an irreducible primitive polynomial (see ugxProblemFiniteUtility ). The polynomial has one variable that is ``anonymous'': it displays as a question mark.
Let's look at a random element from this field.
8.11.5 Normal Basis Representations
finite field field:finite:extension of basis:normal normal basis
Let be a finite extension of degree of the finite field and let have elements. An element of is said to be normal over if the elements
form a basis of as a vector space over . Such a basis is called a normal basis.This agrees with the general definition of a normal basis because the distinct powers of the automorphism constitute the Galois group of .
If is normal over , its minimal polynomial:minimal polynomial is also said to be normal over . minimal polynomial There exist normal bases for all finite extensions of arbitrary finite fields.
In FiniteFieldNormalBasis (abbreviation FFNB), the elements of the finite field are represented by coordinate vectors with respect to a normal basis.
You provide a prime and an extension degree .
Axiom uses the prime field PrimeField(p), here PrimeField 3, and it chooses a normal polynomial of degree , 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 Axiom and is typically something like or . 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
You can calculate in using .
FiniteFieldNormalBasisExtension (abbreviation FFNBX) is similar to FiniteFieldNormalBasis except that the groundfield for 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.
FiniteFieldNormalBasisExtensionByPolynomial (abbreviation FFNBP) is similar to FiniteFieldNormalBasis and FiniteFieldNormalBasisExtension but is more general. For 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.
We use a utility operation to generate an irreducible normal polynomial (see ugxProblemFiniteUtility ). The polynomial has one variable that is ``anonymous'': it displays as a question mark.
Let's look at a random element from this field.
8.11.6 Conversion Operations for Finite Fields
Let be a finite field.
An extension field of degree over is a subfield of an extension field of degree over if and only if divides .
FiniteFieldHomomorphisms provides conversion operations between different extensions of one fixed finite ground field and between different representations of these finite fields.
Let's choose and ,
build the field extensions,
and pick two random elements from the smaller field.
Since divides , is a subfield of .
Therefore we can convert the elements of into elements of .
To check this, let's do some arithmetic.
There are also conversions available for the situation, when and are represented in different ways (see ugxProblemFiniteExtensionFinite ). For example let's choose where the representation is 0 plus the cyclic multiplicative group and with a normal basis representation.
Check the arithmetic again.
8.11.7 Utility Operations for Finite Fields
FiniteFieldPolynomialPackage (abbreviation FFPOLY) provides operations for generating, counting and testing polynomials over finite fields. Let's start with a couple of definitions:
- A polynomial is primitive if its roots are primitive polynomial:primitive elements in an extension of the coefficient field of degree equal to the degree of the polynomial.
- A polynomial is normal over its coefficient field polynomial:normal if its roots are linearly independent elements in an extension of the coefficient field of degree equal to the degree of the polynomial.
In what follows, many of the generated polynomials have one ``anonymous'' variable. This indeterminate is displayed as a question mark (``?'').
To fix ideas, let's use the field with five elements for the first few examples.
You can generate irreducible polynomials of any (positive) degree polynomial:irreducible (within the storage capabilities of the computer and your ability to wait) by using createIrreduciblePolycreateIrreduciblePolyFiniteFieldPolynomialPackage.
Does this polynomial have other important properties? Use primitive? to test whether it is a primitive polynomial.
Use normal? to test whether it is a normal polynomial.
Note that this is actually a trivial case, because a normal polynomial of degree must have a nonzero term of degree . We will refer back to this later.
To get a primitive polynomial of degree 8 just issue this.
This polynomial is not normal,
but if you want a normal one simply write this.
This polynomial is not primitive!
This could have been seen directly, as the constant term is 1 here, which is not a primitive element up to the factor ( ) raised to the degree of the polynomial.Cf. Lidl, R. & Niederreiter, H., 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. The existence of such polynomials is proved in Lenstra, H. W. & Schoof, R. J., Primitive Normal Bases for Finite Fields, Math. Comp. 48, 1987, pp. 217-231.
If you really need one use either createPrimitiveNormalPolycreatePrimitiveNormalPolyFiniteFieldPolynomialPackage or createNormalPrimitivePolycreateNormalPrimitivePolyFiniteFieldPolynomialPackage.
If you want to obtain additional polynomials of the various types above as given by the create... operations above, you can use the next... operations. For instance, nextIrreduciblePolynextIrreduciblePolyFiniteFieldPolynomialPackage 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.
This is the field with five elements.
Our first example irreducible polynomial, say of degree 3, must be ``greater'' than this.
You can generate it by doing this.
Notice that this polynomial is not the same as the one createIrreduciblePolycreateIrreduciblePolyFiniteFieldPolynomialPackage.
You can step through all irreducible polynomials of degree 8 over the field with 5 elements by repeatedly issuing this.
You could also ask for the total number of these.
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 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.
The package 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.
Take these,
and then we have:
What happened?
Well, for the ordering used in nextPrimitivePolynextPrimitivePolyFiniteFieldPolynomialPackage we use as first criterion a comparison of the constant terms of the polynomials. Analogously, in nextNormalPolynextNormalPolyFiniteFieldPolynomialPackage we first compare the monomials of degree 1 less than the degree of the polynomials (which is nonzero, by an earlier remark).
We don't have to restrict ourselves to prime fields.
Let's consider, say, a field with 16 elements.
We can apply any of the operations described above.
Axiom also provides operations for producing random polynomials of a given degree
or with degree between two given bounds.
FiniteFieldPolynomialPackage2 (abbreviation FFPOLY2) exports an operation rootOfIrreduciblePoly for finding one root of an irreducible polynomial polynomial:root of in an extension field of the coefficient field. The degree of the extension has to be a multiple of the degree of . It is not checked whether actually is irreducible.
To illustrate this operation, we fix a ground field
and then an extension field.
We construct an irreducible polynomial over .
We compute a root of .
and check the result
8.12 Primary Decomposition of Ideals
Axiom provides a facility for the primary decomposition ideal:primary decomposition of primary decomposition of ideal polynomial ideals over fields of characteristic zero. The algorithm is discussed in \cite{gtz:gbpdpi} and works in essentially two steps:
- the problem is solved for 0-dimensional ideals by ``generic'' projection on the last coordinate
- a ``reduction process'' uses localization and ideal quotients to reduce the general case to the 0-dimensional one.
The Axiom constructor PolynomialIdeals represents ideals with coefficients in any field and supports the basic ideal operations, including intersection, sum and quotient. 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.
First consider the ideal generated by (which defines a circle in the -plane) and the ideal generated by (corresponding to the straight lines and .
We find the equations defining the intersection of the two loci. This correspond to the sum of the associated ideals.
We can check if the locus contains only a finite number of points, that is, if the ideal is zero-dimensional.
We can find polynomial relations among the generators ( and are the parametric equations of the knot).
We can compute the primary decomposition of an ideal.
We can intersect back.
We can compute the radical of every primary component.
Their intersection is equal to the radical of the ideal of .
8.13 Computation of Galois Groups
As a sample use of Axiom's algebraic number facilities, group:Galois we compute Galois:group the Galois group of the polynomial .
We would like to construct a polynomial such that the splitting field:splitting field splitting field of is generated by one root of . First we construct a polynomial such that one root of generates the field generated by two roots of the polynomial . (As it will turn out, the field generated by two roots of is, in fact, the splitting field of .)
From the proof of the primitive element theorem we know that if and are algebraic numbers, then the field is equal to for an appropriately chosen integer . In our case, we construct the minimal polynomial of , where and are two roots of . We construct this polynomial using resultant. The main result we need is the following: If is a polynomial with roots and is a polynomial with roots , then the polynomial is a polynomial of degree with roots .
For we use the polynomial . For we use the polynomial . Thus, the polynomial we first construct is .
The roots of are . Of course, there are five pairs with , so is a 5-fold root of .
Let's get rid of this factor.
Factor the polynomial .
We see that has two irreducible factors, each of degree . (The fact that the polynomial has two factors of degree is enough to show that the Galois group of is the dihedral group of order .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 is FR POLY INT, that is, Factored Polynomial Integer. Factored This is a special data type for recording factorizations of polynomials with integer coefficients.
We can access the individual factors using the operation nthFactornthFactorFactored.
Consider the polynomial . This is the minimal polynomial of the difference of two roots of . Thus, the splitting field of contains a subfield of degree . We show that this subfield is, in fact, the splitting field of by showing that factors completely over this field.
First we create a symbolic root of the polynomial . (We replaced by in the polynomial so that our symbolic root would be printed as .)
We next tell Axiom to view as a univariate polynomial in with algebraic number coefficients. This is accomplished with this type declaration.
Factor over the field . (This computation will take some time!)
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 factor, the field over which the polynomial is factored is the field generated by
- the algebraic numbers that appear in the coefficients of the polynomial, and
- the algebraic numbers that appear in a list passed as an optional second argument of the operation.
In our case, the coefficients of are all rational integers and only appears in the list, so the field is simply .
It was necessary to give the list 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.
We have shown that the splitting field of has degree . Since the symmetric group of degree 5 has only one transitive subgroup of order , we know that the Galois group of must be this group, the dihedral group group:dihedral of order . Rather than stop here, we explicitly compute the action of the Galois group on the roots of .
First we assign the roots of as the values of five root variables.
We can obtain an individual root by negating the constant coefficient of one of the factors of .
We can obtain a list of all the roots in this way.
The expression
- coefficient(nthFactor(algFactors, i), 0)}
is the -th root of and the elements of are the -th roots of as ranges from to .
Assign the roots as the values of the variables .
Next we express the roots of as polynomials in . We could obtain these roots by calling the operation factor: factors over . However, this is a lengthy computation and we can obtain the roots of as differences of the roots of . Only ten of these differences are roots of and the other ten are roots of the other irreducible factor of . We can determine if a given value is a root of by evaluating at that particular value. (Of course, the order in which factors are returned by the operation factor is unimportant and may change with different implementations of the operation. Therefore, we cannot predict in advance which differences are roots of and which are not.)
Let's look at four examples (two are roots of and two are not).
Take one of the differences that was a root of and assign it to the variable .
For example, if returned , you would enter this.
Of course, if the difference is, in fact, equal to the root , you should choose another root of .
Automorphisms of the splitting field are given by mapping a generator of the field, namely , to other roots of its minimal polynomial. Let's see what happens when is mapped to .
We compute the images of the roots under this automorphism:
Of course, the values are simply a permutation of the values .
Let's find the value of (execute as many of the following five commands as necessary).
Proceeding in this fashion, you can find the values of . You have represented the automorphism as a permutation of the roots . If you wish, you can repeat this computation for all the roots of and represent the Galois group of as a subgroup of the symmetric group on five letters.
Here are two other problems that you may attack in a similar fashion:
- Show that the Galois group of is the dihedral group of order eight. group:dihedral (The splitting field of this polynomial is the Hilbert class field Hilbert class field of field:Hilbert class the quadratic field .)
- Show that the Galois group of has order 6 and is isomorphic to the symmetric group on three letters. group:symmetric (The splitting field of this polynomial is the splitting field of .)
8.14 Non-Associative Algebras and Modelling Genetic Laws
Many algebraic structures of mathematics and Axiom have a multiplication operation * that satisfies the associativity law associativity law for all , and . The octonions are a well known exception. There are many other interesting non-associative structures, such as the class of Lie algebra Lie algebras.Two Axiom implementations of Lie algebras are LieSquareMatrix and FreeNilpotentLie. Lie algebras can be used, for example, to analyse Lie symmetry algebras of symmetry partial differential differential equation:partial equations. partial differential equation In this section we show a different application of non-associative algebras, non-associative algebra the modelling of genetic laws. algebra:non-associative
The Axiom library contains several constructors for creating non-associative structures, ranging from the categories Monad, NonAssociativeRng, and FramedNonAssociativeAlgebra, to the domains AlgebraGivenByStructuralConstants and GenericNonAssociativeAlgebra. Furthermore, the package AlgebraPackage provides operations for analysing the structure of such algebras. The interested reader can learn more about these aspects of the Axiom 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
The implementation of general algebras in Axiom allows us to Mendel's genetic laws use this as the definition for multiplication in an algebra. genetics Hence, it is possible to study questions of genetic inheritance using Axiom. To demonstrate this more precisely, we discuss one example from a monograph of A. Wörz-Busekros, where you can also find a general setting of this theory. Wörz-Busekros, A., 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 and gives zygotes zygote , which produce new gametes. gamete In classical Mendelian segregation we have . In general, we have
The segregation rates are the structural constants of an -dimensional algebra. This is provided in Axiom by the constructor AlgebraGivenByStructuralConstants (abbreviation ALGSC).
Consider two coupled autosomal loci with alleles , , , and , building four different gametes and { and }. The zygotes produce gametes and with classical Mendelian segregation. Zygote undergoes transition to and vice versa with probability .
Define a list of four four-by-four matrices giving the segregation rates. We use the value for .
Choose the appropriate symbols for the basis of gametes,
Define the algebra.
What are the probabilities for zygote to produce the different gametes?
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 and is described by their product .
This product is commutative only if the gametes are not sex-dependent, as in our example.
In general, it is not associative.
Random mating within a population is described by . The next generation is .
Use decimal numbers to compare the distributions more easily.
To compute directly the gametic distribution in the fifth generation, we use plenaryPower.
We now ask two questions: Does this distribution converge to an equilibrium state? What are the distributions that are stable?
This is an invariant of the algebra and it is used to answer the first question. The new indeterminates describe a symbolic distribution.
Because the coefficient has absolute value less than 1, all distributions do converge, by a theorem of this theory.
The second question is answered by searching for idempotents in the algebra.
Solve these equations and look at the first solution.
Further analysis using the package PolynomialIdeals shows that there is a two-dimensional variety of equilibrium states and all other solutions are contained in it.
Choose one equilibrium state by setting two indeterminates to concrete values.
Verify the result.