Polynomials for Radicals

Given a radical expression, it is not hard to discover a polynomial having it as a root.

by Stan Wagon

An algebraic number is a real or complex number that is a root of a polynomial with integer coefficients. They are most familiar when given in terms of radicals. For example, the following are all algebraic numbers:
                                                                                          √(1 + √17)
                     √2             √(3 +  √(3 + √5/2) )           --------------
                                                                                            21/3 + i
A natural problem in symbolic computation is to obtain, for a given radical expression α, a polynomial p(x) in Z[x] such that p(α) = 0. Of course, it is even better if we can get the minimal polynomial for α (irreducible polynomial). As one application, one can perhaps get a simplification of the original radical expression.

Most of the functionality described in this column is in fact built into version 3 of Mathematica. Still, it is illuminating to see how this problem can be solved from first principles, and the pattern-matching ideas that we discuss will be useful in a variety of situations.

Our goal is to define a function, Polynomial, that returns a polynomial for a given radical expression; here are two examples.

Polynomial[(Sqrt[2]) + (Sqrt[3])]

Polynomial[((Sqrt[1 + (Sqrt[17])])/((2^(1/3)) + I))]

The complete Polynomial code (only about 20 lines) can be pasted together from the pieces in this column, or can be found in the electronic supplement.

Heads

Every Mathematica object has a "head", and understanding them is essential to serious symbol manipulation. Here are several examples.

\(\(\(examples
       = {0, 1, 3\/5, 4.3, I, 5 + I, a + b\ I, 5 + \@2\ I, \@2, 1\/\@2, 3\^4, 
        x\^4, Sin, Plus, Head, p, π, Magritte, "I'm not a string"}; 
    \)
\)\)

Most of the heads seem quite reasonable, but there are some subtleties. For example, Complex[3, 7] is the full form of 3 + 7 I; the arguments are the real and imaginary parts. But note that a + 5 I is not a Complex object because the real and imaginary parts are not known: it is not assumed that a is real. This interpretation extends to symbolic expressions such as 5 + Sqrt[2] I. The heads of these objects are Plus, since they have form Plus[u, v]. Another subtlety concerns Sqrt[2], which is simply an abbreviation of 2^(-1/2); and 1/Sqrt[2] is stored as 2^(-1/2). Understanding these concepts and interpretations is crucial to the proper use of pattern matching on symbolic expressions. Practice with FullForm[] and Head[]. The following example shows the internal interpretation of a + b I; Note that I is really Complex[0, 1].

a + b I // FullForm

A common use of heads is to cause restrictions in function definitions. One reason to do this is to avoid unpleasant crashes if an object of the wrong type is sent to your function by mistake.

f[x_Complex] := 17{f[3], f[1 + I], f[I], f[a + (Sqrt[2]) I]}

The Resultant Algorithm

Every nonatomic radical expression is built up from simpler such expressions. So it is natural to use a recursive procedure to get the desired polynomial. Thus the problem reduces to several steps: the following two are the most difficult:

  If α is a root of p(x) and β is a root of q(x), find a polynomial having α + β as a root.
  If α is a root of p(x) and β is a root of q(x), find a polynomial having α·β as a root.

Of course there are similar steps involving exponents and division, but they are easier.

First, however, let's deal with some base cases. If α has the form a + b I where and b are rational, then it is a simple matter to get an explicit quadratic that α satisfies. But if b = 0, then of course we would prefer x – a as an answer, as opposed to a quadratic in x. The following two cases handle all this. Note the use of |, which means "or" when applied to patterns. And note the use of v_:x to set the default variable in the polynomial to x.

Polynomial[a_Integer | a_Rational, v_:  x] := v - aPolynomial[z_Complex, v_:  x] := (v^2) - 2 Re[z] v + Abs[z]^2{Polynomial[3], Polynomial[3 + 7 I]}

One problem is that a rational input will return an object in Q[x].

Polynomial[3 + ((7 I)/4)]

Of course, denominators can be cleared so that the result is in Z[x], which we prefer. We define a function, fix, to do this and then use fix liberally in our definition of Polynomial. We ask fix to determine the variable (via Variables) and then clear by a common denominator (via LCM). We also make fix clear any common divisor from the coefficients.

