How does this compare to AMVW.jl? https://github.com/andreasnoack/AMVW.jl
On Wednesday, April 27, 2016 at 9:15:35 PM UTC+10, Mosè Giordano wrote: > > Dear all, > > I am happy to announce the first version of CmplxRoots.jl > <https://github.com/giordano/CmplxRoots.jl>, 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.1034 > <http://arxiv.org/abs/1203.1034> > > I first wrote the package as a wrapper around the Fortran implementation > available at http://www.astrouw.edu.pl/~jskowron/cmplx_roots_sg/ but then > I rewrote it entirely in Julia, finding that it's even faster than > ccall'ing functions in the shared object. > > CmplxRoots.jl is a registered Julia package, so you can install it with > the package manager: > > Pkg.add("CmplxRoots") > > 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, but I found it to be a bit slower than roots() > . > > The only mandatory argument is the list of coefficients of the polynomial, > in ascending order. So, if you want to find all the roots of polynomial > > (x - i)*(x - 2)*(x - 3*i)*(x - 4)*(x - 5*i) = > x^5 - (9*i + 6)*x^4 + (54*i - 15)*x^3 + (138 - 57*i)*x^2 - (184 + 90*i)*x + > 120*i > > > you can use > > julia> roots([120im, -(184 + 90im), (138 - 57im), (54im - 15), -(6 + 9im), > 1]) > 5-element Array{Complex{Float64},1}: > -1.55431e-15+5.0im > 4.0+1.77636e-16im > 1.55028e-15+3.0im > -1.04196e-16+1.0im > 2.0-2.00595e-16im > > julia> roots5([120im, -(184 + 90im), (138 - 57im), (54im - 15), -(6 + 9im > ), 1]) > 5-element Array{Complex{Float64},1}: > 5.92119e-16+5.0im > 4.0-1.4803e-16im > 2.0+1.88202e-16im > -1.88738e-15+3.0im > 2.10942e-15+1.0im > > Optional roots argument is the array of starting values that can be used > to find the final roots (it won't be changed in-place). Note that even if > the roots are all real, the output array is always of type Complex{Float64}. > > Also Polynomials.jl package has a root finding function, but according to > some random tests performed with @benchmark (from Benchmarks.jl package) > it's roughly 20 times slower than roots() function provided by > CmplxRoots.jl. Of course do not expect results with machine precision (you > can have catastrophic cancellation even for simple quadratic polynomials), > but some random tests indicate that CmplxRoots.jl is slightly more precise > than Polynomials.jl (just compute abs2() of polynomials at the root). > > The package CmplxRoots.jl is licensed under the 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). > > Happy root-finding, > Mosè >
