On 20 Apr 2015, at 12:41 am, Stephen Linton 
<steve.lin...@st-andrews.ac.uk<mailto:steve.lin...@st-andrews.ac.uk>> wrote:

I do observe though, that the imaginary parts of the roots are very small — 
less that 10^-7, so I imagine that the “correct” roots are very close together 
and rounding error is then shifting them off the real line.

Steve is right, and the given roots are correct up to the limitations of fixed 
precision floating point arithmetic, so this does not imply an error in the 
Float package.

Ultimately, the problem is that fixed precision floating point arithmetic is 
not suitable to compute roots except for low-degree polynomials with small 
coefficients. This is because finding the complex roots of polynomials, 
especially if there are multiple roots, is quite sensitive to rounding errors. 
Although polynomial roots vary continuously with the coefficients, it is easy 
to construct examples with only modest sized coefficients for which normal 
floating point arithmetic will produce totally incorrect answers.

Increasing the precision will improve things, of course, but it is not really a 
satisfactory approach, because you’re only specifying the precision of each 
individual operation, while what you want to know is the precision of the final 
results (i.e the roots themselves).

Fortunately, there are freely-avaiable program that find roots to guaranteed 
OUTPUT precision (barring bugs) by bounding the total accumulated error during 
the computation, which is a much more robust approach. Two such programs that I 
have used with polynomials of degree as high as 1000 are the function 
“polroots” in Pari/GP, and the stand alone program “unisolve” in the package 
MPSolve. Both of these will work much faster if they are told in advance that 
the polynomials have only real roots (or that you are only interested in real 

There are also a few things you can do before even sending it to the root 
finder - factorise the polynomial, and removing repeated roots will all help in 
reducing the size of the coefficients.

For your particular polynomial, we see that

^8-19600*x^7-42220*x^6+25784*x^5+8265*x^4-11560* x^3+3700*x^2-400*x

factorises (over the integers) to give us

(-4 + x) (-2 + x)^2 x (1 - 3 x + x^2)^2 (-1 + x + x^2)^4 (5 + 5 x + x^2)^2

from which the roots are easily deduced by hand.

Hope this helps

Forum mailing list

Reply via email to