Clear[Polynomial]; fix[p_] := 
      Module[{temp, v = Variables[p]}, 
        temp = Together[p LCM@@Denominator/@CoefficientList[p, v]]; 
        Simplify[(temp/GCD)@@CoefficientList[temp, v]]]Polynomial[a_Integer | a_Rational, v_Symbol:  x] := fix[v - a]Polynomial[z_Complex, v_Symbol:  x] := 
      fix[(v^2) - 2 Re[z] v + Abs[z]^2]Polynomial[3 + ((7 I)/4)]

Before turning to sums and products, here is a case that will handle square roots for us. Recall from section 1 that 1/Sqrt[2] has the internal form 2^(-1/2). Also, observe that if αn is a root of p(x), then  αn/m is  is a root of p(xm). With these points in mind, it is easy to implement the rule we want; and since the recursion will always see an element of Z[x], we need not use fix.

\(Polynomial[e_\^Rational[n_, m_], v_:  x] := 
      Polynomial[e\^n, v] /. v -> v\^m\)

Polynomial[(1/(Sqrt[2]))]

Polynomial[(Sqrt[(3/5)])]

Now we turn to the sum question. Readers familiar with the elementary theory of algebraic fields — in particular, the theorem that the algebraic numbers form a field, and so is closed under addition — might recall a method based on linear independence of powers. Indeed, I once implemented that method, but it is very naive from a computational-algebra perspective. Some years ago Joe Buhler pointed out to me that the resultant of the two polynomials does the job, and is easy to compute. Since Resultant is included in Mathematica, that is the way to proceed.

Here is the precise idea. The resultant of two polynomials, p(x) and q(x), is defined to be
Π Π (βj - αi), where α and β denote the roots of p and q, respectively. It turns out that there is a fast way to compute the resultant, and it is built into Mathematica.

Resultant[(x^2) - 2, (x^2) - 5, x]

