Next: `RootOf' versus Radicals. Up: AM 511b Notes Previous: The GCD

# Roots of Polynomials and `RootOf'.

The fundamental theorem of algebra says that a polynomial of degreee n with complex coefficients (which includes real coefficients as a special case, of course) has exactly n complex roots, counting multiplicities. Sometimes one wants to know just what these roots are. More often, one just wants to know something about the location of the roots in the complex plane (are they all in , for example).

Many people only think they want roots of polynomials. For example, to solve an eigenvalue problem they first compute the characteristic polynomial, (which is a bad idea for computing eigenvalues but is the first method taught --- often the only method taught) and then ask the hard computational question `What are the roots of this polynomial?', never dreaming that the computational question `What are the eigenvalues of this matrix?' is far easier. Similarly using polynomial roots to specify equilibria of differential equations turns an easy computational problem into a hard one.

Nevertheless, sometimes you do want polynomial roots, and, after all, the computations can be interesting. The numerical approximation of polynomial roots is pursued farther in Part III, but for now we note that the routine `realroot` provides guaranteed bounds of arbitrary tightness for real roots of polynomials. For complex roots we will do something else later. The routines `sturm` and `sturmseq` let us count the number of roots in any given real interval, also.

```> readlib(realroot):
> with(linalg):
Warning: new definition for   norm
Warning: new definition for   trace
> A := hilbert(5):
> p := charpoly(A,x);   # This polynomial has real roots.
5   563  4    735781  3     852401   2      61501
p := x  - --- x  + ------- x  - --------- x  + ----------- x - 1/266716800000
315      2116800      222264000      53343360000

> st := time():
> realroot(p,Float(1,-10));  # This command gets the real roots to 10 places.
28243       56487         656911      5255289
[[----------, -----------], [----------, -----------],
8589934592  17179869184    2147483648  17179869184

195979213    97989607      895647649   3582590597
[-----------, ----------], [----------, -----------],
17179869184  8589934592    4294967296  17179869184

26921725877  13460862939
[-----------, -----------]]
17179869184   8589934592

> time_taken := time() - st;
time_taken := 7.000

> evalf("");

-5                -5
[[.3287917934*10  , .3287976142*10  ], [.0003058980219, .0003058980801],

[.01140749158, .01140749163], [.2085342186, .2085342187],

[1.567050691, 1.567050691]]
```
The routine `realroot` returns a list of pairs of dyadic rationals (rationals whose denominators are pure powers of two --- computing with these is more efficient than arbitrary rationals). Each pair of rationals represents an interval on the real line containing at least one root of the polynomial. In the above example there were five separate intervals returned, each of width less than , so we have successfully located all roots to ten decimal places absolutely.

Hereafter we regard the problem of the computation of arbitrarily accurate approximations or guaranteed bounds for roots of polynomials as a `kitchen problem' --- that is to say, a solved problem, one that we can carry out at any time we like. This isn't quite true in Maple --- there are high-degree polynomials that give trouble and we will see these later --- but for the purpose of the following discussion assume that this is so.

If we can compute roots at any time, then in particular we can compute them `later' and work with a symbolic representation of the roots. The way mathematicians do this is to use a phrase something like `Let be a root of the polynomial . Then ...' and compute formulae in terms of . This is possible in Maple as well, and often very useful. The syntax is similar to the above mathematical phrase, as is the meaning. Maple uses the `RootOf` construct, as follows.

```> z := RootOf(x^5+x+1,x);

5
z := RootOf(_Z  + _Z + 1)
```
z is any particular root (real or complex) of .
```> z^5;
5          5
RootOf(_Z  + _Z + 1)

> simplify(");
5
- RootOf(_Z  + _Z + 1) - 1
```
Use the definition to simplify expressions containing the RootOf.
```> simplify(1/z);
5          4
- RootOf(_Z  + _Z + 1)  - 1
```
Those expressions are ugly to look at, so we can use `alias` to clean up the presentation.
```> alias(alpha=RootOf(x^5+x+1,x));
I, alpha

> z;
alpha

> simplify(z^5);
- alpha - 1

> simplify(1/z);
4
- alpha  - 1

> simplify((z^3+z)/(1-z-z^5));
3
1/2 alpha  + 1/2 alpha

> simplify(alpha^3/(1-alpha));
4            3            2
1/3 alpha  + 1/3 alpha  - 2/3 alpha  - 2/3 alpha - 1/3
```
The Maple `RootOf` is best thought of as one particular root of a polynomial, but one which is not yet specified. It is not a set of roots, because arithmetic with sets is counterintuitive --- for example adding the set of roots with its negation does not give zero, but rather . In Maple, adding `RootOf(p,x)` to `-RootOf(p,x)` gives zero.

