Implementing An Accurate Cbrt() Function Without Extra Precision

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

Solution 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.

Solution 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

Solution 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++