Newton iteration for cube root without division

Question

It's a fairly well known trick to avoid division in calculating square-roots to apply Newton's method to finding $1/\sqrt{x}$, and probably better known, using Newton's method to find reciprocals without division.

In rescuing a StackOverflow thread Seeding the Newton iteration for cube root efficiently from link rot, the thought came to me that a division-free iteration for cube roots should also be possible.

For example, if we were to solve:

$$ x^{-3} = a^2 $$

then $x = a^{-2/3}$ and $\sqrt[3]{a} = ax$. The Newton iteration for the above equation is simply:

$$ x_{n+1} = x_n - \frac{ x_n^{-3} - a^2 }{-3x_n^{-4}} = \frac{4}{3}x_n - \frac{1}{3}a^2 x_n^4 $$

Again we avoid division operations, at least if the fractional constants are pre-evaluated for FP multiplications.

So something of the sort is possible, but I did not find a specific discussion of such methods in my (admittedly shallow) search of the Web. More to the point, I suspect that a clever person has already discovered a better idea and that one of you treasured colleagues has seen it or thought it through.

Answer 1

Cube roots are not nearly as important as square roots (e.g., for normalizing vectors), so that might be why they are discussed much less.

In general, if you apply Newton's method to $x^\alpha-\beta$, you get the iteration $$ (1-\alpha^{-1})x+\frac\beta\alpha x^{1-\alpha}, $$ so you just need to pick the equation so that $\alpha$ is negative integer, it's not just cube and square roots, this works for other rational powers as well.

Your scheme only needs 7 iterations to converge to double precision on the interval $a\in[\frac12,1]$ starting from constant initial guess $1$:

Another interesting thing to consider is what existing libraries implement:

• MPFR does cube roots using integer cube roots (http://www.mpfr.org/algo.html, around page 36)
• Boost (http://www.boost.org/doc/libs/1_57_0/libs/math/doc/html/math_toolkit/powers/cbrt.html, see source code math/special_functions/cbrt.hpp) uses Halley iteration with division: $$ x_{n+1} = x_n \frac{2b+x_n^3}{b+2x_n^3}. $$ I tried it, and it seems to be slightly (10-25%, depending on which unsafe floating-point optimizations I enable) faster than division-free Newton's method, when running on my laptop. Given that division is much slower than multiplication/addition (by $\sim20$, see https://software.intel.com/sites/landingpage/IntrinsicsGuide/ and look for latencies of div{s,p}{s,d}), there must be some other process that allows it to run faster.

Edit Another way to generate a good guess is to start with a polynomial approximation to $a^{-2/3}$ before doing Newton iterations. Two-term Chebyshev series reduces number of iterations to 4, three-term requires 3 iterations, and 6-term requires two iterations. In my testing, three-term Chebyshev series followed by 3 iterations was about $\frac13$ faster than just Halley. I haven't tested all possible combinations yet, it seems the fastest so far is 6-term Chebyshev series followed by 1 Halley iteration for $x^3-a$; I also only tested the interval $[\frac12,1]$, to which every floating point number can be reduced by separating out the mantissa.