This means that Maple assumes that every single occurrence of a given `RootOf` refers to the same unspecified root. Usually this is the case, but sometimes it isn't and there are bugs due to this interpretation known in Maple. As of this writing there is no satisfactory solution to this bug --- we need an acceptable way of telling the difference between different unspecified roots.

We have already seen some examples of the uses for `RootOf`. We saw it used to represent an algebraic function for the purpose of integration. In fact, `RootOf` can be used to represent algebraic functions more generally than radical expressions can. It can also be used to give concise formulae for eigenvectors in terms of eigenvalues, as follows.

```> with(linalg):
Warning: new definition for   norm
Warning: new definition for   trace
> A := hilbert(3):
> p := charpoly(A,x);
3    23   2   127
p := x  - ---- x  + --- x - 1/2160
15       720

> alias(lambda=RootOf(p,x));
I, lambda

> eye := n -> array(1..n,1..n,identity);
eye := n -> array(1 .. n, 1 .. n, identity)

> B := evalm(A - lambda*eye(3));
[ 1 - lambda       1/2           1/3     ]
[                                        ]
B := [     1/2     1/3 - lambda       1/4     ]
[                                        ]
[     1/3          1/4      1/5 - lambda ]

> z := vector(3,0);
z := [ 0, 0, 0 ]

> eigenvector_corresponding_to_lambda := linsolve(B,z);
eigenvector_corresponding_to_lambda := [ _t[1],

154                                            2
- --- _t[1] + 42 _t[1] lambda - 80/3 _t[1] lambda ,
27

2
50/9 _t[1] - 60 _t[1] lambda + 40 _t[1] lambda  ]

>
>
>
> subs(_t[1]=1,");

154                          2                              2
[ 1, - --- + 42 lambda - 80/3 lambda , 50/9 - 60 lambda + 40 lambda  ]
27
```
In the example above, we have found all the eigenvectors of the matrix A in one formula --- once we choose this formula automatically gives the associated eigenvector.

The following example shows that algebraic functions and perturbation expansions can also be investigated using `RootOf`.

```> p := x^5 + e*x - 1;

5
p := x  + e x - 1

> alias(z=RootOf(p,x));
I, z

> alias(z0 = RootOf(x^5-1,x));
I, z, z0

> diff(z,e);
3                 4                    2                    3
e                 e  z                 z                  e z
- 320 ------------- + 64 ------------- - 625 ------------- - 500 -------------
5                  5                   5                   5
256 e  + 3125      256 e  + 3125       256 e  + 3125       256 e  + 3125

2  4
e  z
- 400 -------------
5
256 e  + 3125

> series(z,e);
2            3  2           4  3     21      5      6
z0 - 1/5 z0  e - 1/25 z0  e  - 1/125 z0  e  + ----- z0 e  + O(e )
15625

> convert(",polynom);
2            3  2           4  3     21      5
z0 - 1/5 z0  e - 1/25 z0  e  - 1/125 z0  e  + ----- z0 e
15625

> subs(e=1/100,");
156250000000021              2              3                 4
--------------- z0 - 1/500 z0  - 1/250000 z0  - 1/125000000 z0
156250000000000
```
This gives approximate formulae for the roots of . These formulae are valid no matter which root we take.
```>
> evalf(subs(z0=1,"));

.9979959920
```

Exercises

1. Use `realroot` to compute tight bounds on the real roots of and compare with the series solution given above.
2. Try to integrate z with respect to e, where .

Next: `RootOf' versus Radicals. Up: AM 511b Notes Previous: The GCD

Robert Corless
Wed Feb 28 18:25:56 EST 1996