Now suppose that p and q are monic polynomials in Z[x]. Consider the function h(x), defined to be Resultant[q(t), p(x - t), t]. The third argument means that the first two are viewed as polynomials in t; thus the result is a function of x. To see that the roots of h(x) are precisely the set of sums of roots of p and q, observe that, for any x, the roots of p(x - t) are x - αi, and so h(x) = Π Π[βj - (x - αi), or Π Π(βj + αi - x). Now it is by no means obvious that h(x) is a polynomial in Z[x], but this is in fact true, and follows from a characterization of the resultant as a certain determinant. So this solves our sum problem. See [Pohst and Zassenhaus 1989] for more details.

A slight complication is that the sum case requires a call to Polynomial[a, v - t]. This means that the variable is to be v - t; but fix will not work in such a case, since it basically assumes its input is a polynomial in one variable. So we instead call Polynomial[a, t] and then replace t by  v - t.

Polynomial[e_ + f_, v_Symbol:  x] := 
      Module[{t}, 
        fix[Expand[
            Resultant[Polynomial[f, t], Polynomial[e, t] /. t -> v - t, 
              t]]]]

Polynomial[(Sqrt[2]) + (Sqrt[3])]

Polynomial[(Sqrt[2]) + (Sqrt[3]) + (Sqrt[5])]

Of course, it never hurts to check. Here we carry out a symbolic verification of the preceding result, but for larger problems, it is more expedient to check numerically, making sure to use an appropriately large number of significant digits (via N[***, largeNumber]).

Expand[% /. x -> (Sqrt[2]) + (Sqrt[3]) + (Sqrt[5])]

The case of products is similar to sums. The key result is that the products of the roots of p (degree n) and q (degree m) are the roots of p(0)m Resultant[tn p(x/t), q(t), t].

Polynomial[e_ f_, v_Symbol:  x] := 
      Module[{t, p = Polynomial[e, v], q}, q = Polynomial[f, t]; 
        fix[Expand[
            (((p /. v -> 0))^Exponent)[q, t] 
              Resultant[(t^Exponent)[p, v] p /. v -> (v/t), q, t]]]]

Polynomial[2 + (Sqrt[(3/5)]) ((((1/9)))^(1/3))]

Note that the algorithm does not necessarily produce an irreducible polynomial. In some cases, such as the one just computed, it is easy to determine the minimal polynomial for the given algebraic number: just factor the output.

Factor[%]

In more complicated cases one might want to choose from the factors. Since version 3.0 produces the minimal polynomial, we will not discuss this further.

Two more cases will complete our algorithm, in that it will handle any radical. First we must deal with positive integer powers; a simple use of Expand does the job.

\(Polynomial[e_\^n_Integer?Positive, v_:  x] := 
      fix[Polynomial[Expand[e\^n]]]\)

Polynomial[(1 + Sqrt[2])^2]

Finally there is the case of negative exponents. Properly implemented, this takes care of division as well.

\(Polynomial[e_\^n_Integer?Negative, v_:  x] := 
      Module[{q = Polynomial[e\^\(-n\), v]}, 
        fix[Expand[v\^Exponent[q, v]\ \((q /. v -> 1\/v)\)]]]\)

Polynomial[((Sqrt[3])/(1 + (1/(((1 + (Sqrt[2]) I))^3))))]

Polynomial[((Sqrt[1 + (Sqrt[17])])/((2^(1/3)) + I))]

Applications

Problem B4 on the 1995 Putnam Competition involved expressing (2207 + √22072 – 4)1/8 in the form a + b√2. Our routine can easily do that. We first find a polynomial for the eighth root.

e = (2207 + Sqrt[2207^2 - 4])/2;
Polynomial[e^(1/8)]

Now we find all the roots of the 16th degree polynomial (fortunately, this can be done!).

x /. Solve[%[[1]] == 0]

Only one of these is a positive real, so the answer must be (3 + √5)/2. It is easy to double-check.

{e, Simplify[Expand[((((1/2) (3 + (Sqrt[5]))))^8)]]}

Here's another radical problem; this one was presented in The Mathematica Journal (2:2, p. 39) as an application of Recognize.m. That calls LatticeReduce, a very sophisticated algorithm. The problem was to simplify Sqrt[(3 + Sqrt[(3 + Sqrt[5])/2])]. We can do it by finding a polynomial that the radical satisfies and then factoring it.

a = (Sqrt[3 + (Sqrt[(1/2) (3 + (Sqrt[5]))])]); Factor[Polynomial[a]]

Now we must determine which factor has a as a root. We try the first quartic first.

Roots[%[[1, 1]] == 0, x]

Then we check numerical values.

N[{a, Map[Last, %]}]

We can therefore conclude that the given radical simplifies to the third root in the radical list, Sqrt[(7+Sqrt[5])/2].

It is remarkable that so much can be accomplished with very few lines of code. I should add that the logic and the order of the cases is a little subtle, and the code presented here took some time to work out. Even though this functionality is built into version 3.0 (see next section), the techniques here should prove useful in other polynomial algebra computations. And, of course, it is nice to have some understanding of the internals of the powerful Mathematica functions.

4. The Future Is Now

In version 3.0 one can accomplish what we have done here very simply using the new RootReduce function

a = (Sqrt[3 + (Sqrt[(1/2) (3 + (Sqrt[5]))])]); RootReduce[a]

This output indicates that a is the fourth root (in order of size) of the polynomial represented by the pure function [Editor's note: see "Algebraic numbers in Mathematica 3.0," by Adam Strzebonski, this issue]. Moreover, it is guaranteed that the polynomial is the minimal polynomial for the radical. We can get a radical expression as follows.

ToRadicals[%]

Now we could, as in section 3, find the roots of this quartic and conclude that a is Sqrt[(7+Sqrt[5])/2].

However RootReduce does more in that, sometimes, it reduces the given radical immediately. Consider the Putnam quadratic of section 3.

a = (2207 + Sqrt[2207^2 - 4])/2;
RootReduce[a]

REFERENCES

Pohst, M. and Zassenhaus, H., Algorithmic Algebraic Number Theory, Cambridge Univ. Press, Cambridge, U.K., 1989

ABOUT THE AUTHOR

Stan Wagon uses Mathematica extensively in his research, teaching, and exposition and is especially appreciative of the new ways in which Mathematica allows us to look at mathematical objects. He is author or coauthor of several books, including Mathematica in Action, Animating Calculus, VisualDSolve, and Which Way Did the Bicycle Go?. He also teaches a high-altitude Mathematica course every summer in the mountains of Colorado.

Stan Wagon
Department of Mathematics
Macalester College
St. Paul, MN 55105
wagon@macalester.edu

ELECTRONIC SUBSCRIPTIONS

Included in the distribution for each electronic subscription in the file PolyRad.nb containing Mathematica code for the material described in this article.