精确立方数的立方根计算为何出现浮点误差？

背景 #

计算整数立方根（对立方根向下取整）时，用Python的NumPy或math模块的cbrt()函数实现会出现特定输入下的错误：

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

多数情况该方法有效，但部分立方数输入会得到比预期小1的结果，比如icbrt(15**3)返回14，原因是np.cbrt(15**3)的计算结果是14.999999999999998。在100000以内的失败案例为[3375, 19683, 27000, 50653]，对应15³、27³、30³、37³。

核心问题 #

15、27、30、37这些数字有何特殊之处，导致cbrt()返回略低于精确值的结果？

相关观察 #

• 切换为Python math模块的cbrt()或C语言标准库的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))

该函数对所有x（测试至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

也可采用纯整数计算实现整数立方根，无需浮点运算，这类方法通常比浮点cbrt()更快。

问题解答：为什么这些数会出错？ #

这本质是IEEE 754双精度浮点数的表示局限性结合cbrt()的数值计算误差导致的：

1. 二进制浮点数的精度限制：计算机用双精度浮点数（64位）存储数值，其尾数只有52位（加隐含位共53位），只能精确表示有限的整数和分数。对于15、27、30、37这些数，它们的立方是整数，但对应的立方根无法用二进制浮点数精确表示——真实的整数立方根在二进制中是无限循环的分数，cbrt()只能返回最接近的可表示值，而这个值恰好略小于真实值。

2. cbrt()的数值算法偏移：C标准库的cbrt()通常用迭代近似算法（如牛顿-拉夫逊法）实现，在处理这些特定数时，迭代的收敛结果会向负方向产生微小偏移。结合浮点数的精度限制，最终输出的结果会比真实整数小一个极小的epsilon（例如14.999999999999998），经过int()截断后就得到了比预期小1的结果。

3. 平方根无类似问题的原因：对于完全平方数，其平方根要么是整数（可被双精度浮点数精确表示），要么是有限二进制分数；同时sqrt()函数的实现针对完全平方数做了优化，能确保返回精确的浮点值。双精度浮点数的53位尾数足够覆盖极大范围的整数平方根，因此不会出现这类截断错误。

内容的提问来源于stack exchange，提问作者jmd_dk