无需额外精度即可实现准确的 cbrt() 函数

Question

在 JavaScript 中，没有可用的原生 cbrt 方法。理论上，你可以使用这样的方法：

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

但是，这失败了，因为数学中的恒等式不一定适用于浮点运算。例如，1/3 不能用二进制浮点格式准确表示。

失败的例子如下：

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

随着数量的增加，情况会变得更糟：

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

是否有任何算法能够在几个 ULP（如果可能的话，最好是 1 个 ULP）内计算立方根值？

这个问题类似于Computing a correctly rounded / an almost correctly rounded floating-point cubic root，但请记住，JavaScript 没有任何更高精度的数字类型可以使用（JavaScript 中只有一种数字类型），也没有内置的 @ 987654327@ 函数开头。

Answer 1

您可以将现有实现（例如 this one in C）移植到 Javascript。该代码有两种变体，一种更准确的迭代代码和一种非交互代码。

Ken Turkowski 的实现依赖于将 radicand 拆分为尾数和指数，然后重新组装它，但这仅用于通过强制在 -2 之间的二进制指数将其带入 1/8 和 1 之间的范围以进行第一次近似和 0。在 Javascript 中，您可以通过反复除以或乘以 8 来做到这一点，这不会影响准确性，因为它只是一个指数移位。

论文中所示的实现对于单精度浮点数是准确的，但 Javascript 使用的是双精度数。再添加两次牛顿迭代会产生良好的准确性。

这是所描述的cbrt算法的Javascript端口：

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

我没有对它进行过广泛的测试，尤其是在定义不明确的极端情况下，但我所做的测试和与pow 的比较看起来还不错。性能可能不是那么好。

Answer 2

Math.cbrt 已添加到 ES6 / ES2015 规范中，因此至少首先检查它是否已定义。它可以像这样使用：

Math.cbrt(64); //4

而不是

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

Answer 3

1. 您可以使用pow computation的公式

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

   • log,exp 的基数可以是任何但2 直接在大多数 FPU 上实现
   • 现在你除以 3 ... 并且 3.0 由 FP 准确表示。

2. 或者你可以使用位搜索

   1. 求输出(e= ~1/3 of integer part bit count of x)的指数
   2. 创建适当的固定数字 y（mantissa=0 和 exponent=e）

   3. 从y的MSB位开始bin搜索

      • 将位切换为 1
      • 如果(y*y*y>x) 将位切换回零

   4. 循环 #3 下一位（在 LSB 之后停止）

   二进制搜索的结果尽可能精确（没有其他方法可以击败它）......为此，您需要尾数位计数迭代。您必须使用 FP 进行计算，因此将 y 转换为 float 只是复制尾数位并设置指数。

见pow on integer arithmetics in C++