# How to find and express all polynomial roots with greatest precision?

I've already plowed through the documentation and a bunch of other Questions on this, but I've only found partial or tangentially related answers.

I need numerical results for all of the complex roots of high-order polynomials with real coefficients. But Solve[] give really bad results, for example:

z[n_, c_] := If[n > 1, z[n - 1, c]^2 + c, c];
expr = z[8, c] // Expand; (*Example polynomial*)
numericRoots = Solve[expr == 0, c] // N;
expr /. numericRoots // Abs // Max


The result should be zero, but the numerical solutions numericRoots are so far off that the result is:

(* 1.5218*10^34 *)


Because of the high orders, I need to MMa to return the most precise results possible and express them as precisely as possible. I thought this was the way:

r = SetPrecision[NSolve[1 + 2*c^2 + 3*c^3 + 3*c^4 + 3*c^5 + c^6 == 0, c, WorkingPrecision -> Infinity], Infinity]


But that doesn't actually solve the polynomial for some reason:

(* {{c -> Root[1 + 2 #1^2 + 3 #1^3 + 3 #1^4 + 3 #1^5 + #1^6 &, 1]}, {c -> Root[1 + 2 #1^2 + 3 #1^3 + 3 #1^4 + 3 #1^5 + #1^6 &, 2]}, {c -> Root[1 + 2 #1^2 + 3 #1^3 + 3 #1^4 + 3 #1^5 + #1^6 &, 3]}, {c -> Root[1 + 2 #1^2 + 3 #1^3 + 3 #1^4 + 3 #1^5 + #1^6 &, 4]}, {c -> Root[1 + 2 #1^2 + 3 #1^3 + 3 #1^4 + 3 #1^5 + #1^6 &, 5]}, {c -> Root[1 + 2 #1^2 + 3 #1^3 + 3 #1^4 + 3 #1^5 + #1^6 &, 6]}} *)


I can get numerical results by doing this:

c /. r // N


But I'm not sure that all those high-precision instructions I made in the NSolve[] command above actually made it into these numbers. This isn't a big issue for 6th order polynomials, but it becomes one for 1000th order or so.

So, why didn't my NSolve[] command solve what I wanted it to solve? And how can I get the highest-precision numerical roots I want?

The code in my first example (which is supposed to return 0) is a good test of whether or not a proposed solutions is working right.

Thanks!

• As I said, there is no highest-possible finite precision, except for system limitations (i.e. memory). For instance, M will calculate ten million digits of the first root in about 7-8 sec. on my machine; to get the next twenty million takes nearly 20 sec. and only 8MB. I think you have to tell M how many you want. (I used N[Root[1 + 2 #1^2 + 3 #1^3 + 3 #1^4 + 3 #1^5 + #1^6 &, 1], 2 10^7]; //AbsoluteTiming) Oct 31 '15 at 22:59
• (1) As @JohnDoty noted in response, the roots are numerically fine. Only the residuals are seemingly amiss. In general the roots can be fine to machine precision (as these are) and you might still get bad residuals from cancellation error, even when using a nicer form of the function. Nov 1 '15 at 2:31
• (2) The Root objects are returned by Solve, not by NSolve. Nov 1 '15 at 2:32
• (3) the Root objects in the indicated result are in no sense "unresolved". They encode numeric approximations, isolating intervals, and can be refined (via N[...,...]) to arbitrarily high precision, subject only to computer and bignum numerics limitations. Nov 1 '15 at 2:34
• (4) The remark "The code in my first example (which is supposed to return 0) is a good test of whether or not a proposed solutions is working right." is some combination of correct, incorrect, and seriously incorrect. One must account for scaling, along with cancellation and truncation error, when testing residuals for being zero. The general rule is that you should not compare approximate values to zero unless you know good and well what the result indicates. In this instance it indicates only that truncation and maybe cancellation error are rearing their heads.... Nov 1 '15 at 2:44

Your numericRoots are OK. The core of your problem is that you use Expand to put your polynomial into power series form. High order power series are numerically unstable. The form your z function yields by itself is much better:

z[8, c] /. numericRoots // Abs // Max

(* 9.87841*10^-11 *)


"I can see that what you're saying is correct, but it leaves me even more confused than before. z[8,c] and z[8,c]//Expand are exactly the same because the coefficients are integers."

But once you apply one of your numericRoots rules, the expressions are not exactly the same. Your rules substitute "machine numbers" for c. When you do mixed arithmetic with integers and machine numbers, you wind up with machine numbers. Machine arithmetic is sensitive to the order of operations. Mathematically equivalent expressions that are different in form will generally yield different results in machine arithmetic.

terms = expr /. Plus -> List

(* large expression that looks like an alluvial fan *)


You can see that while the coefficients are indeed integers, some are very large, too large for accurate representation as machine numbers. Then, even though the roots are modest, you're raising them to high powers, so the numbers you're adding are very large indeed:

terms /. numericRoots[[9]] // Abs // Max

(* 3.02339*10^49 *)


When you're adding a bunch of numbers of this magnitude in machine precision (~16 decimal digits), a residual of 10^34 is about as close to zero as you can reasonably expect.

"if simply changing the polynomial to a different but equivalent form creates a residual that humongous, how can I have any confidence in the root Mathematica gave me?"

Computation using machine numbers is fast, but (like other programming systems) Mathematica does not track the resulting precision. Sometimes, as in your formulation, the result contains no significant digits. There are, however, ways you can check.

One way is to use exact roots:

roots = Solve[expr == 0, c];


Don't look at them: they are mostly "algebraic numbers" expressed as complicated Root objects. But they are true, precise numbers. They are expressed in notation that is unfamiliar to most, but there is no general way to express algebraic numbers in familiar notation. Be happy that Mathematica can do this: in some situations it's very powerful. Mathematica is a bit lazier about computing with general algebraic numbers than with other kinds of numbers, but it can here if you coax it, and wait a few minutes:

Simplify[expr /. roots] // Abs // Max

(* 0 *)


So Mathematica asserts that its results are perfect ;-)

Another approach is to use arbitrary precision instead of machine precision:

numericRoots = N[Solve[expr == 0, c], 50];
expr /. numericRoots // Abs // Max // FullForm

(* 0-2.568868432106492 *)


In arbitrary precision arithmetic, Mathematica tracks the precision of the result. Here, the result is 0, but with an accuracy of -2.56... decimal digits, so it might be 300 or so. That's a lot better than 10^34, and Mathematica informs you of the accuracy. You can do much better than a precision of 50 for your roots, of course:

numericRoots = N[Solve[expr == 0, c], 1000];
expr /. numericRoots // Abs // Max // FullForm

0947.4311315678933


The residuals are zero to better than 947 decimal places.

• I can see that what you're saying is correct, but it leaves me even more confused than before. $z[8,c]$ and $z[8,c]//Expand$ are exactly the same because the coefficients are integers. So how on earth can one have a residual of order 10^-11 and the other 10^+34? And if simply changing the polynomial to a different but equivalent form creates a residual that humongous, how can I have any confidence in the root Mathematica gave me? And what if I only had the polynomial to work with, not the un-Expanded z[n_,c_] form? Sorry my response isn't more coherent, but your Answer kinda blew my mind. Nov 1 '15 at 5:44
• Also, it seems like something much deeper has to be wrong for me to get a residual of order 10^+34. That's 20 orders of magnitude larger than the coefficients of z[8,c], much to large to be caused by a small error in the root. Nov 1 '15 at 5:47
• @Jerry, people have gotten into unluckier situations than yours due to catastrophic cancellation in the expanded form. Your polynomial is... specially structured; if you don't mind a bit of a wait, I'll write a detailed answer later. For now: instead of Expand[], try using HornerForm[]. Nov 1 '15 at 6:06
• @Jerry, "z[8, c] /. numericRoots // Abs // Max returns the best results": it's not the expanded form with the risk of cancellation looming over it, so I'd expect that. But in that case, you definitely will want to wait for my post later, if Horner doesn't do the job... Nov 1 '15 at 6:42
• Just for fun, I took the decimal machine number representation of Mathematica's ninth root (the toughest, c->-1.9997740486937272) and fed it to Python versions of z[8,c] and expr. Python does a little better on the residuals here, but not a lot. Using z[8,c], the residual is -1.2048140263232199e-12, but for expr it's 7.03031399207685e+33. Nov 1 '15 at 23:17