Followup: How do I optimize this Java cube root function for BigInteger?

Question

Followup to How do I optimize this Java cube root function for BigInteger?

So I've tried several implementations of this algorithm. The version using only BigInteger sometimes results in a never-ending cycle of candidates:

public static final BigInteger _3=BigInteger.valueOf(3);

public static BigInteger cbrt(BigInteger num)
{
 BigInteger root=BigInteger.ZERO.setBit(num.bitLength()/3),temp;
 do {
  temp=root;
  root=temp.add(temp).add(num.divide(temp.multiply(temp))).divide(_3);
 } while(!root.equals(temp));
 return root;
}

This hangs the program for certain numbers (5 and 26 are prime examples).

Next I tried using BigDecimal to give me more accurate division, but I was afraid that converting back and forth between BigInteger and BigDecimal might be slowing it down:

public static final BigDecimal _3=BigDecimal.valueOf(3);
public static final int UP=BigDecimal.ROUND_HALF_UP;

public static BigInteger cbrt(BigInteger num)
{
 BigDecimal numD=new BigDecimal(num),temp;
 BigInteger root=BigInteger.ZERO.setBit(num.bitLength()/3);
 do {
  temp=new BigDecimal(root);
  root=temp.add(temp).add(numD.divide(temp.multiply(temp),UP)).divide(_3,UP).toBigInteger();
 } while(!root=temp.toBigInteger());
 return root;
}

On a suggestion from the comments, I tried modifying the first function above with the algorithm ending with the root less than or equal to the temp. But that returns the wrong root for some numbers:

public static final BigInteger _3=BigInteger.valueOf(3);

public static BigInteger cbrt(BigInteger num)
{
 BigInteger root=BigInteger.ZERO.setBit(num.bitLength()/3),temp;
 do {
  temp=root;
  root=temp.add(temp).add(num.divide(temp.multiply(temp))).divide(_3);
 } while(root.compareTo(temp)>0);
 return root;
}

Try 1367631. The return value is 113, not 111.

Then there's my solution: multiple temporary values:

public static final BigInteger _3=BigInteger.valueOf(3);

public static BigInteger cbrt(BigInteger n)
{
 BigInteger root=BigInteger.ZERO.setBit(n.bitLength()/3),t1,t2,t3;
 t1=t2=t3=BigInteger.ZERO;
 do {
  t3=t2;t2=t1;t1=root;
  root=t1.add(t1).add(n.divide(t1.multiply(t1))).divide(_3);
 } while(!root.equals(t1.add(t2).add(t3).divide(_3)));
 return root;
}

I'm trying to decide between the BigDecimal version and the multiple-temp version. Any suggestions how I can make either algorithm better?

Edit: Followup can be found here: https://codereview.stackexchange.com/questions/56862

Answer 1

I think your final solution is on the right path.

Just to be clear, I'm assuming you want the integral part of the cube root, i.e. for a given \$n\$, cbrt(n) will return the integer \$m\$ such that \$m^3 \leq n < (m + 1)^3\$. I'm also assuming you only care about positive inputs.

Here is the code I posted in a comment on the previous question, with an improved initial estimate. Inspired by There are Only Four Billion Floats -- So Test Them All!, I have tested it successfully on all integers in the range [1, Integer.MAX_VALUE) (took about 40m, by the way).

private static final BigInteger THREE = BigInteger.valueOf(3);

private static BigInteger cubeRoot(BigInteger n) {
    // Using Newton's method, we approximate the cube root
    // of n by the sequence:
    // x_{i + 1} = \frac{1}{3} \left( \frac{n}{x_i^2} + 2 x_i \right).
    // See http://en.wikipedia.org/wiki/Cube_root#Numerical_methods.
    //
    // Implementation based on Section 1.7.1 of
    // "A Course in Computational Algebraic Number Theory"
    // by Henri Cohen.
    BigInteger x = BigInteger.ZERO.setBit(n.bitLength() / 3 + 1);
    while (true) {
        BigInteger y = x.shiftLeft(1).add(n.divide(x.multiply(x))).divide(THREE);
        if (y.compareTo(x) >= 0) {
            break;
        }

        x = y;
    }

    return x;
}

The test code:

for (int i = 1; i < Integer.MAX_VALUE; i++) {
    // Report progress.
    if (i % 1000000 == 0) {
        System.out.printf("%d%n", i);
    }

    BigInteger n = BigInteger.valueOf(i);
    BigInteger m = cubeRoot(n);

    BigInteger lower = m.pow(3);
    BigInteger upper = m.add(BigInteger.ONE).pow(3);

    if (lower.compareTo(n) <= 0 && n.compareTo(upper) < 0) {
        continue;
    }

    System.err.printf("Error for input %s: Got %s%n", n, m);
    System.err.printf("Expected m^3 <= %s < (m + 1)^3%n", n);
    System.err.printf("But m^3 = %s, (m + 1)^3 = %s%n", lower, upper);
}