Qantas 94 Heavy Qantas 94 Heavy - 4 months ago 8
Javascript Question

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

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.

Comments