提问人：jmd_dk 提问时间：11/24/2022 更新时间：11/29/2022 访问量：306
精确三次输入的立方根中的浮点误差
Floating-point errors in cube root of exact cubic input
问：
我发现自己需要计算“整数立方根”，即整数的立方根，四舍五入到最接近的整数。在 Python 中，我们可以使用 NumPy 浮点函数：cbrt()
import numpy as np
def icbrt(x):
return int(np.cbrt(x))
虽然这在大多数时候都有效，但它在某些输入时失败了，结果是比预期的少一个。例如，，之所以出现，是因为 .下面查找前 100,000 个此类故障：xicbrt(15**3) == 14np.cbrt(15**3) == 14.999999999999998
print([x for x in range(100_000) if (icbrt(x) + 1)**3 == x])
# [3375, 19683, 27000, 50653] == [15**3, 27**3, 30**3, 37**3]
问题：，，，，...有什么特别之处，使回报率略低于确切结果？我找不到这些数字的明显潜在模式。15273037cbrt()
几点观察：
- 如果我们从 NumPy 切换到 Python 模块，或者如果我们从 Python 切换到 C，情况是一样的（这并不奇怪，因为我相信最终两者都是从 C 数学库委托给的）。
cbrt()mathnumpy.cbrt()math.cbrt()cbrt() - 在 C 中替换为 （） 会导致更多的失败情况。让我们坚持.
cbrt(x)x**(1/3)pow(x, 1./3.)cbrt() - 对于平方根，不会出现类似的问题，这意味着返回所有结果的正确结果（测试高达 100,000,000）。测试代码：
import numpy as np def isqrt(x): return int(np.sqrt(x))xprint([x for x in range(100_000) if (y := np.sqrt(x))**2 != x and (y + 1)**2 <= x])
额外
由于上述方法似乎仅在三次输入上失败，因此我们可以通过添加修复来纠正偶尔的错误，如下所示：icbrt()
import numpy as np
def icbrt(x):
y = int(np.cbrt(x))
if (y + 1)**3 == x:
y += 1
return y
另一种解决方案是坚持精确整数计算，在不使用浮点数的情况下实现。例如，在这个 SO 问题中对此进行了讨论。这种方法的另一个好处是，它们比使用浮点快（或可能快）。icbrt()cbrt()
需要明确的是，我的问题不是关于如何写出更好的，而是关于为什么在某些特定输入上失败。icbrt()cbrt()
答：
2赞
Eric Postpischil
11/24/2022
#1
此问题是由 的错误实现引起的。它不是由浮点运算引起的，因为浮点运算不是计算立方根的障碍，当完全正确的结果可以以浮点格式表示时，无法返回完全正确的结果。cbrt
例如，如果使用整数算术来计算 80 的五分之九，我们期望正确结果为 144。如果计算一个数字的五分之九的例程被实现为 ，我们会责怪该例程执行不正确，而不是责怪整数算术没有处理分数。同样，如果一个例程使用浮点运算来计算不正确的结果，而正确的结果是可表示的，我们责怪例程，而不是浮点运算。int NineFifths(int x) { return 9/5*x; }
有些数学函数很难计算，我们接受其中的一些误差。事实上，对于数学库中的一些例程，人类还没有弄清楚如何在已知的有界执行时间内用正确的舍入来计算它们。因此，我们接受并非每个数学例程都是正确舍入的。
然而，当函数的数学值可以用浮点格式精确表示时，可以通过忠实的舍入而不是正确的舍入来获得正确的结果。因此，这是数学库函数的理想目标。
正确四舍五入意味着计算结果等于通过将精确的数学结果四舍五入到最接近的可表示值而获得的数字。1 四舍五入表示计算结果与精确的数学结果相差不到 1 个 ULP。ULP 是最低精度的单位，即两个相邻可表示数字之间的距离。
正确地舍入函数可能很困难，因为通常，函数可以任意接近舍入决策点。对于四舍五入到最接近的数字，这是两个相邻的可表示数字之间的中间位置。考虑两个相邻的可表示数字 a 和 b。它们的中点是 m = （a+b）/2。如果某个函数 f（x） 的数学值刚好低于 m，则应将其四舍五入为 a。如果它刚好在上面，则应四舍五入为 b。当我们在软件中实现 f 时，我们可能会用一些非常小的误差 e 来计算它。当我们计算 f（x） 时，如果我们的计算结果位于 [m-e， m+e]，并且我们只知道误差界是 e，那么我们就无法判断 f（x） 是低于 m 还是高于 m。而且，一般来说，函数 f（x） 可以任意接近 m，这总是一个问题：无论我们计算 f 多么精确，无论我们使误差界 e 多么小，我们的计算值都有可能非常接近中点 m，比 e 更近，因此我们的计算不会告诉我们是向下舍入还是向上舍入。
对于某些特定函数和浮点格式，已经进行了研究并撰写了关于函数接近此类舍入决策点的接近程度的证明，因此某些函数（如正弦和余弦）可以在计算时间的已知边界下通过正确的舍入来实现。到目前为止，其他功能尚未得到证明。
相比之下，忠实舍入更容易实现。如果我们计算一个误差范围小于 1/2 ULP 的函数，那么我们总是可以返回一个忠实的舍入结果，该结果与精确数学结果在一个 ULP 范围内。一旦我们计算了一些结果 y，我们将其四舍五入到最接近的可表示值2 并返回它。从 y 的误差小于 1/2 ULP 开始，四舍五入可能会增加 1/2 ULP 的误差，因此总误差小于 1 个 ULP，即忠实的四舍五入。
忠实舍入的一个好处是，当确切的结果可表示时，函数的忠实舍入实现总是会产生确切的结果。这是因为下一个最接近的结果距离一个 ULP，但忠实舍入的误差始终小于 1 个 ULP。因此，忠实舍入的函数在可表示时返回精确的结果。cbrt
15、27、30、37 ...，使回报率略低于确切结果，有什么特别之处？我找不到这些数字的明显潜在模式。
cbrt()
错误的实现可能会通过将参数减少到 [1， 8） 或类似区间中的值，然后应用预先计算的多项式近似值来计算多维数据集根。该多项式中的每次加法和乘法都可能引入舍入误差，因为每次运算的结果都以浮点格式四舍五入到最接近的可表示值。此外，多项式具有固有误差。舍入误差的行为有点像随机过程，有时向上舍入，有时向下舍入。当它们在多次计算中累积时，它们可能会碰巧向不同的方向四舍五入并取消，或者它们可能会向同一方向四舍五入并加强。如果错误恰好在计算结束时取消，则从 中得到一个确切的结果。否则，您可能会从 获得不正确的结果。cbrtcbrtcbrt
脚注
1 In general, there is a choice of rounding rules. The default and most common is round-to-nearest, ties-to-even. Others include round-upward, round-downward, and round-toward-zero. This answer focuses on round-to-nearest.
2 Inside a mathematical function, numbers may be computed using extended precision, so we may have computed results that are not representable in the destination floating-point format; they will have more precision.
评论
(2*3*4*5*6*7*8*9)/7.0 == 51840.0sqrt(123456789.0**2) == 123456789.0cbrtn**(1/3)1/3cbrt1/3cbrtcbrt