Build Status | Code Coverage |
---|---|
PolynomialRoots.jl
is a library for finding roots of complex univariate
polynomials, written in Julia.
This is an implementation in Julia of the General Complex Polynomial Root Solver, written in Fortran, by Jan Skowron and Andrew Gould. All the credits goes to them for the original algorithm.
The root finding algorithm employed in this library is described in
- J. Skowron & A. Gould, 2012, "General Complex Polynomial Root Solver and Its Further Optimization for Binary Microlenses", arXiv:1203.1034
This algorithm aims to be fast and precise, more than the well known zroots
procedure described in Numerical Recipes books, whose implementations in C and
Fortran are not available as free software, according to the
definition of the Free Software
Foundation.
PolynomialRoots.jl
can also take advantage of native arbitrary precision
capabilities of Julia and achieve more precise results.
Note: the adopted algorithm can give highly inaccurate results for polynomials of order above ~30. This can be mitigated by using multiple precision calculations (see example below).
PolynomialRoots.jl
is available for Julia 1.0 and later versions, and can be
installed with
Julia built-in package manager.
In a Julia session run the command
julia> Pkg.add("PolynomialRoots")
You may need to update your package list with Pkg.update()
in order to get the
latest version of PolynomialRoots.jl
.
Older versions are available also for Julia 0.4-0.7.
After installing the package, run
using PolynomialRoots
or put this command into your Julia script.
PolynomialRoots.jl
provides two functions to find the roots of complex
polynomials
roots(polynomial[, roots, polish=true, epsilon=1e-20])
roots5(polynomial[, roots, epsilon=1e-20])
roots
can be used to solve polynomials of any degree, roots5
is tailored to
(and works only for) polynomials of degree 5. This function exists because
one of the methods to calculate
gravitational microlensing
amplification by a binary lens requires solving a fifth-order complex
polynomial. roots5
should be more robust than roots
for this class of
polynomials.
The mandatory argument for both functions is:
polynomial
, the vector holding coefficients of the polynomial you want to solve, in ascending order, from the lowest to the highest. Coefficients can be complex and real
Optional arguments are:
roots
: if you roughly know in advance the position of the roots, you can pass the vector with the known roots to speed up convergence.roots
vector must be one element shorther thanpolynomial
. Inroots5
, the roots will be only polished. Elements of the vector can be complex and realpolish
(boolean keyword, only forroots
): if set totrue
, after all roots have been found by dividing original polynomial by each root found, all roots will be polished using full polynomial. Default isfalse
epsilon
(floating point optional keyword): this is used to determine the stopping criterion described in Skowron & Gould paper. If not set, it defaults to machine precision ofpolynomial
(androots
) argument(s). This is not the precision with which the roots will be calculated.
The functions return in output the Complex
vector with all roots of
polynomial
. Note: if roots
optional argument is provided, it is not
changed in-place.
Find the roots of
(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
This is a fifth-order polynomial, so we can find its zeros with both roots
and
roots5
:
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
PolynomialRoots.jl
handles polynomials with high-multiplicity roots as well.
For example, consider
(x + 1)^5 = x^5 + 5x^4 + 10x^3 + 10x^2 + 5x + 1
You can find its roots with
julia> roots([1, 5, 10, 10, 5, 1])
5-element Array{Complex{Float64},1}:
-1.0+0.0im
-1.0+0.0im
-1.0+0.0im
-1.0+0.0im
-1.0+0.0im
julia> roots5([1, 5, 10, 10, 5, 1])
5-element Array{Complex{Float64},1}:
-1.0+0.0im
-1.0+0.0im
-1.0+0.0im
-1.0+0.0im
-1.0+0.0im
Due to limited precision of Float64
type, extraction of roots of polynomials
can give inaccurate results, even for low-order polynomials. This is caused,
i.e., by
catastrophic cancellation
in computation of discriminant Δ = sqrt(b^2 - 4ac) of second-order polynomials.
For example, the actual roots
of
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. Julia natively supports multiple precision
calculations, so what you have to do is only to pass BigFloat
numbers to
roots
function (note: in order to correctly initialize an arbitrary
precision floating point constant it is better to use the big
string literal,
see http://docs.julialang.org/en/stable/stdlib/numbers/#Base.BigFloat):
julia> r = roots([big"94906268.375", -189812534, big"94906265.625"]);
julia> Float64(r[1])
1.0000000289759583
julia> Float64(r[2])
1.0
Note that in this case there is a trade-off between speed and higher accuracy and precision.
Another Julia package for finding roots of complex polynomials is
Polynomials.jl
, by Jameson Nash and
other contributors. This package does much more than finding roots of
polynomials (among other features, it can integrate and differentiate
polynomials). In order to solve the polynomial, Polynomials.jl
calculates
eigenvalues of its companion matrix, but PolynomialRoots.jl
is usually faster
by up to an order of magnitude and often slightly more precise. In addition,
Polynomials
cannot extract roots in arbitrary precision. If you are after
speed and precision, PolynomialRoots.jl
can still be a better option.
The package is developed at https://github.com/giordano/PolynomialRoots.jl. There you can submit bug reports, make suggestions, and propose pull requests.
The ChangeLog of the package is available in NEWS.md file in top directory.
The PolynomialRoots.jl
package is licensed under version 2.0 of the Apache
License or the GNU Lesser General Public License version 3 or any later version,
at your option. These are the same licenses used by the General Complex
Polynomial Root Solver.
A custom in the scientific comunity is (regardless of the licence you chose to use or distribute this software under) that if this code was important in the scientific process or for the results of your scientific work, we kindly ask you for the appropriate citation of the paper Skowron & Gould 2012 (http://arxiv.org/abs/1203.1034), and we would be greatful if you pass the information about the proper citation to anyone whom you redistribute this software to.
The authors of the General Complex Polynomial Root Solver, the original Fortran
library (http://www.astrouw.edu.pl/~jskowron/cmplx_roots_sg/) from which
PolynomialRoots.jl
has been translated, are Jan Skowron, Andrew Gould.
The original author of PolynomialRoots.jl
is Mosè Giordano.