3

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.

Community
  • 1
  • 1
Qantas 94 Heavy
  • 15,750
  • 31
  • 68
  • 83
  • Your question is very close to http://stackoverflow.com/questions/18063755/computing-a-correctly-rounded-an-almost-correctly-rounded-floating-point-cubic . I can only reiterate the trick I discovered while thinking about it that while 1/3 is not available to represent as a double, 3 is, so you can `pow(candidate, 3)` and check that against the original value (assuming `pow()` is of good quality) (ideally you'd have a correctly rounded `pow_ru()` and one bit of extra precision). – Pascal Cuoq May 01 '14 at 08:32
  • @PascalCuoq: it's indeed rather close. The only thing about that question is that it seems to concern C99, which has options such as long doubles to be able to get higher precision if I'm not wrong. The thing about JS is that it doesn't have that extra bit of precision owing to its single number type, which makes things a bit more difficult. – Qantas 94 Heavy May 01 '14 at 08:39

3 Answers3

2

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.

M Oehm
  • 28,726
  • 3
  • 31
  • 42
1

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

bkaid
  • 51,465
  • 22
  • 112
  • 128
0

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

Community
  • 1
  • 1
Spektre
  • 49,595
  • 11
  • 110
  • 380
  • 1- There are in the literature descriptions of okay-ish implementations of `pow()` as the composition of `exp()` and `log()`, but these descriptions point out that you need additional precision for the intermediate results if you intend the end result to be reasonably accurate. In fact, I think that the description I saw of such a `pow()` implementation was in a discussion of the 8087 extended-double format with its 64 bits of precision, that makes such an approach appropriate for double-precision `pow()`. – Pascal Cuoq May 01 '14 at 08:43
  • 2- You don't need to start from a rough approximation of the result where only the exponent is correct. You can use `pow(x, 1.0/3.0)` as starting point. – Pascal Cuoq May 01 '14 at 08:46
  • @Pascal Cuoq i8087 works on 80bits internally if I remember correctly and yes you can use that as start point – Spektre May 01 '14 at 10:21
  • 2
    80 bits is the total size of the 8087 extended format. Of these, 15 are used for the exponent and 1 for the sign, leaving 64 bits of precision. Anyway, the question is in the context of Javascript, which does not provide access to any extended format. – Pascal Cuoq May 01 '14 at 11:02
  • @PascalCuoq: "You can use pow(x, 1.0/3.0) as starting point". Only for positive `x`. For example, -27^(1/3.0) is `nan`. – J D Nov 16 '15 at 00:38