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 diagramssolve
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 lowdimensional 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 rootfinding 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:
 Given a function and some particular interval, or an initial guess, find a value in the interval/close to the initial guess for which .
 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 mathfunctions
package is probably your best bet; it implements both Ridders’ method and the NewtonRaphson 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 (NewtonRaphson 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 diagramssolve
. Here are my notes:

The
dsp
package has a modulePolynomial.Roots
containing an implementation of Laguerre’s method, which finds all the (complex) roots of a polynomial. However, thedsp
package is a rather heavyweight dependency to pull in just for the rootfinding 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, floatingpoint computation, etc.

The
hmatrixgsl
package has theNumeric.GSL.Polynomials
module, which has an interface to a root finding algorithm from GSL (apparently using something called “balancedQR 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 JenkinsTraub algorithm is another widely used method for polynomial rootfinding, and often preferred over Laguerre’s method. However, it seems rather complex, and the only Haskell implementation of JenkinsTraub 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?
It probably doesn’t apply to your case, but just In case you want to find all roots of an arbitrary function, there is a MATLAB package called `chebfun` that offers interpolation by Chebychev polynomials, which I found highly useful. The roots of Chebychev polynomials can be found by means of a Fourier transform, but that is also hidden away in the package.
Yeah, that doesn’t really help me right now, but good to know for future reference!
chebfun is indeed a great package – I don’t think it would be too hard to implement some of it in Haskell (I even started doing this several years ago) – the clever idea is that because of the precision of floating point you can approximate *any* (realvalued) function to this precision by a (Chebyshev) polynomial.
I don’t know nearly enough Haskell to do it myself, but https://en.wikipedia.org/wiki/Durand–Kerner_method looks promising for finding the entire root set and easier to implement than Laguerre.
But honestly what I think I’d do in that situation is bite the bullet about foreign dependencies and bring in hmatrix (NOT hmatrixgsl). Then to find the roots of a polynomial, build the companion matrix : https://en.wikipedia.org/wiki/Companion_matrix and use hmatrix to compute the eigenvalues, which are exactly the roots.
Yeah, I tried implementing DurandKerner once but never really got it to work. I might try again.
Thanks for the advice re: hmatrix. I may end up doing that!
I also implemented DurandKernerWeierstraß when I was in highschool. The implementation was in C and had to implement complex arithmetics and polynomials including input and output of them. Altogether about 300 lines. I couldn’t find any polynomial that wouldn’t very rapidly yield all roots (maybe I was just lucky). The code is lost, of course, long ago; it didn’t seem savingworthy since it was no more than two hours of work. I would give DKW a shoot.
Converting to an eigenvalue problem is now standard practice. A monomial basis representation leads to companion matrices. If you have an orthogonal basis representation then it is better not convert to monomials and use colleague/comrade matrices (this is the approach used in the Chebfun package already mentioned). Some names to look up: IJ Good, JP Boyd, V Noferini, DS Watkins, R Vandebril, J Aurentz.
I found an article in the Monad.Reader that discusses a Haskell implementation of the Jacobi method for finding the eigenvalues of a matrix (using the Repa package) [1]. It seems an implementation is available [2].
The other standard algorithms for finding all eigenvalues are QR and divideandconquer, but getting all the details right can be quite involved.
[1] https://themonadreader.files.wordpress.com/2013/03/issue214.pdf
[2] https://github.com/felipeZ/Haskellabinitio/tree/master/Science/QuantumChemistry/NumericalTools
Just reading through the comments, you might be able to do some quick experiments using hmatrix. It can find eigenvalues (probably via LAPACK), Chebyshev polynomials (⚠️ – I never tried this) and roots. The latter two are in hmatrixgsl which is a binding to GSL (GNU Scientific Library) which comes with its GPL!
The implementation of Laguerre’s method in the DSP package is quite compact [1,2].
Maybe the original author would be open to having it in a separate package (possibly using a more permissive license).
[1] https://hackage.haskell.org/package/dsp0.2.1/docs/src/PolynomialRoots.html
[2] http://hackage.haskell.org/package/dsp0.2.4.1/docs/src/Polynomial.Basic.html
After a bit of digging I found that it already is in a separate package: http://hackage.haskell.org/package/numericquest . Still GPL of course. No idea whether the original author (Jan Sibinski) can even be contacted, the code is from 20 years ago (and was rescued via the Wayback Machine) and I haven’t been able to find them anywhere online.
I thought the Alberth method looked quite elegant, so I wrote a basic implementation. This is not meant to be robust and it is not difficult to come up with situations where this basic implementation fails.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
alberth_method.hs
hosted with ❤ by GitHub
Oops. I thought that would be a link, but it seems WordPress embedded the whole Gist. Feel free to remove or change the comment.
Nah, no worries, it’s fine. Thanks for sharing!
Pingback: What’s the right way to QuickCheck floatingpoint routines?  blog :: Brent > [String]