Qantas 94 Heavy - 9 months ago 22

Javascript Question

In JavaScript, there is no native

`cbrt`

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

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.