## 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?

Assistant Professor of Computer Science at Hendrix College. Functional programmer, mathematician, teacher, pianist, follower of Jesus.
This entry was posted in haskell, math and tagged , , , , , . Bookmark the permalink.

### 14 Responses to Finding roots of polynomials in Haskell?

1. Heinrich Apfelmus says:

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.

• Brent says:

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* (real-valued) function to this precision by a (Chebyshev) polynomial.

2. Jake says:

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 hmatrix-gsl). 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.

• Brent says:

Yeah, I tried implementing Durand-Kerner 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 Durand-Kerner-Weierstraß 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 saving-worthy since it was no more than two hours of work. I would give DKW a shoot.

3. Jan Van lent says:

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 divide-and-conquer, but getting all the details right can be quite involved.

4. 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 hmatrix-gsl which is a binding to GSL (GNU Scientific Library) which comes with its GPL!

5. Jan Van lent says:

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).

• Brent says:

After a bit of digging I found that it already is in a separate package: http://hackage.haskell.org/package/numeric-quest . 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.

6. Jan Van lent says:

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.

 — Alberth method for finding all roots of a (complex) polynomial — The calculations have to use complex numbers. — Some optimizations are possible if the coefficients are real, — This is not implemented here. — Implementation was prompted by — This is not meant to be robust implementation, just a proof of concept. — For analysis and implementation of this specific method, — see the work by D Bini and others. — — Jan Van lent — Started: 2019-02-20 Wed — Last updated: 2019-02-21 Thu import Data.Complex — |Given [p0, p1, …, pn] x, evaluate polynomial p0 + p1*x + … + pn*x^n polyval :: Num a => [a] -> a -> a polyval p x = rec p 1 0 where rec [] z y = y rec (p0:p) z y = rec p (x*z) (y + p0*z) — alternative: Horner scheme, sparse polynomials, FFT, etc. diff :: Num a => [a] -> [a] diff p = zipWith (*) (map fromIntegral [1..]) (tail p) scales :: Fractional a => [a] -> [a] scales zs = [ sum [ 1/(zk–zj) | (zj, j) <- zip zs [1..], j /= k ] | (zk, k) <- zip zs [1..] ] offset :: Fractional a => a -> a -> a offset q s = –q/(1–q*s) offsets :: Fractional a => [a] -> [a] -> [a] -> [a] offsets fs gs ss = ws where qs = zipWith (/) fs gs ws = zipWith offset qs ss solve' :: RealFloat a => [Complex a] -> [Complex a] -> [Complex a] -> Int -> a -> [Complex a] solve' p p' zs maxit tol | maxit <= 0 = zs — | maxf <= tol = zs | maxw <= tol = newzs | otherwise = solve' p p' newzs (maxit–1) tol where fs = map (polyval p) zs maxf = maximum (map magnitude fs) gs = map (polyval p') zs ss = scales zs ws = offsets fs gs ss maxw = maximum (map magnitude ws) newzs = zipWith (+) zs ws solve :: RealFloat a => [Complex a] -> [Complex a] -> Int -> a -> [Complex a] solve p zs maxit tol = solve' p (diff p) zs maxit tol rootBound :: RealFloat a => [Complex a] -> a rootBound p = r where revp = reverse p num = maximum (map magnitude (tail revp)) den = magnitude (head revp) r = 1 + num / den — choose initial points equally spaced on a circle with radius based on bounds — for the absolute values of the roots — for a polynomial with real coefficients it would make sense to start — with conjugate pairs, since this is preserved by the iteration — I suspect the literature on this method will have much more to say — about how to choose good initial values circleInit :: RealFloat a => [Complex a] -> [Complex a] circleInit p = [ mkPolar (r/2) (((fromIntegral i)+0.5)*step) | i <- [0..n–1] ] where step = 2*pi/(fromIntegral n) r = rootBound p n = (length p) – 1 phi = ((sqrt 5) + 1)/2 goldenAngle = 2 * pi / phi**2 — similar to circleInit — the golden angle is used to generate a quasi-uniform distribution — this is still deterministic, but a bit more like random goldenCircleInit :: [Complex Double] -> [Complex Double] goldenCircleInit p = [ mkPolar (r/2) ((fromIntegral i)*goldenAngle) | i <- [0..n–1] ] where r = rootBound p n = (length p) – 1 roots :: [Complex Double] -> Int -> Double -> [Complex Double] —roots p tol = solve p (circleInit p) tol roots p maxit tol = solve p (goldenCircleInit p) maxit tol rootpoly' :: Num a => [a] -> [a] -> [a] rootpoly' [] p = p rootpoly' (z:zs) p = rootpoly' zs (zipWith (-) (0:p) ((map (z*) p)++[0])) — |find the coefficients of a polynomial with given roots rootpoly :: Num a => [a] -> [a] rootpoly zs = rootpoly' zs [1] test1 = (z0s, ss, ws) where —p = [3:+0, 2:+0, 1:+0] p = [3, 2, 1] z0s = goldenCircleInit p p' = diff p ss = scales z0s fs = map (polyval p) z0s gs = map (polyval p') z0s ws = offsets fs gs ss test2 = zs where p = rootpoly [1:+0, 2:+0, 3:+0] zs = roots p 30 1e-3 test3 n = zs where p = rootpoly [ i:+i^2 | i <- [1..n] ] zs = roots p 30 1e-3

• Jan Van lent says:

Oops. I thought that would be a link, but it seems WordPress embedded the whole Gist. Feel free to remove or change the comment.

• Brent says:

Nah, no worries, it’s fine. Thanks for sharing!

This site uses Akismet to reduce spam. Learn how your comment data is processed.