sometimes returns bogus duplicate roots

[stevengj]

As noted on discourse, the roots function seems to return bogus roots, which are all nearly identical, for some polynomials. It happens pretty often for me if you just randomly try different polynomials of degree 6. Here is one example:

c = [0.1767237967771577 - 0.8611076733845967im, 0.00017422934461020398 - 9.453901157858016e-6im, 2.8796166645520005e-5 - 6.235942236501559e-5im, 0.00011096228622067326 + 6.305229743814886e-5im, -4.383872407211953e-5 - 3.677479055394419e-5im, -5.464453142543401e-6 + 6.577470932911938e-5im, 1.0 + 0.0im]
r = roots(c, polish=true)

returns 6 nearly duplicate "roots"

6-element Vector{ComplexF64}:
 1.0126353489009172 + 7.1699666737369805im
  1.012635472445031 + 7.169966573453136im
 1.0126355118156847 + 7.169966616995025im
 1.0126353647489426 + 7.169966609253118im
 1.0126354389431256 + 7.169966613009901im
  1.012635426566087 + 7.169966602733274im

which are not roots at all:

julia> evalpoly.(r, Ref(c))
6-element Vector{ComplexF64}:
 -96023.81180290483 + 107521.12470406065im
  -96023.7931286698 + 107521.12824563321im
 -96023.79406644785 + 107521.13519457991im
 -96023.80469749296 + 107521.12117898901im
 -96023.79932061117 + 107521.12823827742im
 -96023.79933709523 + 107521.12631673364im

The correct roots are returned by Polynomials.roots, and they are completely different from the values returned above:

julia> import Polynomials

julia> r2 = Polynomials.roots(Polynomials.Polynomial(c))
6-element Vector{ComplexF64}:
   -0.936358238274595 - 0.28505582219111214im
   -0.714980912681819 + 0.6683037772456318im
 -0.22129525444596665 - 0.9534008044046407im
  0.22128062807035184 + 0.9533676310467337im
   0.7150642724349877 - 0.6683776196201094im
   0.9362949693501847 + 0.2850970632141694im

julia> evalpoly.(r2, Ref(c))
6-element Vector{ComplexF64}:
   8.881784197001252e-16 + 1.4432899320127035e-15im
   5.218048215738236e-15 - 3.3306690738754696e-15im
 -2.7755575615628914e-17 - 4.9960036108132044e-15im
  -2.248201624865942e-15 + 2.9976021664879227e-15im
  1.7763568394002505e-15 - 2.55351295663786e-15im
   4.440892098500626e-16 + 2.55351295663786e-15im

[giordano]

Trying the original library (code of the roots! and root functions based on https://github.com/giordano/PolynomialRoots.jl/blob/194c4ffa4d1e94e72b734b971283100ebb093b2b/src/CmplxRoots.jl):

import Downloads

src = Downloads.download("http://www.astrouw.edu.pl/~jskowron/cmplx_roots_sg/cmplx_roots_sg.f90")
libcmplxroots = joinpath(tempdir(), "libcmplxroots.so")
run(`gfortran -x f95 $(src) -shared -fPIC -o $(libcmplxroots)`)

c = [
    0.1767237967771577 - 0.8611076733845967im,
    0.00017422934461020398 - 9.453901157858016e-6im,
    2.8796166645520005e-5 - 6.235942236501559e-5im,
    0.00011096228622067326 + 6.305229743814886e-5im,
    -4.383872407211953e-5 - 3.677479055394419e-5im,
    -5.464453142543401e-6 + 6.577470932911938e-5im,
    1.0 + 0.0im,
]

function roots!(roots::Vector{ComplexF64}, poly::Vector{ComplexF64},
                degree::Integer, polish::Bool, use_roots::Bool)
    ccall((:cmplx_roots_gen_, libcmplxroots), Ptr{Cvoid},
          (Ptr{ComplexF64}, # roots
           Ptr{ComplexF64}, # poly
           Ptr{Cint}, # degree
           Ptr{Cint}, # polish_roots_after
           Ptr{Cint}),# use_roots_as_starting_points
          roots, poly, Ref{Cint}(degree),
          Ref{Cint}(polish), Ref{Cint}(use_roots))
    return roots
end

function roots(poly::Vector{<:Number}; polish::Bool=false)
    degree = length(poly) - 1
    roots!(Array{ComplexF64}(undef, degree), float(complex(poly)), degree, polish, false)
end

r = roots(c, polish=true)

gives

julia> evalpoly.(r, (c,))
6-element Vector{ComplexF64}:
   7.216449660063518e-16 + 5.551115123125783e-16im
                     0.0 + 4.440892098500626e-16im
 -1.1102230246251565e-16 - 2.220446049250313e-16im
  -4.440892098500626e-16 - 1.1102230246251565e-16im
                     0.0 + 3.3306690738754696e-16im
  -2.498001805406602e-16 + 1.1102230246251565e-16im

I guess there are small discrepancies in the implementation with the upstream library.

[stevengj]

The good news is hopefully that you should be able to trace through the execution in Julia vs Fortran for this particular example and see where they start to disagree in order to find the bug, if the algorithms are supposed to be identical.