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

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

1 Answer

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);
}

• That works, but why did you use while(true) and break? Wouldn't it work just as well as a do-while loop? – Brian J. Fink Jul 12 '14 at 1:06
• @BrianJ.Fink I wrote it that way as a direct translation of the algorithm I referenced. Also the assignment x = y makes a do-while version clumsy. You are of course free to rewrite it and test that. – mjolka Jul 12 '14 at 3:21
• Thanks for making my algorithm more accurate; unfortunately I also wanted to make it faster. – Brian J. Fink Jul 12 '14 at 15:40
• I rewrote it as a do loop with the assignment at the beginning instead of the end. It works just fine. – Brian J. Fink Jul 12 '14 at 18:11
• – Brian J. Fink Jul 12 '14 at 18:44