Find real roots of a polynomial

Question

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.

Answer 1

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?