Open stevengj opened 1 year ago
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.
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.
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:returns 6 nearly duplicate "roots"
which are not roots at all:
The correct roots are returned by
Polynomials.roots
, and they are completely different from the values returned above: