Methods to find zeros (roots) of polynomials over given domains
This package provides the method poly_roots
to find roots of
univariate polynomial functions over the complex numbers, the real
numbers, the rationals, the integers, or Z_p. (A "root" is the name
for a "zero" of a polynomial function.) The package takes advantage of
other root-finding packages for polynomials within Julia (e.g.,
PolynomialRoots
for numeric solutions over the complex numbers and
PolynomialFactors
for exact solutions over the rationals and integers).
The basic interface is
poly_roots(f, domain)
The polynomial, f
, is specified through a function, a vector of
coefficients ([p0, p1, ..., pn]
), or as a Poly{T}
object, from the
the Polynomials.jl
package. The domain is specified by Over.C
(the
default), Over.R
, Over.Q
, Over.Z
, or over.Zp{p}
, with variants
for specifying an underlying type.
Examples:
julia> poly_roots(x -> x^4 - 1, Over.C) # uses `roots` from `PolynomialRoots.jl`
4-element Array{Complex{Float64},1}:
0.0+1.0im
1.0-0.0im
0.0-1.0im
-1.0+0.0im
julia> poly_roots(x -> x^4 - 1, Over.R)
2-element Array{Float64,1}:
1.0
-1.0
julia> poly_roots(x -> x^4 - 1, Over.Q) # uses `PolynomialFactors.jl`
2-element Array{Rational{Int64},1}:
-1//1
1//1
julia> poly_roots(x -> x^4 - 1, Over.Z) # uses `PolynomialFactors.jl`
2-element Array{Int64,1}:
-1
1
julia> poly_roots(x -> x^4 - 1, Over.Zp{5}) # uses `PolynomialFactors.jl`
4-element Array{Int64,1}:
4
1
3
2
Domains can also have their underlying types specified. For example, to solve
over the BigFloat
type, we have:
poly_roots(x -> x^4 - 1, Over.CC{BigFloat}) # `CC{BigFloat}` not just `C`
using DoubleFloats # significantly faster than BigFloat
poly_roots(x -> x^4 - 1, Over.CC{Double64})
4-element Array{Complex{DoubleFloat{Float64}},1}:
1.0 + 0.0im
-1.0 + 0.0im
-7.851872429582108e-35 - 1.0im
-7.851872429582108e-35 + 1.0im
There are other methods for Over.C
. This will use the AMVW method:
julia> poly_roots(x -> x^4 - 1, Over.C, method=:amvw) # might differ slightly
4-element Array{Complex{Float64},1}:
-2.1603591655723396e-16 - 0.9999999999999999im
-2.1603591655723396e-16 + 0.9999999999999999im
1.0 + 0.0im
-1.0 + 0.0im
This method is useful for high-degree polynomials (cf. FastPolynomialRoots):
n = 5000
rs = poly_roots(randn(n+1))
sum(isreal, rs) # 0
sum(!isnan, rs) # 1 (should be n)
rs = poly_roots(randn(n+1), method=:amvw)
sum(isreal, rs) # 6 ~ 6.04... = 2/π*log(n) + 0.6257358072 + 2/(n*π)
sum(!isnan, rs) # 5000 = n
This package uses:
The PolynomialRoots
package to find roots over the complex
numbers. The Roots
package can also be used. As well, an
implementation of the
AMVW algorithm can
be used. The default seems to be faster and as accurate as the others,
but for very high degree polynomials, the :amvw
method should be
used, as it will be faster and more reliable.
The PolynomialFactors
package to return roots over the
rationals, integers, and integers modulo a prime.
As well, it provides an algorithm to find the real roots of polynomials.
The main motivation for this package was to move the polynomial
specific code out of the Roots
package. This makes the Roots
package have fewer dependencies and a more focused task. In addition,
the polynomial specific code could use some better implementations of
the underlying algorithms.
In the process of doing this, making a common interface to the other root-finding packages seemed to make sense.
The package also provides
PolynomialZeros.AGCD.agcd
for computing an approximate GCD of
polynomials p
and q
over Poly{Float64}
. (This is used to
reduce a polynomial over the reals to a square-free
polynomial. Square-free polynomials are needed for the
algorithm used. This algorithm can become unreliable for degree 15
or more polynomials.)
PolynomialZeros.MultRoot.multroot
for finding roots of p
in
Poly{Float64}
over Complex{Float64}
which has some advantage if
p
has high multiplicities. The roots
function from the
Polynomials
package will find all the roots of a polynomial. Its
performance degrades when the polynomial has high
multiplicities. The multroot
function is provided to handle this
case a bit better. The function follows algorithms due to Zeng,
"Computing multiple roots of inexact polynomials", Math. Comp. 74
(2005), 869-903.
x = variable(Float64)
p = (x-1)^4 * (x-2)^3 * (x-3)^2 * (x-4)
q = polyder(p)
gcd(p,q) # should be (x-1)^3 * (x-2)^2 * (x-3), but is a constant
u,v,w,resid = PolynomialZeros.AGCD.agcd(p,q)
LinearAlgebra.norm(u - (x-1)^3*(x-2)^2*(x-3), Inf) ~ 3.8e-6
poly_roots(p) # 2 real, 8 complex
PolynomialZeros.MultRoot.multroot(p) # ([4.0,3.0,2.0,1.0], [1, 2, 3, 4])
01/26/2017
about 1 year ago
51 commits