A couple of years ago I wrote a blog post about finding the closest point on a cubic bezier to a given point in which I showed how to reduce the problem to solving a quintic equation and provided an implementation in C# that finds the roots of the quintic using Laguerre’s Method. I also mentioned at the end of that post that there is a fancier algorithm for solving specifically quintics described in Doyle and McMullen’s “Solving the Quintic by Iteration”  and speculate that it might be the fastest way to find the closest point to a bezier.

Recently I spent the time to actually implement the Doyle and McMullen algorithm in C++ as a numeric algorithm. I don’t know of any other implementations. My code is here.

It is an interesting algorithm. Its existence implies that although there is no general formula, akin to the quadratic formula, that you can plug into to get the roots of a 5th degree equation, and in fact although the roots of most 5th degree equations cannot even be represented as expressions composed of ordinary arithmetic operations plus radicals — you can represent numbers with radical expressions that are arbitrarily close to these roots. Each of these the roots can be defined exactly in terms of the limit a recursive function converges on as the depth of recursion approaches infinity; each step along the way approximates with increasing accuracy the root of an unsolvable quintic expressed using arithmetic and radicals.

Above I specify my implementation is “a numeric algorithm”, I mean as opposed to an implementation using a symbolic math package; I use ordinary floating point numbers. This is an important point:  is from the wilds of number theory, not from a CS text book, and it thus defines an algorithm over real numbers which of course have arbitrary precision. It was an open question to me whether the algorithm has value as a numeric algorithm. The issue I saw is that it is a solution to quintics in the single parameter Brioschi form. The Brioschi form, by magic, collapses every general quintic equation defined by six complex coefficients into a single complex number; it is impossible for such a reduction to not be limited by finite precision. That is, to some extent the finite precision of floating point numbers must determine to which general quintics a numeric implementation of  will return sensible results.

I will leave an analysis of this question to someone who is an actual mathematician or computer scientist but empirically the Doyle and McMullen algorithm does seem to me to have merit as a numeric algorithm at double precision. On randomly generated general quintics with coefficients uniformly distributed between -1000 and 1000, my implementation of solving the quintic by iteration returns results that look to me to be about as good or better than the implementation of Laguerre’s Method from Data Structures and Algorithms in C++ but is about 500 times faster. Both algorithms perform better in terms of speed and correctness when the coefficients are lower. When they are between -10 and 10 my implementation is more accurate and about 100 times faster than the Data Structures and Algorithms code.

My implementation works as follows:

1. Given a general quintic, convert to principal form. I do this conversion as explained here. I found the resultant and solved for the c₁ and c₂ being canceled out using SageMath.
2. Convert the principal quintic to Brioschi form. I followed the explanation here.
3. Find two solutions to the Brioschi form quintic via iteration as described in  pages 32 to 33. The only difficulty here was that I needed to find the derivative to the function g(Z,w). I did this symbolically again via SageMath; however, a speed optimization would be to replace use of an explicit C++ function for g’ and instead evaluate g(Z,w) and g'(Z,w) simultaneously as described in the chapter in Numerical Recipes: The Art of Scientific Computing on polynomials. Also just as a note to anyone else who may want to write an implementation of this algorithm, the wikipedia article here has a mistake in the definition of h(Z,w). Use the original paper. (Or better yet cut-and-paste from Peter Doyle’s macsyma output and overload C++ such that the caret operator is exponentiation, which is what I did)
4. Convert the two roots back to the general quintic.
5. Test both roots. If the error is less than threshold k pass along both roots. If one is less than k pass along the one good root. If both roots yield errors greater than k perform n iterations of Halley’s Method and retest for one or two good roots. If neither root has error less than k pass both along anyway.
6. If you have two good roots v₁ and v₂, perform synthetic division by (zv₁ )(zv₂ ) yielding a cubic and solve the cubic via radicals. If you have only one good root divide the quintic by (z-v) and solve the resulting quartic. I’m using the cubic solving procedure as described in Numerical Recipes and the quartic formula as described here.

My implementation is a header-only C++17 library (“quintic.hpp” in the github repo i linked to above) parametrized on the specific floating point type you want to use. Single precision is not good enough for this algorithm. Double precision works. I didn’t test on long doubles because Visual Studio does not support them.