Fast and backward stable computation of roots of polynomials in Julia
8 Stars
Updated Last
2 Months Ago
Started In
August 2014

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

Travis Coverage Status

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.

Example 1: Speed up roots

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)

Example 2: Roots of a polynomial of degree 10,000

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)

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