Fix type stability of roots

[jwmerrill]

Currently, roots(p::Poly{Int64}) can return Array{Int64,1}, Array{Float64,1}, or Array{Complex{Float64},1} depending on the value of p (see examples for roots in README).

The convention in base julia is to have the return type be a function of the argument types only (not the argument values) wherever possible.

I think it would be better to return only real roots for real polynomials, and to always promote the results to floats.

If you want complex roots, you can promote your polynomial to the complex field before calling roots.

[simonster]

I think we'd want to always return a complex vector. For even very simple polynomials with only real roots, roots gives small complex parts:

julia> roots(Poly([1, 3, 3, 1]))
3-element Array{Complex{Float64},1}:
         -1.00001-0.0im
 -0.999995-8.53946e-6im
 -0.999995+8.53946e-6im

[jwmerrill]

This may be a fundamental problem in general, but I don't see why, in principle, it wouldn't be possible to get the roots correct to many more digits in this case.

[jwmerrill]

It looks like this issue is related to repeated roots. For an n-fold repeated root at x, the error in the computed root (using the current algorithm) is typically eps(x)^(1/n).

[vtjnash]

While it has been awhile since I coded this algorithm, Prof Edelman told me this was the best method for finding arbitrary roots, despite the perception of low precision results.

For better accuracy, perhaps we want a two stage approach. First adding a realroots function that returns all real roots (Wikipedia suggests this can be done most efficiently and with arbitrarily high precision using the Vincent–Collins–Akritas algorithm), then compute the remaining polynomial and use this formula to approximate the solution to the remaining complex roots. (with a fixed return type of Complex{Float64}).

Regardless, a realroots function would be quite nice to incorporate if someone wants to do some math research on the 'best' algorithm and make a pull request.

[vtjnash]

I now directly return the results of the eigenvalue calculation without trying to trim the imaginary part.