search

giordano / PolynomialRoots.jl

Fast complex polynomial root finder, with support for arbitrary precision calculations
Other
53 stars 15 forks source link

roots hangs forever with input #2

Closed dfagnan closed 7 years ago

dfagnan commented 7 years ago

@giordano Is roots guaranteed to converge? Should there be a maximum number of iterations or a timeout?

I have attached a jld file with coefficients that are likely numerically unfriendly. roots hangs unless epsilon is significantly larger than 1E-8 or so. Let me know if I can help debug, I would love to learn more about this awesome implementation!

badPolynomial.zip

giordano commented 7 years ago

Thanks for the report. I guess it should converge somewhere, never seen it hanging, even with polynomials of much higher degree. For large order polynomials it usually converges to the wrong solutions, but decreasing the epsilon might help.

This may be an implementation bug, it would be great if you could debug it, at least find where it hangs. I don't think I'll have the time to do it for the next few weeks.

I would love to learn more about this awesome implementation!

You can read the paper at https://arxiv.org/abs/1203.1034 ;-)

giordano commented 7 years ago

I tried the original Fortran routines, they seem to work with the coefficients you provided, so it may well be a bug in my implementation.

dfagnan commented 7 years ago

Ok, great thanks for letting me know, I will try to debug when I have time!

giordano commented 7 years ago

A hint: the solver gets stuck in laguerre2newton, where there is an infinite loop, maybe an exit condition is missing. Thank you!

giordano commented 7 years ago

There was an error with an iteration index in laguerre2newton function. Now the function gets the right solution with an error at most of 2e-11:

julia> poly = [63.22256105356723, 271.8182248226112, 144.3588342039991, -25.60790629850817,
               952.6388846106129, -32.65159275777219, 62.19331611327388, 44.70637211946786,
               -28.265078307398895, -37.63653902029289, -72.26102355751738, -25.501990478720046,
               -47.40236121905153, -71.92379520244637, -57.977452001749555];

julia> Base.Test.@test zeros(length(poly)-1) ≈ 
           PolynomialRoots.evalpoly(roots(poly), poly) atol = 2e-11
Test Passed

Please tell me if you need a tagged version.

Thank you so much for the bug report (and the coefficients!).

dfagnan commented 7 years ago

Awesome turnaround, thanks!