This looks interesting. It is in my opinion a bit unfortunate that you chose to use the GPL license for this since it means that it will not be usable by most packages in Julia which are under MIT license.
On Wednesday, April 27, 2016 at 1:15:35 PM UTC+2, 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è >
