精确三次输入的立方根中的浮点错误

Question

我发现自己需要计算；整数立方根"；，意思是整数的立方根，四舍五入到最接近的整数。在Python中，我们可以使用NumPy浮点cbrt()函数：

import numpy as np
def icbrt(x):
return int(np.cbrt(x))

尽管这在大多数情况下都有效，但它在特定的输入x处失败，结果比预期少一个。例如，icbrt(15**3) == 14，这是因为np.cbrt(15**3) == 14.999999999999998。以下发现了前100000个此类故障：

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]

问题：15、27、30、37…有什么特别之处。。。，使得CCD_ 9的返回比精确结果稍微低一点？我找不到这些数字的明显潜在模式。

一些观察：

• 如果我们从NumPy的cbrt()切换到Python的math模块，或者如果我们从Python切换到C，情况是一样的(这并不奇怪，因为我相信numpy.cbrt()和math.cbrt()最终都从C数学库委托给cbrt())
• 用x**(1/3)(C中的pow(x, 1./3.))代替cbrt(x)会导致更多的故障情况。让我们坚持cbrt()
• 对于平方根，不会出现类似的问题，这意味着

  import numpy as np
  def isqrt(x):
  return int(np.sqrt(x))
  

  返回所有CCD_ 19(最多测试100000000)的正确结果。测试代码：

  print([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

另一种解决方案是坚持精确整数计算，在不使用浮点数的情况下实现icbrt()。例如，在本SO问题中对此进行了讨论。这种方法的一个额外好处是，它们比使用浮点cbrt()更快。

需要明确的是，我的问题不是如何编写更好的icbrt()，而是为什么cbrt()在某些特定输入上失败。

Answer 1

此问题是由cbrt的错误实现引起的。它不是由浮点运算引起的，因为浮点运算并不是很好地计算立方根的障碍，当完全正确的结果可以用浮点格式表示时，它不能返回完全正确的结论。

例如，如果使用整数运算来计算80的五分之九，我们将期望144的正确结果。如果一个计算数字五分之九的例程被实现为int NineFifths(int x) { return 9/5*x; }，我们会责怪该例程实现错误，而不是责怪整数运算没有处理分数。类似地，如果一个例程使用浮点运算来计算不正确的结果，而正确的结果是可表示的，我们会责怪例程，而不是浮点运算。

有些数学函数很难计算，我们接受其中的一些误差。事实上，对于数学库中的一些例程，人类还没有弄清楚如何在已知的有限执行时间内通过正确的舍入来计算它们。因此，我们承认并不是每个数学例程都是正确的四舍五入。

然而，当函数的数学值可以精确地用浮点格式表示时，可以通过忠实的舍入而不是正确的舍入来获得正确的结果。因此，这是数学库函数的理想目标。

正确取整意味着计算结果等于通过将精确的数学结果取整到最接近的可表示值而获得的数字¹完全四舍五入意味着计算结果与精确的数学结果相比小于一个ULP。ULP是最小精度的单位，即两个相邻可表示数字之间的距离。

正确舍入函数可能很困难，因为通常情况下，函数可以任意接近舍入决策点。对于四舍五入，这是两个相邻可表示数字之间的中间值。考虑两个相邻的可表示数a和b。它们的中点是m=(a+b)/2。如果某个函数f(x)的数学值刚好低于m，它应该四舍五入到a。如果它刚好高于m，它就应该四舍五入到b。当我们在软件中实现f时，我们可能会用一些非常小的误差e来计算它。当我们计算f(x，一般来说，函数f(x)可以任意接近m，这总是一个问题：无论我们计算f有多准确，无论我们使误差界e有多小，我们的计算值都有可能非常接近中点m，比e更接近，因此我们的计算不会告诉我们是向下取整还是向上取整。

对于一些特定的函数和浮点格式，已经进行了研究并证明了函数与这种舍入决策点的接近程度，因此某些函数(如正弦和余弦)可以在计算时间上以已知的边界通过正确的舍入来实现。到目前为止，其他功能还没有得到证实。

相比之下，忠实取整更容易实现。如果我们计算一个误差界小于½ULP的函数，那么我们总是可以返回一个忠实的四舍五入结果，一个在精确数学结果的一个ULP内的结果。一旦我们计算了某个结果y，我们就将其四舍五入到最接近的可表示值²，并返回该值。从y的误差小于½ULP开始，四舍五入可能会增加½ULP的误差，因此总误差小于一个ULP，这是忠实的四舍五进。

忠实取整的一个好处是，当精确的结果是可表示的时，函数的忠实取整实现总是产生精确的结果。这是因为下一个最接近的结果是一个ULP，但忠实舍入的误差总是小于一个ULP。因此，当它们是可表示的时，忠实取整的cbrt函数返回精确的结果。

15、27、30、37…有什么特别之处。。。，使得CCD_ 28的返回值比精确结果略低？我找不到这些数字的明显潜在模式。

坏的cbrt实现可能通过将参数减少为[1，8)或类似区间，然后应用预先计算的多项式近似。该多项式中的每次加法和乘法都可能引入舍入误差，因为每次运算的结果都被舍入到浮点格式中最接近的可表示值。此外，多项式具有固有误差。舍入误差的行为有点像随机过程，有时会向上舍入，有时会下降。当它们在几次计算中累积时，它们可能会碰巧向不同的方向取整并取消，也可能向同一方向取整以加强。如果在计算结束时错误恰好取消，则可以从cbrt中获得准确的结果。否则，您可能会从cbrt中得到不正确的结果。

脚注

¹一般来说，可以选择舍入规则。默认和最常见的是四舍五入，平局为偶数。其他包括向上取整、向下取整和向零取整。这个答案侧重于四舍五入。

²在数学函数中，可以使用扩展精度来计算数字，因此我们可能会得到无法用目标浮点格式表示的计算结果；它们将具有更高的精度。