Find real roots of a polynomial
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}]&amp;,Table[(0.9+0.1*I)^i,{i,l}],2*Max[Table[Abs[#1[[i]]-#2[[i]]],{i,l}]]&gt;t&amp;,2],Abs[Im[#]]&lt;t^.5&amp;]]

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}]&amp;:
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:
&gt; p[x_] := x^2 - 2
&gt; r[p, .01]
{1.41421, -1.41421}

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

&gt; p[x_] := x^20 - 210 x^19 + ... + 2432902008176640000 (Wilkinson's)
&gt; 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 &quot;from one iteration to next the guesses changed by no more than epsilon&quot; then the algorithm never stops. I guess I should just use Durand-Kerner as input to Newton's method?