Implementing an accurate cbrt() function without extra precision

Question

In JavaScript, there is no native cbrt method available. In theory, you could use a method like this:

function cbrt(x) {
  return Math.pow(x, 1 / 3);
}

However, this fails because identities in mathematics don't necessarily apply to floating point arithmetic. For example, 1/3 cannot be accurately represented using a binary floating point format.

An example of when this fails is the following:

cbrt(Math.pow(4, 3)); // 3.9999999999999996

This gets worse as the number gets larger:

cbrt(Math.pow(165140, 3)); // 165139.99999999988

Is there any algorithm which is able to calculate a cube root value to within a few ULP (preferably 1 ULP if possible)?

This question is similar to Computing a correctly rounded / an almost correctly rounded floating-point cubic root, but keep in mind that JavaScript doesn't have any higher-precision number types to work with (there is only one number type in JavaScript), nor is there a built-in cbrt function to begin with.

Answer 1

You can port an existing implementation, like this one in C, to Javascript. That code has two variants, an iterative one that is more accurate and a non-interative one.

Ken Turkowski's implementation relies on splitting up the radicand into mantissa and exponent and then reassembling it, but this is only used to bring it into the range between 1/8 and 1 for the first approximation by enforcing a binary exponent between -2 and 0. In Javascript, you can do this by repeatedly dividing or multiplying by 8, which should not affect accuracy, because it is just an exponent shift.

The implementation as shown in the paper is accurate for single-precision floating-point numbers, but Javascript uses double-precision numbers. Adding two more Newton iterations yields good accuracy.

Here's the Javascript port of the described cbrt algorithm:

Math.cbrt = function(x) 
{    
    if (x == 0) return 0;
    if (x < 0) return -Math.cbrt(-x);

    var r = x;
    var ex = 0;

    while (r < 0.125) { r *= 8; ex--; }
    while (r > 1.0) { r *= 0.125; ex++; }

    r = (-0.46946116 * r + 1.072302) * r + 0.3812513;

    while (ex < 0) { r *= 0.5; ex++; }
    while (ex > 0) { r *= 2; ex--; }

    r = (2.0 / 3.0) * r + (1.0 / 3.0) * x / (r * r);
    r = (2.0 / 3.0) * r + (1.0 / 3.0) * x / (r * r);
    r = (2.0 / 3.0) * r + (1.0 / 3.0) * x / (r * r);
    r = (2.0 / 3.0) * r + (1.0 / 3.0) * x / (r * r);

    return r;
}

I haven't tested it extensively, especially not in badly defined corner cases, but the tests and comparisons with pow I have done look okay. Performance is probably not so great.

Answer 2

Math.cbrt has been added to ES6 / ES2015 specification so at least first check to see if it defined. It can be used like:

Math.cbrt(64); //4

instead of

Math.pow(64, 1/3); // 3.9999999999999996

Answer 3

1. You can use formula for pow computation

   x^y     = exp2(y*log2(x))
   x^(1/3) = exp2(log2(x)*1/3)
           = exp2(log2(x)/3)
   

   • base for log,exp can be any but 2 is directly implemented on most FPU's
   • now you divide by 3 ... and 3.0 is represented by FP accurately.

2. or you can use bit search

   1. find the exponent of output (e= ~1/3 of integer part bit count of x)
   2. create appropriate fixed number y (mantissa=0 and exponent=e)

   3. start bin search from MSB bit of y

      • toggle bit to one
      • if (y*y*y>x) toggle bit back to zero

   4. loop #3 with next bit (stop after LSB)

   The result of binary search is as precise as it can be (no other method can beat it) ... you need mantissa-bit-count iterations for this. You have to use FP for computation so conversion of your y to float is just copying mantissa bits and set exponent.

See pow on integer arithmetics in C++