11-2 Integer Cube Root

For cube roots, Newton's method does not work out very well. The iterative formula is a bit complex:

and there is of course the problem of getting a good starting value x₀.

However, there is a hardware algorithm, similar to the hardware algorithm for square root, that is not too bad for software. It is shown in Figure 11-5.

Figure 11-5 Integer cube root, hardware algorithm.

int icbrt(unsigned x) {
   int s; 
   unsigned y, b; 
 
   s = 30; 
   y = 0; 
   while(s >= 0) {            ?// Do 11 times. 
      y = 2*y; 
      b = (3*y*(y + 1) + 1) << s; 
      s = s - 3; 
      if (x >= b) {
         x = x - b; 
         y = y + 1; 
      } 
   } 
   return y; 
} 

The three add's of 1 can be replaced by or's of 1, because the value being incremented is even. Even with this change, the algorithm is of questionable value for implementation in hardware, mainly because of the multiplication y * (y + 1).

This multiplication is easily avoided by applying the compiler optimization of strength reduction to the y-squared term. Introduce another unsigned variable y2 that will have the value of y-squared, by updating y2 appropriately wherever y receives a new value. Just before y = 0 insert y2 = 0. Just before y = 2*y insert y2 = 4*y2. Change the assignment to b to b = (3*y2 + 3*y + 1) << s (and factor out the 3). Just before y = y + 1, insert y2 = y2 + 2*y + 1. The resulting program has no multiplications except by small constants, which can be changed to shift's and add's. This program has three add's of 1, which can all be changed to or's of 1. It is faster unless your machine's multiply instruction takes only two or fewer cycles.

Caution: [GLS1] points out that the code of Figure 11-5, and its strength-reduced derivative, do not work if adapted in the obvious way to a 64-bit machine. The assignment to b can then overflow. This problem can be avoided by dropping the shift left of s from the assignment to b, inserting after the assignment to b the assignment bs = b << s, and changing the two lines if (x >= b) {x = x - b ?to if (x >= bs && b == (bs >> s)) {x = x - bs ?