Finding roots of polynomials in Haskell?

tl;dr: There are several Haskell packages one can use to find an individual root of a function on a certain interval. But I’ve had less luck finding a suitable package for finding all the roots of a given polynomial. This blog post consists of my notes and a plea for help.

The diagrams-solve package contains miscellaneous numeric solving routines used in diagrams, including tridiagonal and cyclic tridiagonal linear system solvers (used for generating cubic splines and emulating Metafont paths) as well as functions for finding roots of low-dimensional polynomials (quadratic, cubic, and quartic). Solving quadratics is used in a bunch of places; solving cubics is needed specifically for doing interior/exterior testing on closed loops built from cubic Bezier curves; thankfully we have never needed to solve a quartic or higher.

Unfortunately, the polynomial root solvers are pretty naive: I simply transcribed some formulas (which is why it only goes up to quartics). This works OK a lot of the time, but the formulas are very numerically unstable, and it’s not hard to come up with example inputs where the returned roots are off by quite a bit. In fact, people regularly run into this when running the test suite. I am not specifically aware of any diagrams bugs that have arisen in actual practice due to the cubic solutions being off, but it’s probably just a matter of time.

So I decided it was finally time to look into better root-finding methods. This blog post is both a plea for help and a place to write down some of the things I’ve learned so far.

There’s root finding, and then there’s root finding

The first thing to get out of the way is that when you talk about “root finding” there are (at least!) two pretty different things you could mean:

1. Given a function \(f\) and some particular interval, or an initial guess, find a value \(x\) in the interval/close to the initial guess for which \(f(x) = 0\).
2. Given a polynomial with {real, complex} coefficients, find all its {real, complex} roots.

If you want to do (1), there are several nice Haskell packages you could use. The Numeric.RootFinding module from the math-functions package is probably your best bet; it implements both Ridders’ method and the Newton-Raphson method, which both attempt to find a single root of a function on a given interval. They both work on any continuous Double -> Double function, not just polynomials (Newton-Raphson also needs to know the first derivative). But they don’t work well if you don’t already have an interval in mind to search; and if you want to find all the roots you have to call these multiple times (somehow coming up with an appropriate interval each time).

As for (2), I haven’t been able to find anything that would work well for diagrams-solve. Here are my notes:

• The dsp package has a module Polynomial.Roots containing an implementation of Laguerre’s method, which finds all the (complex) roots of a polynomial. However, the dsp package is a rather heavyweight dependency to pull in just for the root-finding code; it’s also licensed under the GPL and I’d like to avoid having to “infect” the entire diagrams ecosystem with the GPL.

• Laguerre’s method seems like it should be fairly easy to implement myself—but writing my own solver from scratch is how I got here in the first place; I’d really like to avoid it if possible. I am far from being an expert on numerical analysis, floating-point computation, etc.

• The hmatrix-gsl package has the Numeric.GSL.Polynomials module, which has an interface to a root finding algorithm from GSL (apparently using something called “balanced-QR reduction”), but I’d like to avoid pulling in a C library as as dependency, and also, again, GPL.

• From the Wikipedia page for Laguerre’s method I learned that the Jenkins-Traub algorithm is another widely used method for polynomial root-finding, and often preferred over Laguerre’s method. However, it seems rather complex, and the only Haskell implementation of Jenkins-Traub I have been able to fnid is this one which seems to be just a toy implementation; I don’t even know if it works correctly.

If you know of a good place where I can find polynomial solving code in Haskell, can you point me to it? Or if you know more about numerics than me, could you maybe whip up a quick implementation of Laguerre’s method and put it up on Hackage?