为何x=10^9时，1/[(x²-1)^0.5 -x]会触发ZeroDivisionError？

问题描述 #

练习中计算x=10^{9时，表达式出现异常，但改用`np.float128`类型却能正常运行。原以为是存储范围不足，但64位有符号整数最大值约9×10}18（大于x²=10^{18），且64位浮点数的绝对值范围覆盖10}-308到10^308，范围足够，因此疑惑问题根源。

失败代码（Python原生float） #

x = 10 ** 9
f00 = 1 / ((x ** 2 - 1) ** 0.5 - x)
print(f00)
f01 = -(x ** 2 - 1) ** 0.5 - x
print(f01)

成功代码（np.float128） #

import numpy as np

x = np.float128(10 ** 9)
f00 = 1 / ((x ** 2 - 1) ** 0.5 - x)
print(f00)
f01 = np.float128(-(x ** 2 - 1) ** 0.5 - x)
print(f01)

问题根源：浮点数精度丢失，而非范围不足 #

这是双精度浮点数的精度限制导致的，和存储范围无关：

• Python原生float是64位双精度浮点数，仅保留约15-17位有效数字。当x=10^{9时，x²=10}18，此时双精度浮点数的最小可分辨间隔（ULP，单位最后位置）为256。x²-1与x²的差值为1，远小于256，因此双精度会将x²-1和x²视为同一个值。
• 这会导致(x**2 -1)**0.5计算结果等于x，分母变为x - x = 0，触发除以零错误。
• 而np.float128是128位扩展精度浮点数，有效位数约33-34位，足以分辨x²和x²-1的差异，因此能正确计算出平方根，分母不为零，运行正常。

另外，x=10^7时x²=1014，此时双精度浮点数的ULP约0.022，x²-1与x²的差值1大于ULP，因此能正常区分，不会出现异常。

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