High-order polynomial root finder

[Mathnerd314]

Many interesting functions in graphics are implemented by solving high-order polynomials. For example, finding the point of a cubic Bezier that is closest to a given point (the "nearest-point" problem) is implemented by isolating the roots of a 6th-order polynomial.

Wikipedia lists a few algorithms. It appears most systems use several of these in sequence, something like the following:

1. The roots are isolated/bracketed into intervals, with at most one root per interval. Sturm sequences appear to be the main method, although most CAS systems have switched to an algorithm based on Vincent's Theorem for speed.
2. The intervals are refined by bisection until sufficiently small, with various algorithms.
3. An approximate root is calculated through some form of iteration, e.g. Newton's method.

[Mathnerd314]

And misc/DKSolve.hs is an implementation of the Durand-Kerner method, which falls into (3) above.

[bergey]

roots seems to provide many of the iterative methods, though probably not the interval splitting algorithms, which are specific to polynomials.

[kuribas]

An other, more recent way to find roots is to use bezier clipping. It has good numerical stability and it already works in the bernstein basis, which could be a good choice for any CAGD program. It's described in chapter 9 http://tom.cs.byu.edu/~557/text/cagd.pdf
I have a version of the algorithm implemented in my cubicbezier library, which doesn't implement the full algorithm yet, but it probably wouldn't be much work to do so. I could make it into a separate package, or you can copy paste the code.

[bergey]

Thanks for the link! In general, I think it's good to get algorithms like this into libraries for reuse, separate from our representations of points and paths and such. I know that's tricky with geometric algorithms, where you at least need some point / vector type.