[GAP Forum] Issue with Float package

[Gordon Royle]

On 20 Apr 2015, at 12:41 am, Stephen Linton <steve.linton at st-andrews.ac.uk<mailto:steve.linton at 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 roots).

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

x^20-40*x^18+610*x^16-8*x^15-4640*x^14+40*x^13+ 19475*x^12-560*x^11-46380*x^10+5560*x^9+61610*x ^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

Gordon