Pas*_*uoq 9 algorithm floating-point ieee-754
假设可以使用正确舍入的标准库函数,例如CRlibm中的函数.那么如何计算双精度输入的正确圆角立方根？
引用FAQ时,这个问题不是"[我]面临的实际问题".这有点像家庭作业.但立方根是经常发现的操作,可以想象这个问题是某人面临的实际问题.
由于"最好的Stack Overflow问题中有一些源代码",这里有一些源代码:
y = pow(x, 1. / 3.);
上面没有计算正确的圆形立方根,因为1/3不能完全表示为a double.
补充说明:
的文章介绍如何计算浮点立方根,但推荐牛顿迭代算法的最后一次迭代(S)将不得不做更高精度的算法来计算正确舍入的双精度立方根.这可能是计算它的最佳方式,但我仍在寻找一种利用现有正确舍入标准化功能的快捷方式.
C99包含一个cbrt()函数,但不能指望所有编译器都能正确舍入甚至忠实.CRlibm的设计者本可以选择包含cbrt()在提供的功能列表中,但他们没有.可以参考其他正确舍入的数学函数库中可用的实现.
恐怕我不知道如何保证正确舍入的双精度立方根，但可以提供一个非常接近正确舍入的双精度立方根，正如问题中也提到的那样。换句话说，最大误差似乎非常接近 0.5 ulp。
\n\nPeter Markstein，“IA-64 和基本函数：速度和精度”（Prentice-Hall 2000）
\n\n提出了基于 FMA 的有效技术，用于正确舍入倒数、平方根和倒数平方根，但在这方面它没有涵盖立方根。一般来说，Markstein 的方法需要在最终舍入序列之前获得精确到 1 ulp以内的初步结果。我没有数学资源来将他的技术扩展到立方根的舍入，但在我看来，这原则上应该是可能的，这是一个有点类似于倒数平方根的挑战。
\n\n按位算法可以轻松地计算具有正确舍入的根。由于 IEEE-754 舍入模式到最接近或什至的平局情况不会发生，因此只需进行计算，直到产生所有尾数位加上一位舍入位。基于二项式定理的平方根逐位算法在非恢复和恢复变体中都是众所周知的，并且一直是硬件实现的基础。通过二项式定理的相同方法适用于立方根，并且有一篇鲜为人知的论文列出了非恢复实现的细节：
\n\nH. Peng，\xe2\x80\x9c 提取平方根和立方根的算法，\xe2\x80\x9d 第五届 IEEE 国际计算机算术研讨会论文集，第 121-126 页，1981 年。
\n\n从实验中我可以看出，这对于从整数中提取立方根来说足够有效。由于每次迭代仅产生一个结果位，因此速度并不快。对于浮点运算中的应用程序，它的缺点是使用几个簿记变量，这些变量大约需要最终结果位数的两倍。这意味着需要使用 128 位整数运算来实现双精度立方根。
\n\n下面我的 C99 代码基于Halley 有理数立方根方法，该方法具有三次收敛性，这意味着初始近似值不必非常准确，因为每次迭代中有效数字的数量都会增加三倍。计算可以以多种方式安排。一般来说，将迭代方案安排为在数值上是有利的
\n\n新的猜测 := 旧的猜测 + 修正
\n\n因为对于足够接近的初始猜测，correction明显小于old_guess。这导致立方根的以下迭代方案：
x := x - x * (x 3 - a) / (2*x 3 + a)
\n\n卡汉关于立方根的注释中也列出了这种特殊的排列。它的另一个优点是自然地适合使用FMA（融合乘加）。一个缺点是 2*x 3的计算可能会导致溢出，因此至少部分双精度输入域需要一种参数缩减方案。在我的代码中，我只是将参数减少应用于所有非异常输入。
\n\n区间 [0.125,1) 用作主要近似区间。使用多项式极小极大近似来返回 [0.5,1] 中的初始猜测。窄范围有利于使用单精度算术来进行计算的低精度部分。
\n\n我无法证明我的实现的错误界限，但是，使用 2 亿个随机测试向量针对参考实现（精确到约 200 位）进行测试，发现总共 277 个不正确舍入的结果（因此错误率约为 1.4 ppm）最大误差为 0.500012 ulps。
\n\ndouble my_cbrt (double a)\n{\n double b, u, v, r;\n float bb, uu, vv;\n int e, f, s;\n\n if ((a == 0.0) || isinf(a) || isnan(a)) {\n /* handle special cases */\n r = a + a;\n } else {\n /* strip off sign-bit */\n b = fabs (a);\n /* compute exponent adjustments */\n b = frexp (b, &e);\n s = e - 3*342;\n f = s / 3;\n s = s - 3 * f;\n f = f + 342;\n /* map argument into the primary approximation interval [0.125,1) */\n b = ldexp (b, s);\n bb = (float)b;\n /* approximate cube root in [0.125,1) with relative error 5.22e-3 */\n uu = 0x1.2f32c0p-1f;\n uu = fmaf (uu, bb, -0x1.62cc2ap+0f);\n uu = fmaf (uu, bb, 0x1.7546e0p+0f);\n uu = fmaf (uu, bb, 0x1.5d0590p-2f);\n /* refine cube root using two Halley iterations w/ cubic convergence */\n vv = uu * uu;\n uu = fmaf (fmaf (vv, uu, -bb) / fmaf (vv, 2.0f*uu, bb), -uu, uu);\n u = (double)uu;\n v = u * u; // this product is exact\n r = fma (fma (v, u, -b) / fma (v, 2.0*u, b), -u, u);\n /* map back from primary approximation interval by jamming exponent */\n r = ldexp (r, f);\n /* restore sign bit */\n r = copysign (r, a);\n }\n return r;\n}\n小智 6
我最近编写了一个正确舍入的立方根，在https://github.com/mockingbirdnest/Principia/blob/master/numerics/cbrt.cpp中实现，并在https://github.com/mockingbirdnest/Principia/blob中记录/master/documentation/cbrt.pdf。\n该实现采用 C++ 语言并使用 Intel 内在函数，但应该很容易将其适应一个\xe2\x80\x99 选择的语言和库。
\n下面是实际计算的简要总结；这里的参考文献使用该文档中的书目代码，因为我显然声誉不佳，无法在此时链接许多内容。对 \xe2\x80\x9cpart I+\xe2\x80\x9d 或 \xe2\x80\x9cappendix [A-Z]\xe2\x80\x9d 的引用是 cbrt.pdf。
这可以通过使用几乎正确的忠实方法来完成，如果结果太接近平局，则使用一轮经典的逐位算法（类似于长除法的算法）来获得第54位；由于该算法计算向 0 舍入的结果的位，并且由于立方根不存在中间情况，因此该位足以确定舍入到最接近的值。
\n有可能获得接近正确的方法，该方法比大多数现有的立方根实现更快，并且足够正确，使得平均成本基本上不受校正步骤的影响；参见附录 D 和 F。
\n\xc2\xab\xc2\xa0 几乎正确的忠实方法\xc2\xa0\xc2\xbb 部分可能更有趣；如果没有 FMA，基本轮廓与 Kahan\xe2\x80\x99s [KB01] 中的相同：
\n在步骤2中使用巧妙的求根方法，可以在保持或提高性能的同时降低误舍率。\n在这种情况下我挖出了一种不合理的方法，
\n\n\n\xe2\x86\xa6 \xc2\xbd + \xe2\x88\x9a(\xc2\xbc\xc2\xb2 + /(3)) 近似 \xe2\x88\x9b(\xc2\xb3+)，
\n
摘自 Thomas Fantet de Lagny 所著的 17 世纪书籍，[Fan92]；也许令人惊讶的是，因为它同时具有除法和平方根，因此在适当重写时（以便尽早安排除法，并避免平方根和除法之间的串行依赖性，请参阅附录 D），其性能类似于更广为人知的有理方法（现在通常称为 Halley\xe2\x80\x99s，但这两种方法均由 Halley 考虑，并且在应用于立方根时均归因于 Lagny；请参阅第一部分中的讨论）。这是因为有理方法的除数较高，因此无法提前安排其除法。\n该无理方法的误差是有理方法的一半。
\n以Stephen Canon的推文（两个推特线程，[Can18a] 和 [Can18b]）启发的方式优化系数以最小化步骤 2 结束时的误差，在该非理性方法的误差上又获得两位成本。
\n通过步骤 4 中的 5 阶有理方法，该方法实现了每百万 4.564(68) 的误舍率，这很容易以基本上没有平均成本的方式进行纠正（缓慢的 \xe2\x80\x9c 潜在误舍入中的通过率\ xe2\x80\x9d 路径为 2.6480(52)\xe2\x8b\x8510 \xe2\x88\x924，慢速路径的延迟小于快速路径的十倍）。
\n关键思想在于njuffa \xe2\x80\x99s 对这个问题的回答：步骤 4 可以仅使用精确的平方来执行，而不需要精确的立方，因此步骤 2 的目标应该是，步骤 3 应该四舍五入到一半的精度而不是三分之一。
\n与我自始至终使用的 njuffa 相反double。在步骤 2 中，我选择了 5 阶无理方法（通过推广 Lagny\xe2\x80\x99s 方法获得；请参阅第一部分和附录 B），在步骤 4 中，选择 4 阶而不是 3 阶有理方法。
所得方法的误舍率为每十亿分之 6.10(25)，而 \xe2\x80\x9c 潜在误舍\xe2\x80\x9d 路径中的通过率为 3.05(18)\xe2\x8b\x8510 \xe2 \ x88\x927。
\n鉴于曲线x = y ^ 3上存在大量易于计算的有理点,我很想减小s ^ 3~x,s理性且只有几位宽.然后你有:
cbrt(x) = s * cbrt(1 + (x - s^3)/s)
显而易见的是使用您最喜欢的系列近似来评估校正项,并通过头尾FMA算法计算残差,以便在需要时将结果上调或下调一个ulp(大部分时间不需要完整计算) ,显然).
这不完全是问题的精神,但绝对可以使用,并且很容易以这种方式证明必要的界限.希望其他人可以提出更聪明的建议(我已经用尽了这个月的聪明才智).
| 归档时间： |
|
| 查看次数： |
721 次 |
| 最近记录： |