Quintic roots via Halley

In the previous post, I posted a couple of C++ template functions to approximate n^th roots via a bit hack. The double case can be optimized using a method such as Kahan's, but both the 32-bit and 64-bit cases appear to be working.

This time, I'm going to post a little function that uses Halley's method to calculate the quintic root of a real number.

inline double halley_quint1d(double x, double R)
{
// iterative quintic root approximation using Halley's method (double)const double x5 = x*x*x*x*x;
return x*(3*R + 2*x5) / (2*R + 3*x5);
}

Given a guess 'x' for the quintic root of R, this function will return an approximation that is closer to the quintic root of R than 'x', and the function will exhibit cubic convergence.

You can combine calls to this function with the previously posted nth_roota( ) function to generate an approximation of a quintic root. The iteration, halley_quint1d( ), can be called repeatedly to a desired degree of accuracy.

For example:

double quintica(double R)
{
return halley_quint1d(nth_roota<5>(R), R);
}

Of course one could also use Newton-Raphson, but I'll leave that as an exercise to the reader.

Lazily, I used wxMaxima to work out the simplified version of Halley's method applied to quintic roots. If you'd like to verify it yourself, fire up wxMaxima (an excellent free program I highly recommend) and enter the following:

(%i1) x^5-R;
(%i2) x - 2*%o1*diff(%o1,x)/(2*diff(%o1,x)^2-%o1*diff(%o1,x,2));
(%i3) ratsimp(%);