Newton iteration for cube root without division

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.

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$:
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.
• Thanks, I appreciate the discussion of Halley's method, and the consideration of starting values. However I think for interval reduction the cube root would naturally require $[1/8,1]$ or similar multiplied-by-8 interval (since a binary exponent can be reduced by a multiple of 3). – hardmath Feb 5 '15 at 13:43
• @hardmath If the exponent is not divisible by $3$, you end up multiplying the result by a known constant like $2^{1/3}$. Also, a longer interval might require more iterations to guarantee convergence everywhere in it, so it needs testing. Plus, $[\frac12,1]$ is very natural for floating-point numbers. – Kirill Feb 5 '15 at 22:55