Algebraic Numbers

The following sections give a high-level introduction to the way Mathematica handles algebraic numbers. It is not essential to understanding the Mathematica equation solvers and is included for interested readers. The text is adapted from a talk presented by Robert Knapp at the Japanese Society for Synmbolic and Algebraic Computing.

For higher-order polynomials, it is not always possible to represent the exact solution in closed form. Currently, the highest polynomial that a closed form solution has been found for is the Quinoid or 5th order polynomial. However, the expression represents an exact mathematical quantity, which is based on the abstraction of algebraic numbers, represented in Mathematica as Root objects.

The method Mathematica uses internally to calculate roots of polynomials is the well-established Jenkins-Traub method. It works by finding successive roots with what is essentially a numerically stable variant of the Newton-Raphson method. Then the polynomial is deflated by dividing out by a factor involving each root that is found.  This is the engine behind algebraic numbers.  What makes it possible to have an implicit representation that is reasonably efficient are mathematical theorems, which allow roots to be isolated, and Mathematica's significance arithmetic.

Let's take a closer look at the roots of a higher order equation.

Solve[(x^5) - (x^3) +  σ x  +  1  ==  0, x]

This gives the solution to the equation.  This is given in terms of replacement rules so that you can see how the variables relate to the solution.  In this case it indicates that there are five possible solutions.  In cases where you just want to work with the roots, you can extract them into a list.

polyroots  =  x  /.  %

A plot of the root positions in the complex plane as a function of the parameter value helps to visualize the numbers represented here implicitly. First, an auxiliary package is loaded in order to highlight the roots with different colors.

Needs["Graphics`Colors`"];

ParametricPlot[
        Evaluate[{Re[#], Im[#]} &  /@  roots], {σ, -2, 2}, 
        PlotStyle -> {Red,  Blue,  Green,  Black,  Purple}];

It is interesting to note that two roots come together on the real axis at the point σ = -1.  The lines going across the plot are an artifact of the canonical order Mathematica gives to the roots.  You may be asking how Mathematica distinguishes between and orders the roots.  The process is an elegant example of numerical and symbolic mathematics working together.  We will work with the last root to illustrate this. To make things a bit more convenient, we define the root as a function of the parameter.

lastroot[σ_]  =  Last[roots];

When the parameter σ is given a numerical value, Mathematica does the computations necessary to determine an approximate root that is precise enough to be isolated from the other roots of the polynomial.

lastroot[0]

First, root approximations at machine precision are found using the Jenkins-Traub algorithm.  The isolation is possible because, for a given root approximation, it is possible to find an enclosing circle in which an exact root of the polynomial is guaranteed to lie.  Isolation makes sure that none of the circles around the approximate values overlap. Optionally, Mathematica allows you to use an exact root isolation algorithm that is more complicated and time consuming.  See [Strzebonski, 1996] for details on both methods of isolation.  When σ = 0, the roots are reasonably far apart from each other, so that 16 digits is sufficient to ensure isolation.  Once the roots are isolated, they are ordered based on a canonical ordering.  Real roots are ordered first along the real line, then complex ones based on complex absolute value.  This is why you see the lines across the plot above at σ = -1.  The pairs of roots switch to maintain the canonical ordering.

When σ = -1, the polynomial is reducible in terms of radicals. Mathematica recognizes this, so the polynomial is given in its simplest form.

lastroot[-1]

However, for σ sightly greater than -1, the polynomial cannot be reduced, and so it is given in terms of a root object.

lastroot[-1 + (10^(-30))]

Notice that the polynomial is recalculated to be the minimal polynomial with integer coefficients.  Internally, Mathematica recognizes that the two roots are very close together and adapts to using higher precision (in this case, 32 digits) to do the root isolation for these.  Then, when you request a numerical value, the positive imaginary part is indicated correctly.

N[%]

Hopefully this example has given you a rough idea of how the fusion of implicit representations with numerical approximations using adaptive precision yields a powerful computational tool.  Mathematica provides a framework where objects like algebraic numbers fit naturally. You can similarly use Mathematica to set up implicit representations of your own objects by defining numeric and symbolic rules that apply to them. Implicit representations can be thought of, in some sense, as a transparent interface to a combination of powerful numerical algorithms and mathematical information without overburdening a user with the implementation details.