-
-
Notifications
You must be signed in to change notification settings - Fork 15
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Wrong results #19
Comments
I think this is going haywire because of the magnitude difference: julia> using PolynomialRoots
julia> d = 20; c = [randn(ComplexF64); randn(ComplexF64, d - 1) * 1e-4; 1];
julia> PolynomialRoots.evalpoly(roots(c), c)
20-element Array{Complex{Float64},1}:
9.081074789250304e16 + 1.2232077690076664e17im
9.081074789250294e16 + 1.2232077690076717e17im
9.081074789250203e16 + 1.2232077690076765e17im
9.081074789250253e16 + 1.2232077690076682e17im
9.081074789250288e16 + 1.2232077690076715e17im
9.081074789250278e16 + 1.2232077690076686e17im
9.081074789250277e16 + 1.223207769007676e17im
9.081074789250288e16 + 1.2232077690076754e17im
9.081074789250253e16 + 1.2232077690076669e17im
9.081074789250283e16 + 1.2232077690076715e17im
9.081074789250288e16 + 1.2232077690076699e17im
9.081074789250286e16 + 1.223207769007674e17im
9.081074789250362e16 + 1.2232077690076762e17im
9.081074789250314e16 + 1.2232077690076688e17im
9.081074789250256e16 + 1.2232077690076763e17im
9.081074789250288e16 + 1.2232077690076726e17im
9.081074789250272e16 + 1.2232077690076722e17im
9.081074789250307e16 + 1.2232077690076648e17im
9.08107478925015e16 + 1.2232077690076675e17im
9.081074789250309e16 + 1.2232077690076734e17im It's interesting that all the points are being evaluated more or less at the same, wildly wrong, value 🤔 When I don't scale the central elements, I get more reasonable results: julia> d = 20; c = [randn(ComplexF64); randn(ComplexF64, d - 1); 1];
julia> PolynomialRoots.evalpoly(roots(c), c)
20-element Array{Complex{Float64},1}:
8.326672684688674e-16 - 7.549516567451064e-15im
8.881784197001252e-15 + 4.9960036108132044e-15im
3.3306690738754696e-15 + 7.105427357601002e-15im
1.8189894035458565e-11 + 0.0im
9.103828801926284e-15 + 3.552713678800501e-15im
4.440892098500626e-16 + 3.6637359812630166e-15im
1.0658141036401503e-14 + 0.0im
1.4210854715202004e-14 + 0.0im
1.3655743202889425e-14 - 5.495603971894525e-15im
-1.5987211554602254e-14 - 1.687538997430238e-14im
-9.947598300641403e-14 + 5.684341886080802e-14im
7.66053886991358e-15 - 1.1102230246251565e-15im
1.3322676295501878e-15 + 1.2434497875801753e-14im
7.105427357601002e-15 - 3.552713678800501e-15im
1.2212453270876722e-14 - 1.7763568394002505e-15im
0.0 - 9.992007221626409e-15im
5.329070518200751e-15 - 3.552713678800501e-15im
-8.881784197001252e-16 + 5.329070518200751e-15im
2.6645352591003757e-15 - 3.552713678800501e-15im
4.440892098500626e-16 + 1.942890293094024e-16im I need to check whether the original Fortran implementation has the same troubles. |
Scaling the central values is exactly the point of the example :) |
Sure, I just wanted to make sure the root finder didn't get broken recently 🙂 |
This has something to do with how the machine epsilon is used:
The idea was to vary the number of bits in Turns out that precision is actually relatively good at around 20 bits of precision, but somewhere between 20 and 30 bits (depending on the input value), the results suddenly become completely wrong. Also note that the maximum error stabilizes, at some point it stops depending on the I think the machine epsilon is basically the most important thing that changes between, e.g., 22 and 23 bits of precision, so the bug is probably related to |
Overriding |
Following similar approach as #25 (comment): import Downloads, PolynomialRoots
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)`)
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 julia> d = 20; c = [randn(ComplexF64); randn(ComplexF64, d - 1) * 1e-4; 1];
julia> evalpoly.(roots(c), (c,)) # With upstream library
20-element Vector{ComplexF64}:
6.795511325941872e34 - 1.5437795762995263e35im
6.795511325941971e34 - 1.543779576299517e35im
6.795511325941929e34 - 1.5437795762995083e35im
6.79551132594187e34 - 1.5437795762995149e35im
6.795511325941931e34 - 1.5437795762995116e35im
6.795511325941875e34 - 1.5437795762995158e35im
6.795511325941925e34 - 1.5437795762995188e35im
6.795511325941927e34 - 1.5437795762995103e35im
6.795511325941969e34 - 1.5437795762995164e35im
6.795511325941953e34 - 1.5437795762995027e35im
6.79551132594192e34 - 1.543779576299509e35im
6.795511325941883e34 - 1.54377957629951e35im
6.795511325941962e34 - 1.5437795762995177e35im
6.795511325941907e34 - 1.543779576299515e35im
6.795511325941883e34 - 1.5437795762995092e35im
6.795511325941816e34 - 1.5437795762995114e35im
6.795511325941931e34 - 1.5437795762995142e35im
6.795511325941916e34 - 1.5437795762995162e35im
6.795511325941931e34 - 1.5437795762995142e35im
6.795511325942111e34 - 1.5437795762995188e35im
julia> evalpoly.(PolynomialRoots.roots(c), (c,)) # With PolynomialRoots.jl
20-element Vector{ComplexF64}:
3.423287799691176e15 - 8.646704914995811e15im
3.423287799691144e15 - 8.646704914995869e15im
3.423287799691134e15 - 8.646704914995846e15im
3.423287799691135e15 - 8.646704914995831e15im
3.423287799691126e15 - 8.64670491499579e15im
3.423287799691176e15 - 8.646704914995856e15im
3.4232877996911025e15 - 8.646704914995833e15im
3.423287799691136e15 - 8.646704914995835e15im
3.423287799691091e15 - 8.646704914995843e15im
3.423287799691131e15 - 8.64670491499586e15im
3.4232877996911195e15 - 8.646704914995844e15im
3.423287799691112e15 - 8.646704914995826e15im
3.423287799691128e15 - 8.646704914995828e15im
3.423287799691126e15 - 8.646704914995835e15im
3.423287799691089e15 - 8.646704914995832e15im
3.423287799691184e15 - 8.646704914995863e15im
3.423287799691118e15 - 8.646704914995827e15im
3.4232877996911175e15 - 8.64670491499583e15im
3.423287799691123e15 - 8.646704914995827e15im
3.423287799691152e15 - 8.646704914995844e15im They both get it wildly wrong (but results are different, perhaps because of #25. |
@nsajko perhaps looking at #25 might be easier since that one has somewhat less pathological values, and if we're lucky solving that issue might solve the issue here with BigFloat 🙂 🤞 If that's of any help, original Fortran source code, which is giving the correct result for #25, is at https://github.com/giordano/cmplx-roots-sg/blob/472b4d25a42e32a30a6817b937881de8516aa2f2/cmplx_roots_sg.f90, I tried to keep the Julia code as faithful to the original one as possible (but there are some limited changes), also comments should be the same, to make it easier to follow the code. |
I just briefly looked at the Fortran source code but it looks like polynomials of degree larger than 4 get deflated (i.e. you find one zero, and divide the polynomial by |
I did some test computations and found the following troubling result:
I am not sure whether this is a problem with the algorithm or the specific implementation.
The text was updated successfully, but these errors were encountered: