Generic function for solving n-order polynomial roots in Julia

Question

All,

I've just been starting to play around with the Julia language and am enjoying it quite a bit. At the end of the 3rd tutorial there's an interesting problem: genericize the quadratic formula such that it solves for the roots of any n-order polynomial equation.

This struck me as (a) an interesting programming problem and (b) an interesting Julia problem. Has anyone out there solved this one? For reference, here is the Julia code with a couple toy examples. Again, the idea is to make this generic for any n-order polynomial.

Cheers,

Aaron

function derivative(f)
    return function(x)
        # pick a small value for h
        h = x == 0 ? sqrt(eps(Float64)) : sqrt(eps(Float64)) * x

        # floating point arithmetic gymnastics
        xph = x + h
        dx = xph - x

        # evaluate f at x + h
        f1 = f(xph)

        # evaluate f at x
        f0 = f(x)

        # divide the difference by h
        return (f1 - f0) / dx
    end
end


function quadratic(f)

    f1 = derivative(f)

    c = f(0.0)

    b = f1(0.0)

    a = f(1.0) - b - c

    return (-b + sqrt(b^2 - 4a*c + 0im))/2a, (-b - sqrt(b^2 - 4a*c + 0im))/2a
end

quadratic((x) -> x^2 - x - 2)
quadratic((x) -> x^2 + 2)

Answer 1

The package PolynomialRoots.jl provides the function roots() to find all (real and complex) roots of polynomials of any order. The only mandatory argument is the array with coefficients of the polynomial in ascending order.

For example, in order to find the roots of

6x^5 + 5x^4 + 3x^2 + 2x + 1

after loading the package (using PolynomialRoots) you can use

julia> roots([1, 2, 3, 4, 5, 6])
5-element Array{Complex{Float64},1}:
           0.294195-0.668367im
 -0.670332+2.77556e-17im
           0.294195+0.668367im
          -0.375695-0.570175im
          -0.375695+0.570175im

The package is a Julia implementation of the root-finding algorithm described in this paper: http://arxiv.org/abs/1203.1034

PolynomialRoots.jl has also support for arbitrary precision calculation. This is useful for solving equation that cannot be solved in double precision. For example

julia> r = roots([94906268.375, -189812534, 94906265.625]);

julia> (r[1], r[2])
(1.0000000144879793 - 0.0im,1.0000000144879788 + 0.0im)

gives the wrong result for the polynomial, instead passing the input array in arbitrary precision forces arbitrary precision calculations that provide the right answer (see https://en.wikipedia.org/wiki/Loss_of_significance):

julia> r = roots([BigFloat(94906268.375), BigFloat(-189812534), BigFloat(94906265.625)]);

julia> (Float64(r[1]), Float64(r[2]))
(1.0000000289759583,1.0)