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