Hi, a couple of days ago I released version 0.0.2 of PolynomialRoots.jl <https://github.com/giordano/PolynomialRoots.jl> (formerly known as CmplxRoots.jl, some people suggested to change the name to something more descriptive), a fast root finder of real and complex polynomials.
This is a Julia implementation of the algorithm described in - J. Skowron & A. Gould, 2012, "General Complex Polynomial Root Solver and Its Further Optimization for Binary Microlenses", arXiv:1203.103 <http://arxiv.org/abs/1203.1034> PolynomialRoots.jl is a registered Julia package, so you can install it with the package manager: Pkg.add("PolynomialRoots") Two functions are exposed to the users: roots(polynomial[, roots, polish=true]) roots5(polynomial[, roots]) The first function can be used to solve polynomials of any order, the second one is specialized and optimized for (and works only for) polynomials of fifth-order This new version of PolynomialRoots.jl features support for arbitrary precision calculations. This is very useful for calculating roots of polynomials of very large degree (some hundreds) that couldn't be solved using double precision calculations. Actually, this cures inaccurate results also for low order polynomials, like second-order ones, where catastrophic cancellation is a problem. For example, the actual roots of (see https://en.wikipedia.org/wiki/Loss_of_significance) 94906265.625*x^2 - 189812534*x + 94906268.375 are x_1 = 1.000000028975958 x_2 = 1.000000000000000 but when you try to calculate them in double-precision you get julia> r = roots([94906268.375, -189812534, 94906265.625]); julia> r[1] 1.0000000144879793 - 0.0im julia> r[2] 1.0000000144879788 + 0.0im If you are interested in double-precision accuracy, you can work around this problem by calculating the roots with higher precision and then transforming the result to double-precision. What you have to do is only to pass BigFloat numbers to roots function: julia> r = roots([BigFloat(94906268.375), BigFloat(-189812534), BigFloat( 94906265.625)]); julia> Complex128(r[1]) 1.0000000289759583 - 0.0im julia> Complex128(r[2]) 1.0 + 0.0im as expected. This package is licensed under Apache License 2.0 or GNU Lesser General Public License version 3 or any later version, as well as under a "customary scientific license", which implies that if this code was important in the scientific process or for the results of your scientific work, you are asked for the appropriate citation of the paper Skowron & Gould 2012 (http://arxiv.org/abs/1203.1034). Cheers, Mosè
