Algebraic numbers in Mathematica 3.0

Complex algebraic numbers form a new class of "infinite precision" numbers available in Mathematica 3.0.

by AdamW.Strzebonski

Version 3.0 of Mathematica represents algebraic numbers using "Root objects," for example:

Root[7 #(1^5) + #(1^3) - 11&, 1]

The first element of a Root object is a pure function representing the minimal polynomial multiplied by the least common denominator of its coefficients. The second element is the root number : an integer between one and the degree of the minimal polynomial specifying which of the roots we have in mind.

To determine the ordering of roots of the minimal polynomial Mathematica isolates all the roots of the polynomial (i.e. finds disjoint subsets of the complex plane, each containing exactly one root), and then orders the isolating sets. The information about the isolating sets is stored with a Root object in its third, hidden, element.

Mathematica uses two types of isolating sets. The type used in current computation depends on the setting of the option ExactRootIsolation in Root.With the default setting ExactRootIsolation -> False, the isolating sets are circles in the complex plane represented by arbitrary-precision number approximations to the roots. To find the isolating circles Mathematica uses a validated numeric method based on the Jenkins-Traub algorithm  [Jenkins 1969], [Ralston and Rabinowitz 1978].

When ExactRootIsolation -> True the isolating sets used are open intervals with rational endpoints for real roots, and open rectangles in the complex plane (Cartesian products of open intervals with rational endpoints) for complex roots. To obtain the isolating sets Mathematica avoids methods involving computations with inexact numbers. For real root isolation we use the method described in [Akritas, Bocharov and Strzebonski, 1994], based on the Descartes' Rule of Signs (see also [Arkitas 1980], [Akritas and Collins 1976]). Our complex root isolation algorithm is based on [Geddes, Czapor and Labahn 1992].

Root isolation

The root isolation procedures used by Root when ExactRootIsolation -> True are available after loading the package Algebra`RootIsolation`.

<< "Algebra`RootIsolation`"

The function RealRootIntervals isolates the real roots of 7 (x^9) - 3 x + 1. There are three real roots, one in each of the open intervals (-1, 0), (0, (1/2)), and ((1/2), 1).

RealRootIntervals[7*x^9 - 3*x + 1]

The function ComplexRootIntervals isolates the complex roots.

ComplexRootIntervals[7*x^9 - 3*x + 1]

We can use the function showIsolatingSets to display the complex roots of 7 (x^9) - 3 x + 1. The isolating sets are the open intervals and open rectangles. The red dots mark approximations of roots computed with NSolve. (In the graphics below, the size of the dots is exaggerated for display purposes).

showIsolatingSets[f_] := Show[Graphics[({RGBColor[Random[]/2, Random[], Random[]], 
        Thickness[0.01], Line[{{#1[[1]], 0}, {#1[[2]], 0}}]} & ) \
/@ Select[ComplexRootIntervals[f], FreeQ[#1, Complex] & ]], 
   Prolog -> {({RGBColor[Random[]/2, Random[], Random[]], 
         Rectangle[{Re[#1[[1]]], Im[#1[[1]]]}, 
          {Re[#1[[2]]], Im[#1[[2]]]}]} & ) /@ 
      Select[ComplexRootIntervals[f],  !FreeQ[#1, Complex] & ]}, 
   Epilog -> {RGBColor[1, 0, 0], PointSize[0.03],
     Point[{Re[x], Im[x]}] /. NSolve[f == 0, x]}, 
   Axes -> True, PlotRange -> {{-2, 2}, {-2, 2}},
   AspectRatio -> 1];

showIsolatingSets[7*x^9  -  3*x  +  1]

We can simultaneously isolate roots of more then one polynomial (i.e. get disjoint sets such that each contains exactly one root of exactly one of the polynomials). The result is a pair of lists of equal length. The first is a list of isolating intervals, the second is a list of numbers showing which root of each input polynomial is contained in the corresponding interval.

RealRootIntervals[{(x^2) - 2, (x^2) - 3}]

The function CountRoots can be used to show the number of roots of a polynomial within a designated open rectangle. For example, there are five 21-st roots of unity in the open rectangle
 (0, 1)×(0, 1).

CountRoots[(x^21) - 1, {x,  0, 1 + I}]

We can see these five roots—in fact, all the roots—of the polynomial using our function showIsolatingSets.

showIsolatingSets[x^21  -  1];

Ordering the Roots

The ordering of roots of the minimal polynomial is defined as follows. For ExactRootIsolation -> False, we order the approximations so that the real numbers come first in natural order, and then the pairs of complex conjugate numbers ordered lexicographically with respect to the roots with negative imaginary part. For ExactRootIsolation -> True, the real roots come first in natural order, and then the pairs of complex conjugate roots ordered lexicographically with respect to the coordinates of the lower left corner of the isolating rectangle of the root with the negative imaginary part. The ordering used guarantees that the real roots always have the lowest root numbers and are ordered in accordance with their actual values. The ordering of complex roots cannot be well-described in terms of their actual values, however the conjugate roots are always adjacent.

Root[f, k] gives the k-th root of the polynomial f. To order roots of polynomials that are not irreducible, Mathematica first orders factors (in the order given by FactorList) and then orders roots of each factor as described above.

Root[(((#(1^6) - 3 #1 + 15))^2) (((#(1^5) - 7 #(1^2) - 11))^3)&, 17]

Roots of all polynomials of degree less than five can be represented in terms of radicals. Algebraic numbers with minimal polynomials of degree less then three are automatically converted to their radical representation.

Root[#(1^2) + #1 + 1&, 2]

ToRadicals converts algebraic numbers with minimal polynomials of degree three and four into their radical representations. For algebraic numbers with minimal polynomials of higher degree ToRadicals attempts to find a radical representation, but is not guaranteed to succeed, even if such representation exists.

Root[4 + #(1^5) + (((-3 + #(1^5)))^3), 1]

ToRadicals[%]

Computation with algebraic numbers

Algebraic numbers form an algebraically closed field containing the field of rational numbers. This means that sums and products of algebraic numbers, and multiplicative inverses of nonzero algebraic numbers are algebraic numbers. Moreover, every univariate polynomial of degree n with algebraic number coefficients has exactly n roots that are algebraic numbers.

Since the degree of the minimal polynomial of a sum or product of algebraic numbers grows very quickly (it is usually equal to the product of degrees of the minimal polynomials of the algebraic numbers), Mathematica does not transform sums, products, and rational powers of 
algebraic numbers into single Root objects automatically. To do this we can use the function RootReduce.

Root[#(1^5) + 3 #1 - 7&, 3]/2 + (Sqrt[11])

RootReduce[%]

((Root[#(1^7) - 3 #(1^2) + 9&, 5]^2)^(1/3))

RootReduce[%]

RootReduce represents algebraic numbers in a canonical way. Algebraic numbers with minimal polynomials of degrees higher than two are represented by single Root objects. Algebraic numbers of degree ≤2 are given a uniques representation in radicals. Here we use RootReduce to show that a rather complicated radical expression is in fact equal to (Sqrt[2]).

-(Sqrt[2 (185 + 21 (Sqrt[21]))]) + 
        (Sqrt[2 (186 + 21 (Sqrt[21]) + 2 (Sqrt[185 + 21 (Sqrt[21])]))]) // RootReduce

Let us now briefly describe the algorithms used by RootReduce. For a more detailed description see [Strzebonski 1996]. If algebraic numbers a and b have minimal polynomials f and g, then the following table gives polynomials divisible by the minimal polynomials of a+b, a-b, ab, a/b, (a^(p/q)) :

a + b     Resultant[f(x-y), g(y), y]
a - b     Resultant[f(x+y), g(y), y]
a · b     Resultant[(y^(Degree(f)))f((x/y)), g(y), y]
a ÷ b     Resultant[f(x·y), g(y), y]
(a^(p/q))       Resultant[f(y), (x^q)  -  (y^p), y]

To decide which root of the polynomial given in the table is the result c of the arithmetic operation, RootReduce factors the polynomial and isolates roots of factors. Knowing isolating sets of a and b,  we can find a region in the complex plane that contains c.  By making all the isolating sets smaller if necessary, we make sure that only one of the isolating sets of the roots of factors intersects this region. This root must be c.

The complex conjugate, the real and imaginary parts and the absolute value of an algebraic number are algebraic numbers. While the complex conjugate of an algebraic number is a root of the same minimal polynomial as an adjoining root number, the  real and imaginary parts and the absolute value of an algebraic number usually have much more complicated minimal polynomials. If the degree of minimal polynomial of a is n, then usually the minimal polynomial of Re(a) has degree n(n-1)/2 and the minimal polynomials of Im(a), and Abs(a) have degree n(n-1).

a = Root[#(1^5) - 7 #1 - 3&, 4];

Conjugate[a]

Re[a]

Im[a]

Abs[a]

Solving polynomial equations

All roots of polynomials with rational number coefficients are algebraic numbers, so Mathematica can solve arbitrary degree polynomial equations with rational number coefficients.

Solve[(x^5) - 11 x + (7/2) == 0, x]

Roots of polynomials with algebraic number coefficients are also algebraic numbers, so Mathematica can solve arbitrary degree polynomial equations with algebraic number coefficients.

Solve[(x^5) - (Sqrt[3]) (x^3) + 2 Root[#(1^3) - 2 #1 + 7&, 1] == 0, x]

If a system of multivariate polynomial equations has a finite number of solutions, then the solutions can be represented using algebraic numbers. Here are the last 3 (out of the total of 30) solutions of a system of three equations.

Solve[{x + (y^2) == 1, y + (z^3) == 1, z + (x^5) == 1}, {x, y, z}]
        [[3]]

Polynomial operations over algebraic extensions of rationals

The result of the Mathematica polynomial operations: PolynomialGCD, PolynomialLCM, Together, Cancel, FactorSquareFree, FactorSquareFreeList, Factor, and FactorList depends on the coefficient field over which the computation is performed. The coefficient field used depends on the value of  the option Extension. With the default value None the coefficient field is taken to be the field of rational numbers, and algebraic numbers are treated as independent variables.

PolynomialGCD[(x^5) - 2 x + 7, x - Root[#(1^5) - 2 #1 + 7&, 3]]

Cancel[(((x^2) - 2 (Sqrt[5]) x + 5)/((x^2) - 5))]

Factor[7 (x^2) - (Sqrt[7]) (x^3) - (Sqrt[7]) (y^2) + x (y^2)]

With Extension -> Automatic the coefficient field used is the field of rational numbers extended by the algebraic numbers that appear in the coefficients of the input.

PolynomialGCD[(x^5) - 2 x + 7, x - Root[#(1^5) - 2 #1 + 7&, 3], 
      Extension -> Automatic]

Cancel[(((x^2) - 2 (Sqrt[5]) x + 5)/((x^2) - 5)), 
      Extension -> Automatic]

Factor[7 (x^2) - (Sqrt[7]) (x^3) - (Sqrt[7]) (y^2) + x (y^2), 
      Extension -> Automatic]

When Extension is set to a list of algebraic numbers, the coefficient field used is the filed of rational numbers extended by the algebraic numbers that appear in the coefficients of the input or on the list. Extending the coefficient field by algebraic numbers that do not appear in the coefficients of the input may change only the result of Factor and FactorList.

Factor[7 (x^2) - (Sqrt[7]) (x^3) - (Sqrt[7]) (y^2) + x (y^2), 
      Extension -> {(7^(1/4))}]

To compute the polynomial greatest common divisor (g.c.d.) over algebraic number fields Mathematica uses the Euclidean algorithm [Geddes, Czapor, and Labahn 1992]. When computing the g.c.d. of polynomials of relatively low degree over an extension of rationals by many algebraic numbers, we use a zero test algorithm working over a multiple algebraic extension each time we need to decide whether a coefficient is zero. Otherwise we switch to computations over a single algebraic extension by computing a primitive element of the extension [Loos 1982]. To compute the factorization over an algebraic extension of rationals, Mathematica first finds a primitive element for the extension, and then uses the algorithm described in [Geddes, Czapor, and Labahn 1992].

Algebraic functions

When the coefficients of the polynomial f in Root[f, k] contain variables, Root[f, k] represents an "algebraic function": the k-th root of f.

r = Root[#(1^5) - 7 #1 + a&, 3]

r /. a -> 0

Replacing variables with rational numbers turns an algebraic function into an algebraic number. Here, following the replacement, the minimal polynomial became reducible, and we get a root of one of the factors.

r /. a -> (Sqrt[3])

If we replace variables with algebraic numbers, an algebraic function also turns into an algebraic number which is automatically computed as a single Root object.

Symmetric functions of roots

Often a result from a computation (for example an integration or an infinite sum) can be represented as a sum over all roots of a polynomial  function evaluated at each root. In Mathematica such sums can be represented as RootSum objects.

\(RootSum[#1\^5 - 3\ #1 - 7&, E\^#1&]\)

This represents a sum of exponents of all roots of the polynomial #(1^5) - 3 #1 - 7.

N[%]

Here is the numeric value of the sum.

RootSum[#(1^3) - a #1 + 9&, ArcTan[x - #1]&]

Both the polynomial and the function to sum over the roots of the polynomial may contain symbolic parameters.

(∂_x) %

The derivative of ArcTan[x - #1]with respect to x is a rational function of #1. Whenever the second argument of a RootSum is a rational function, the sum is computed explicitly using the following fact.

If ((f(x))/(g(x)))is a rational function, p(x) is a polynomial, and R(z) = (Subscript[a, 0]) (z^n) + (Subscript[a, 1]) (z^(n - 1)) + ... + (Subscript[a, n]) is the resultant of g(x) z - f(x) and p(x) with respect to x, then (∑_(p(r) = 0 ))((f(r))/(g(r)))= -((Subscript[a, 1])/(Subscript[a, 0])).

RootSum also computes an explicit result in some logarithmic cases.

RootSum[#(1^7) - 3 #(1^2) - 11&, Log[2 #1 + 3]&]

In the following examples  an integral and an infinite sum are represented using RootSum.

∫((x (E^x))/((x^5) - 3 x + 1))\[DifferentialD]x

(∑_(k = 1)^∞) ((k^3)/((k^7) - (k^2) - 1))

Symmetric functions of roots of the same polynomial can be represented in terms of coefficients of the minimal polynomial. This follows from Vieta's formulas and the following Fundamental Theorem of Symmetric Polynomials [Becker and Weispfenning 1993].

Every symmetric polynomial in (Subscript[x, 1]),..., (Subscript[x, n]) can be represented as a polynomial in elementary symmetric polynomials:

(Subscript[s, 1]) = (Subscript[x, 1]) + (Subscript[x, 2]) +...+ (Subscript[x, n])
(Subscript[s, 2]) = (Subscript[x, 1]) (Subscript[x, 2]) + (Subscript[x, 1]) (Subscript[x, 3]) +...+ (Subscript[x, n - 1]) (Subscript[x, n])
...
(Subscript[s, n]) = (Subscript[x, 1]) (Subscript[x, 2])·...·(Subscript[x, n])

When the ordering of variables is fixed, every polynomial can be uniquely represented as a sum of its "symmetric part" and the "remainder": f((Subscript[x, 1]),..., (Subscript[x, n])) = p((Subscript[s, 1]),..., (Subscript[s, n])) + q((Subscript[x, 1]),..., (Subscript[x, n])). The polynomial f is symmetric if and only if the remainder q is zero. The uniqueness of this representation is guaranteed by the requirement that q does not contain descending monomials, where a monomial c (Subscript[(Subscript[x, 1]^e), 1])·...·(Subscript[(Subscript[x, n]^e), n]) is called descending iff (Subscript[e, 1])>=...>=(Subscript[e, n]).

Simplify and RootReduce use the Fundamental Theorem of Symmetric Polynomials and Vieta's formulas to simplify expressions containing Root objects.

Simplify[
      Root[#(1^3) + a #1 + b&, 1] Root[#(1^3) + a #1 + b&, 2] + 
        Root[#(1^3) + a #1 + b&, 2] Root[#(1^3) + a #1 + b&, 3] + 
        Root[#(1^3) + a #1 + b&, 3] Root[#(1^3) + a #1 + b&, 1] - 
        2 Root[#(1^3) + a #1 + b&, 1] - 2 Root[#(1^3) + a #1 + b&, 2]]

Functions for generating elementary symmetric polynomials, and finding the symmetric part and the remainder of a polynomial are available after loading the package Algebra`SymmetricPolynomials`.

<< Algebra`SymmetricPolynomials`

\(\(	SymmetricPolynomial[{x, y, z, t}, 3]\)\)

Here is the third elementary symmetric polynomial in four variables.

SymmetricReduction[
      (x^2) + (y^2) + (z^2) - 3 x y - 3 x z - 5 y z, {x, y, z}]

This gives the symmetric part and the remainder of the polynomial (x^2) + (y^2) + (z^2) - 3 x y - 3 x z - 5 y z.

SymmetricReduction[
      (x^2) + (y^2) + (z^2) - 3 x y - 3 x z - 5 y z, {x, y, z}, {(Subscript[s, 1]), 
        (Subscript[s, 2]), (Subscript[s, 3])}]

Note that we can specify the variables to represent the elementary symmetric polynomials in SymmetricReduction.

References

Akritas, A.G.,  An implementation of Vincent's theorem.  Num.Math. 36, 53-62, 1980.

Akritas, A.G., Bocharov, A.V., Strzebonski, A.W. , Implementation of real root isolation algorithms in Mathematica. Abstracts of the Interval'94 Conference, St.Petersburg, Russia, 23-27, 1994.

Akritas, A.G., Collins, G.E.,  Polynomial real root isolation using Descartes' rule of signs. Proceedings of the ACM Symposium on Symbolic and Algebraic Computations, Yorktown Heights, N.Y., 272-275, 1976.

Becker, T., Weispfenning, V., Groebner Bases. A Computational Approach to Commutative Algebra, Springer-Verlag, 1993.

Collins, G.E., Krandick, W., An efficient algorithm for infallible polynomial complex root isolation. Proceedings of the ACM Symposium on Symbolic and Algebraic Computations, Berkeley, California, 189-194, 1992.

Geddes, K.O., Czapor, S.L., Labahn, G., Algorithms for Computer Algebra. Kluwer Academic Publishers, 1992.

Jenkins, M.A.G. ,  Three-stage variable-shift iterations for the solution of polynomial equations with a posteriori error bounds for the zeros. Ph.D. dissertation, Stanford University, 1969.

Loos, R. , Computing in algebraic extensions. Computing Suppl. 4, 173-187, 1982.

Ralston, A., Rabinowitz, P. (1978). A First Course in Numerical Analysis. McGraw-Hill, Inc, 1978.

Strzebonski, W. A., Algebraic numbers in Mathematica 3.0., The Mathematica Journal, 6:4, 74-80, 1996.

ABOUT THE AUTHOR

Adam Strzebonski is a mathematician with Wolfram Research, Inc., 
the developers of Mathematica, and also a mathematics faculty 
member at the Jagiellonian University (Krakow, Poland). He received 
his M.S. from the Jagiellonian University, and his Ph.D from the 
State University of New York at Buffalo. His research  areas in 
mathematics include real and complex algebraic and analytic geometry.

Michael Strzebonski
Wolfram Research, Inc.
100 Trade Center Drive
Champaign, IL. 61820