algorithm - 计算正确舍入/几乎正确舍入的浮点立方根

Question

假设正确舍入了标准库函数，例如 CRlibm 中的函数可用。那么如何计算 double 输入的正确四舍五入的立方根呢？

^{引用常见问题解答，这个问题不是“[我]面临的实际问题”。这样有点像作业。但是立方根是一种常见的运算，可以想象这个问题是某人面临的实际问题。}

由于“最好的 Stack Overflow 问题中有一些源代码”，这里有一些源代码:

  y = pow(x, 1. / 3.);

上面没有计算出正确舍入的立方根，因为 1/3 不能完全表示为 double。

---

附加说明:

安article描述了如何计算浮点立方根，但推荐的 Newton-Raphson 算法的最后一次迭代必须以更高的精度进行，以便算法计算正确舍入的 double 立方根。这可能是计算它的最佳方法，但我仍在寻找可以利用现有正确舍入的标准化函数的捷径。

C99 包含一个 cbrt() 函数，但不能期望它被正确舍入 or even faithful对于所有编译器。 CRlibm 的设计者本可以选择将 cbrt() 包含在提供的函数列表中，但他们没有这样做。欢迎引用其他库中可用的正确舍入数学函数的实现。

Answer 1

恐怕我不知道如何保证正确四舍五入的 double 立方根，但可以提供一个非常接近正确四舍五入的立方根，正如问题中提到的那样。换句话说，最大误差似乎非常接近 0.5 ulp。

Peter Markstein，“IA-64 和初等函数:速度和精度”(Prentice-Hall 2000)

提供了基于 FMA 的高效技术，用于正确舍入倒数、平方根和倒数平方根，但它没有涵盖这方面的立方根。一般来说，Markstein 的方法需要一个精确到 1 ulp 以内的初步结果。在最后的舍入序列之前。我没有数学基础将他的技术扩展到立方根的四舍五入，但在我看来这在原则上应该是可能的，这是一个有点类似于平方根倒数的挑战。

按位算法很容易计算出正确舍入的根。由于 IEEE-754 舍入模式的连带情况不会发生，因此只需进行计算，直到产生所有尾数位加上一个舍入位。基于二项式定理的逐位平方根算法在非恢复和恢复变体中都是众所周知的，并且已经成为硬件实现的基础。通过二项式定理的相同方法适用于立方根，并且有一篇鲜为人知的论文列出了非恢复实现的细节:

H. Peng，“提取平方根和立方根的算法”，第 5 届 IEEE 国际计算机算术研讨会论文集，第 121-126 页，1981 年。

最好的我可以从试验中看出它对于从整数中提取立方根来说效果很好。由于每次迭代只产生一个结果位，因此速度并不快。对于浮点运算中的应用程序，它的缺点是使用了几个簿记变量，这些变量大约需要最终结果位数的两倍。这意味着需要使用 128 位整数运算来实现 double 立方根。

我下面的 C99 代码基于 Halley's rational method for the cube root它具有三次收敛性，这意味着初始近似值不必非常准确，因为每次迭代中有效数字的数量都会增加三倍。可以以各种方式安排计算。一般来说，将迭代方案安排为在数值上是有利的

new_guess := old_guess + 更正

因为对于足够接近的初始猜测，correction 明显小于 old_guess。这导致立方根的以下迭代方案:

x := x - x * (x³ - a)/(2*x³ + a)

此特殊安排也列在 Kahan's notes on cube root 中.它的另一个优势是可以自然地使用FMA (fused-multiply add)。 .一个缺点是 2*x³ 的计算可能导致溢出，因此至少部分 double 输入域需要参数缩减方案。在我的代码中，基于对 IEEE-754 double 操作数的指数的直接操作，我简单地将参数约简应用于所有 非异常输入。

区间 [0.125,1) 用作主要近似区间。使用多项式 minimax 近似返回 [0.5,1] 中的初始猜测。窄范围有助于在计算的低精度部分使用单精度算术。

我无法证明我的实现的错误范围，但是，针对引用实现(精确到大约 200 位)使用 2 亿个随机测试向量进行测试发现总共有 277 个不正确的舍入结果(因此错误率大约为1.4 ppm)，最大误差为 0.500012 ulps。

double my_cbrt (double a)
{
    double b, u, v, r;
    float bb, uu, vv;
    int e, f, s;

    if ((a == 0.0) || isinf(a) || isnan(a)) {
        /* handle special cases */
        r = a + a;
    } else {
        /* strip off sign-bit */
        b = fabs (a);
        /* compute exponent adjustments */
        b = frexp (b, &e);
        s = e - 3*342;
        f = s / 3;
        s = s - 3 * f;
        f = f + 342;
        /* map argument into the primary approximation interval [0.125,1) */
        b = ldexp (b, s);
        bb = (float)b;
        /* approximate cube root in [0.125,1) with relative error 5.22e-3 */
        uu =                0x1.2f32c0p-1f;
        uu = fmaf (uu, bb, -0x1.62cc2ap+0f);
        uu = fmaf (uu, bb,  0x1.7546e0p+0f);
        uu = fmaf (uu, bb,  0x1.5d0590p-2f);
        /* refine cube root using two Halley iterations w/ cubic convergence */
        vv = uu * uu;
        uu = fmaf (fmaf (vv, uu, -bb) / fmaf (vv, 2.0f*uu, bb), -uu, uu);
        u = (double)uu;
        v = u * u; // this product is exact
        r = fma (fma (v, u, -b) / fma (v, 2.0*u, b), -u, u);
        /* map back from primary approximation interval by jamming exponent */
        r = ldexp (r, f);
        /* restore sign bit */
        r = copysign (r, a);
    }
    return r;
}