8
$\begingroup$

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.

$\endgroup$
8
$\begingroup$

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$:

enter image description here

Another interesting thing to consider is what existing libraries implement:

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.

$\endgroup$
  • $\begingroup$ 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). $\endgroup$ – hardmath Feb 5 '15 at 13:43
  • $\begingroup$ @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. $\endgroup$ – Kirill Feb 5 '15 at 22:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.