search

diagrams / diagrams-lib

Diagrams standard library
https://diagrams.github.io/
Other
135 stars 63 forks source link

High-order polynomial root finder #188

Open Mathnerd314 opened 10 years ago

Mathnerd314 commented 10 years ago

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 commented 10 years ago

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

bergey commented 10 years ago

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

kuribas commented 9 years ago

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 commented 9 years ago

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.