假定可以正确舍入标准库函数，例如CRlibm中的函数。那么，如何计算双精度输入的正确舍入的三次根？
这个问题不是“ [我]面临的实际问题”，引用FAQ。这样有点像功课。但是立方根是一种经常发现的操作，并且可以想象这个问题是某人面临的实际问题。
因为“最佳堆栈溢出问题中包含一些源代码，所以这里有一些源代码：
y = pow(x, 1. / 3.);
上面没有计算出正确四舍五入的立方根，因为1/3不能精确表示为double。
附加说明：
article描述了如何计算浮点立方根，但是推荐的Newton-Raphson算法的最后一次迭代必须以更高的精度完成，才能计算出正确舍入的双精度立方根。这可能是计算它的最佳方法，但是我仍在寻找一种捷径，该捷径可以利用现有的正确四舍五入的标准化函数。
C99包括cbrt()函数，但是对于所有编译器，都不能期望它正确地四舍五入or even faithful。 CRlibm的设计人员可以选择将cbrt()包含在提供的功能列表中，但实际上并没有。欢迎引用其他舍入为正确的数学函数的库中提供的实现。
恐怕我不知道如何保证正确舍入的双精度多维数据集根，但是可以提供在问题中也提到过几乎正确舍入的一个。换句话说，最大误差似乎非常接近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's rational method for the cube root，它具有三次收敛性，这意味着初始逼近不必非常精确，因为有效位数在每次迭代中都增加了三倍。可以以各种方式来安排计算。通常，在以下情况下安排迭代方案在数值上是有利的：
new_guess：= old_guess +更正
由于有足够接近的初始猜测，correction明显小于old_guess。这导致多维数据集根的以下迭代方案：
x：= x-x *（x 3
-a）/（2 * x 3 + y）Kahan's notes on cube root中也列出了这种特殊的安排。它具有进一步自然使用FMA (fused-multiply add)的优势。一个缺点是2 * x 3
的计算可能会导致溢出，因此至少对双精度输入域的一部分需要使用自变量简化方案。在我的代码中，我基于对IEEE-754双精度操作数指数的直接操作，将参数归约简单地应用于all非异常输入。间隔[0.125,1）用作一次近似间隔。使用多项式最小极大逼近，可返回[0.5,1]中的初始猜测。窄范围便于对计算的低精度部分使用单精度算法。
我无法证明我的实现的错误范围，但是，使用2亿个随机测试向量对参考实现（准确度约为200位）进行测试，发现总共有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;
}
[鉴于曲线x = y ^ 3上有许多容易计算的有理点，我很想将s ^ 3〜x减小约s ^ 3〜x，且有理且只有几位宽。然后您有：