# Find real roots of a polynomial

Write a self-contained program which when given a polynomial and a bound will find all real roots of that polynomial to an absolute error not exceeding the bound.

### Constraints

I know that Mathematica and probably some other languages have a one-symbol solution, and that's boring, so you should stick to primitive operations (addition, subtraction, multiplication, division).

There is certain flexibility on the input and output formats. You may take input via stdin or command-line arguments in any reasonable format. You may allow floating point or require that some representation of rational numbers be used. You may take the bound or the reciprocal of the bound, and if you're using floating point you may assume that the bound will not be less than 2 ulp. The polynomial should be expressed as a list of monomial coefficients, but it may be big- or little-endian.

You must be able to justify why your program will always work (modulo numerical issues), although it's not necessary to supply full proofs inline.

The program must handle polynomials with repeated roots.

### Example

x^2 - 2 = 0 (error bound 0.01)


Input could be e.g.

-2 0 1 0.01
100 1 0 -2
1/100 ; x^2-2


Output could be e.g.

-1.41 1.42


but not

-1.40 1.40


as that has absolute errors of about 0.014...

### Test cases

Simple:

x^2 - 2 = 0 (error bound 0.01)

x^4 + 0.81 x^2 - 0.47 x + 0.06 (error bound 10^-6)


Multiple root:

x^4 - 8 x^3 + 18 x^2 - 27 (error bound 10^-6)


Wilkinson's polynomial:

x^20 - 210 x^19 + 20615 x^18 - 1256850 x^17 + 53327946 x^16 -1672280820 x^15 +
40171771630 x^14 - 756111184500 x^13 + 11310276995381 x^12 - 135585182899530 x^11 +
1307535010540395 x^10 - 10142299865511450 x^9 + 63030812099294896 x^8 -
311333643161390640 x^7 + 1206647803780373360 x^6 -3599979517947607200 x^5 +
8037811822645051776 x^4 - 12870931245150988800 x^3 + 13803759753640704000 x^2 -
8752948036761600000 x + 2432902008176640000  (error bound 2^-32)


NB This question was in the Sandbox for approximately 3 months. If you think it needed improving before posting, visit the Sandbox and comment on the other proposed questions before they are posted on Main.

• "You must be able to justify why your program will always work" Just in case someone needs help with that – Dr. belisarius May 16 '13 at 4:08
• @belisarius, ?? – Peter Taylor May 16 '13 at 7:20
• was intended as a joke :( – Dr. belisarius May 16 '13 at 11:57
• I know this is an old challenge, so don't feel obliged to answer if you don't feel like reopening it. (a) Can we write a function, or only a full program? (b) In case we can write a function, can we assume that the input uses some convenient data type, e.g., Python's fractions.Fraction (a rational type)? (c) Do we have to handle polynomials of degree <1? (d) Can we assume that the leading coefficient is 1? – Ell Dec 16 '15 at 9:30
• (e) With regard to polynomials with repeated roots, it's worth making a distinction between roots of odd and even multiplicities (the test cases only have roots of odd multiplicities.) While roots of odd multiplicity are not too hard to deal with, I'm not sure how tangible it is to correctly handle roots of even multiplicity numerically, especially since you only specify an error margin for the values of the roots, not for their existence. (...) – Ell Dec 16 '15 at 9:30

### Mathematica, 223

r[p_,t_]:=Module[{l},l=Exponent[p[x],x];Re@Select[NestWhile[Table[#[[i]]-p[#[[i]]]/Product[If[i!=j,#[[i]]-#[[j]],1],{j,l}],{i,l}]&,Table[(0.9+0.1*I)^i,{i,l}],2*Max[Table[Abs[#1[[i]]-#2[[i]]],{i,l}]]>t&,2],Abs[Im[#]]<t^.5&]]


This solution implements the Durand–Kerner method for solving polynomials. Note that this is not a complete solution (as will be shown below) because I cannot yet handle Wilkinson's Polynomial to the specified precision. First an explanation of what I'm doing: #[[i]]-p[#[[i]]]/Product[If[i!=j,#[[i]]-#[[j]],1],{j,l}]&: Thus function computes for each index i the next Durand-Kerner approximation. Then this line is encapsulated in a Table and applied using a NestWhile to the input points generated by Table[(0.9+0.1*I)^i,{i,l}]. The condition on the NestWhile is that the max change (over all terms) from one iteration to the next is greater than the specified precision. When all terms have changed less than this, the NestWhile ends and Re@Select removes the zeros that do not fall on the real line.

Example output:

> p[x_] := x^2 - 2
> r[p, .01]
{1.41421, -1.41421}

> p[x_] := x^4 - 8 x^3 + 18 x^2 - 27
> r[p, 10^-6]
{2.99999, 3., 3.00001, -1.}

> p[x_] := x^20 - 210 x^19 + ... + 2432902008176640000 (Wilkinson's)
> Sort[r[p, 2^-32]]
{1., 2., 3., 4., 5., 6., 7.00001, 7.99994, 9.00018, 10.0002, 11.0007, \
11.9809, 13.0043, 14.0227, 14.9878, 16.0158, 16.9959, 17.9992, \
19.0001, 20.}


As you can probably see, when the degree grows higher this method begins to bounce around the correct values, never really homing in completely. If I set the stopping condition of my code to be anything stricter than "from one iteration to next the guesses changed by no more than epsilon" then the algorithm never stops. I guess I should just use Durand-Kerner as input to Newton's method?

• Durand-Kerner also has potential problems with multiple roots. (Newton's method might not help much either - Wilkinson's polynomial is specifically chosen to be ill-conditioned). – Peter Taylor Jun 5 '13 at 13:45
• You are quite correct: I abandoned that course of action after zooming on on Wilkinson's near x=17, it's an absolute mess. I worry that I'll have to go for a symbolic solution with a Groebner basis to get much more accuracy. – Kaya Jun 5 '13 at 14:37