Seeding the Newton iteration for cube root efficiently

Question

How can I find the cube root of a number in an efficient way? I think Newton-Raphson method can be used, but I don't know how to guess the initial solution programmatically to minimize the number of iterations.

Answer 1

This is a deceptively complex question. Here is a nice survey of some possible approaches.

Answer 2

In view of the "link rot" that overtook the Accepted Answer, I'll give a more self-contained answer focusing on the topic of quickly obtaining an initial guess suitable for superlinear iteration.

The "survey" by metamerist (Wayback link) provided some timing comparisons for various starting value/iteration combinations (both Newton and Halley methods are included). Its references are to works by W. Kahan, "Computing a Real Cube Root", and by K. Turkowski, "Computing the Cube Root".

metamarist updates the DEC-VAX era bit-fiddling technique of W. Kahan with this snippet, which "assumes 32-bit integers" and relies on IEEE 754 format for doubles "to generate initial estimates with 5 bits of precision":

inline double cbrt_5d(double d) 
{ 
   const unsigned int B1 = 715094163; 
   double t = 0.0; 
   unsigned int* pt = (unsigned int*) &t; 
   unsigned int* px = (unsigned int*) &d; 
   pt[1]=px[1]/3+B1; 
   return t; 
} 

The code by K. Turkowski provides slightly more precision ("approximately 6 bits") by a conventional powers-of-two scaling on float fr, followed by a quadratic approximation to its cube root over interval [0.125,1.0):

/* Compute seed with a quadratic qpproximation */
fr = (-0.46946116F * fr + 1.072302F) * fr + 0.3812513F;/* 0.5<=fr<1 */

and a subsequent restoration of the exponent of two (adjusted to one-third). The exponent/mantissa extraction and restoration make use of math library calls to frexp and ldexp.

Comparison with other cube root "seed" approximations

To appreciate those cube root approximations we need to compare them with other possible forms. First the criteria for judging: we consider the approximation on the interval [1/8,1], and we use best (minimizing the maximum) relative error.

That is, if f(x) is a proposed approximation to x^{1/3}, we find its relative error:

        error_rel = max | f(x)/x^(1/3) - 1 | on [1/8,1]

The simplest approximation would of course be to use a single constant on the interval, and the best relative error in that case is achieved by picking f_0(x) = sqrt(2)/2, the geometric mean of the values at the endpoints. This gives 1.27 bits of relative accuracy, a quick but dirty starting point for a Newton iteration.

A better approximation would be the best first-degree polynomial:

 f_1(x) = 0.6042181313*x + 0.4531635984

This gives 4.12 bits of relative accuracy, a big improvement but short of the 5-6 bits of relative accuracy promised by the respective methods of Kahan and Turkowski. But it's in the ballpark and uses only one multiplication (and one addition).

Finally, what if we allow ourselves a division instead of a multiplication? It turns out that with one division and two "additions" we can have the best linear-fractional function:

 f_M(x) = 1.4774329094 - 0.8414323527/(x+0.7387320679) 

which gives 7.265 bits of relative accuracy.

At a glance this seems like an attractive approach, but an old rule of thumb was to treat the cost of a FP division like three FP multiplications (and to mostly ignore the additions and subtractions). However with current FPU designs this is not realistic. While the relative cost of multiplications to adds/subtracts has come down, in most cases to a factor of two or even equality, the cost of division has not fallen but often gone up to 7-10 times the cost of multiplication. Therefore we must be miserly with our division operations.