9

假设可以使用正确舍入的标准库函数，例如CRlibm中的函数。那么如何计算双精度输入的正确舍入三次根呢？

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

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

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

上面没有计算正确舍入的三次根，因为 1/3 不能完全表示为double.

补充说明：

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

C99 包含一个cbrt()函数，但不能期望它对所有编译器都正确舍入甚至忠实。CRlibm 的设计者本可以选择包含cbrt()在提供的函数列表中，但他们没有。欢迎参考其他正确舍入数学函数库中可用的实现。

4

3 回答 3

5

鉴于曲线 x = y^3 上有许多易于计算的有理点，我很想减少大约 s^3 ~ x，其中 s 有理且只有几位宽。然后你有：

cbrt(x) = s * cbrt(1 + (x - s^3)/s)

显而易见的事情是使用您最喜欢的级数近似来评估校正项，并通过头尾 FMA 算法计算残差，以便在需要时将结果向上或向下提高或降低 ulp（大多数情况下您不需要完整计算， 明显地）。

这不完全符合问题的精神，但绝对可以工作，并且很容易以这种方式证明必要的界限。希望其他人可以提出更聪明的建议（我已经用光了这个月的聪明才智）。

于 2013-08-05T17:30:10.787 回答
3

恐怕我不知道如何保证正确舍入的双精度立方根，但可以提供一个非常接近正确舍入的值，如问题中所述。换句话说，最大误差似乎非常接近 0.5 ulp。

Peter Markstein，“IA-64 和基本函数：速度和精度”（Prentice-Hall 2000）

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

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

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

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

下面我的 C99 代码基于Halley 的有理法立方根，它具有三次收敛性，这意味着初始近似值不必非常准确，因为每次迭代中有效数字的数量都会增加三倍。可以以各种方式安排计算。通常，将迭代方案安排为在数值上是有利的

new_guess := old_guess + 更正

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

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

这种特殊的排列方式也列在Kahan 关于立方根的笔记中。它的另一个优势是自然地使用FMA（融合乘加）。一个缺点是 2*x 3的计算可能导致溢出，因此至少部分双精度输入域需要参数缩减方案。在我的代码中，我只是将参数归约应用于所有非异常输入，基于对 IEEE-754 双精度操作数的指数的直接操作。

区间 [0.125,1) 用作主要近似区间。使用多项式极小极大近似，返回 [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;
}
于 2014-11-19T03:40:29.617 回答
3

我最近编写了一个正确舍入的立方根，在https://github.com/mockingbirdnest/Principia/blob/master/numerics/cbrt.cpp中实现并记录在https://github.com/mockingbirdnest/Principia/blob /master/documentation/cbrt.pdf。该实现是在 C++ 中并使用 Intel 内在函数，但它应该很容易将其调整为一个人选择的语言和库。

下面是实际计算的简要总结；此处的参考文献使用了该文档中的书目代码，因为我显然声名狼藉，无法在此时链接许多东西。对“部分I+”或“附录[A-Z]”的引用是 cbrt.pdf。

正确的四舍五入。

这可以通过使用几乎正确的忠实方法来完成，如果结果太接近平局，则使用一轮经典的逐位算法（类似于长除法）得到第 54 位；由于该算法计算的结果位向 0 舍入，并且由于立方根没有中途情况，因此该位足以确定舍入到最接近的值。

有可能得到几乎正确的方法，它比大多数现有的立方根实现更快，并且足够正确，平均成本基本上不受校正步骤的影响；见附录 D 和 F。

关于没有 FMA 的忠实方法的评论。

“几乎正确的忠实方法”部分可能更有趣；没有 FMA 的基本轮廓与 Kahan 的 [KB01] 中的相同：

  1. 给定浮点数表示的整数运算；
  2. 寻根精度为三分之一；
  3. 四舍五入到精度的三分之一，以获得精确的立方体；
  4. 使用该精确立方体的高阶方法的一个步骤。

在步骤 2 中使用巧妙的求根方法，可以降低错误取整率，同时保持或提高性能。在这种情况下，我挖了一个不合理的方法，

↦ ½ + √(¼² + /(3)) 近似为 ∛(³+)，

摘自 Thomas Fantet de Lagny 的 17 世纪著作，[Fan92]；也许令人惊讶，因为它同时具有除法和平方根，它在适当重写时的性能（以便尽早安排除法，并避免平方根和除法之间的串行依赖，参见附录 D）类似于更广为人知的有理方法（现在通常称为哈雷法，但哈雷考虑了这两种方法，当应用于立方根时，两者都是由于拉格尼；参见第一部分的讨论）。这是因为有理方法有一个更高阶的除数，所以它的除数不能提前安排。这种非理性方法的错误是理性方法的一半。

优化系数以最小化步骤 2 结束时的误差，以斯蒂芬佳能（两个推特线程，[Can18a] 和 [Can18b]）的推文启发的方式，在该非理性方法的错误上再获得两位，在没有成本。

使用第 4 步中的 5 阶合理方法，该方法实现了每百万 4.564(68) 的错误取整率，这很容易在基本上没有平均成本的情况下被纠正（在缓慢的“潜在错误取整”路径中的通过率为 2.6480（ 52)⋅10 -4，慢速路径的延迟小于快速路径的十倍)。

关于 FMA 的忠实方法的评论。

关键思想在于njuffa对这个问题的回答：步骤 4 可以只用一个精确的正方形来执行，而不需要一个精确的立方体，所以步骤 2 应该瞄准，而步骤 3 应该四舍五入到一半的精度，而不是超过三分之一。

与我自始至终使用的 njuffa 相反double。在步骤 2 中，我选择了 5 阶的非理性方法（通过推广 Lagny 的方法获得；参见第一部分和附录 B），在步骤 4 中，我选择了 4 阶而不是 3 阶的有理方法。

所得方法的错误取整率为 6.10(25)/十亿，“潜在错误取整”路径的通过率为 3.05(18)⋅10 -7。

于 2021-07-17T14:57:07.160 回答