FastPolynomialRoots.jl - Fast and backward stable computation of roots of polynomials

This package is a Julia wrapper of the Fortran programs accompanying Fast and Backward Stable Computation of Roots of Polynomials by Jared L. Aurentz, Thomas Mach, Raf Vandebril and David S. Watkins.

The package provides the unexported function FastPolynomialRoots.rootsFastPolynomialRoots(p::Vector{<:Union{Float64,Complex{Float64}}}) which computes the roots of the polynomial p[1] + p[2]*x + p[3]*x^2 + ... + p[k]*x^(k-1). The package also overwrites the roots(::Polynomial) methods in the Polynomials package for Float64 and Complex{Float64} elements with the fast versions provided by this package. See the examples below.

julia> using Polynomials, BenchmarkTools

julia> @btime roots(p) setup=(p = Polynomial(randn(500)));
  223.135 ms (23 allocations: 3.97 MiB)

julia> using FastPolynomialRoots

julia> @btime roots(p) setup=(p = Polynomial(randn(500)));
  30.786 ms (7 allocations: 26.41 KiB)

A computation of this size would not be feasible on a desktop with the traditional method but can be handled by FastPolynomialRoots.

julia> using Polynomials, BenchmarkTools, FastPolynomialRoots

julia> n = 10000;

julia> r = @btime roots(p) setup=(p = Polynomial(randn(n + 1)));
  10.290 s (13 allocations: 508.38 KiB)

julia> sum(isreal, r)
7

julia> 2/π*log(n) + 0.6257358072 + 2/(n*π) # Edelman and Kostlan
6.489284260212659