假设可以正确舍入标准库函数,例如CRlibm中提供的函数。那么,如何计算 double 输入的正确舍入的三次根呢?
引用常见问题解答,此问题不是“[I]面临的实际问题”。这样有点像功课。但是立方根是一种经常发现的操作,并且可以想象这个问题是某人面临的实际问题。
由于“最佳堆栈溢出问题中包含一些源代码”,因此这里有一些源代码:
y = pow(x, 1. / 3.);
上面的方法无法计算出正确四舍五入的立方根,因为1/3不能完全表示为
double
。补充笔记:
article描述了如何计算浮点立方根,但是对于该算法而言,为计算出正确舍入的 double 立方根,必须以更高的精度完成推荐的Newton-Raphson算法的最后一次迭代。那可能是计算它的最佳方法,但是我仍在寻找一种捷径,该捷径可以利用现有的正确四舍五入的标准化函数。
C99包含
cbrt()
函数,但是不能期望所有编译器都将其正确舍入为or even faithful。 CRlibm的设计师可以选择在提供的功能列表中包括cbrt()
,但他们没有。欢迎引用其他舍入为正确的数学函数的库中提供的实现。 最佳答案
恐怕我不知道如何保证正确舍入的 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 *(x3-a)/(2 * x3 + a)
Kahan's notes on cube root中也列出了这种特殊的安排。它具有进一步的优势,可以自然地将自身借给FMA (fused-multiply add)的使用。一个缺点是2 * x3的计算可能导致溢出,因此对于 double 输入域的至少一部分需要使用自变量约简方案。在我的代码中,我基于对IEEE-754 double 操作数指数的直接操作,将参数归约简单地应用于所有非异常输入。
间隔[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;
}
关于algorithm - 计算正确舍入/几乎正确舍入的浮点立方根,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/18063755/