numpy浮点数的精度限制？

Question

给定顺序

1/1, 1/2, 1/3, ... , 1/n

如果我使用 1/i 和 1/i+1，我如何计算在什么时候我将无法以精确度 E 区分两个连续元素numpy.float16 ？即“我”是什么？

其他 np-floats 呢？

最小的 E 是多少？以及如何计算它的'i'？

例如，如果 E = 0.01，我可以区分 1/9 和 1/10，但不能区分 1/10 和 1/11，因为：

1/9 = 0.111
1/10 = 0.100
1/11 = 0.091

0.111 - 0.100 = 0.01  >= E
0.100 - 0.091 = 0.009 < E

i = 10

---

以更抽象的方式，给定 f(i) 在 np.floatXX 中可表示的最大“i”是多少？

---

有趣的是，实践中的精度比计算的要差： /逻辑中断的地方/

for i in range(int(1e3),int(12e6)) : 
   if not np.floatXX(1/i) > np.floatXX(1/(i+1)) :
       print(i); break

float32: 11864338
float16: 1464

Answer 1

不要加 1，而是加倍分母。您可以放心地假设它是一些二进制数。这是一种简单的方法：

one = np.float64(1.0)
two = np.float64(2.0)
n = one
bits = 0
while one + n != one:
    bits += 1
    n /= two

您从bits = 0 开始，否则您将获得使您超过分辨率的位数。

最后，您会得到bits = 53，这是 IEEE-754 编码的 64 位浮点数中的位数。

这意味着对于以有效二进制科学记数法编码的任何数字，ULP（最低精度单位）大约为n * 2**-53。具体来说，n 是四舍五入到最高位的数字。您将无法解决浮点数中较小的相对变化。

奖励：为其他浮点类型运行上述代码：

float16 (half):   11 bits
float32 (single): 24 bits
float64 (double): 53 bits
float96 (sometimes longdouble): 80 bits
float128 (when available): 113 bits

---

您可以修改上面的代码以适用于任何目标号码：

target = np.float16(0.0004883)
one = np.float16(1.0)
two = np.float16(2.0)
n = two**(np.floor(np.log2(target)) - one)
bits = 0
while target + n != target:
    bits += 1
    n /= two

结果 (ULP) 由 n * 2 给出，因为在您失去分辨率后循环停止。这与我们以bits = 0 开头的原因相同。在这种情况下：

>>> n * two
5e-07

如果您预先知道尾数中的位数，您可以完全缩短计算。所以对于float16，其中bits = 11，你可以这样做

>>> two**(np.floor(np.log2(target)) - np.float16(bits))
5e-07

---

在这里阅读更多：

• https://en.wikipedia.org/wiki/IEEE_754
• https://en.wikipedia.org/wiki/Unit_in_the_last_place

Answer 2

我的另一个答案提供了您实际要求的理论，但需要一些重要的解释。这是缺少的步骤：

给定一个整数i，你可以写

1 / i - 1 / (i + 1) =
(i + 1 - i) / (i * (i + 1)) =
1 / (i * (i + 1)) =
1 / (i**2 + i)

要找到i 使得1 / (i**2 + i) 在某些二进制表示中低于1 / i 的ULP，您可以直接使用我的其他答案。

1 / i 的 ULP 由下式给出

ulp = 2**(floor(log2(1 / i)) - (bits + 1))

您可以尝试找到i 这样

1 / (i**2 + i) < 2**(floor(log2(1 / i)) - (bits + 1))
1 / (i**2 + i) < 2**floor(log2(1 / i)) / 2**(bits + 1)
2**(bits + 1) < (i**2 + i) * 2**floor(log2(1 / i))

由于floor 操作和Wolfram Alpha runs out of time，编写起来并不简单。既然我便宜又不想买 Mathematica，我们就做个大概吧：

2**(bits + 1) < (i**2 + i) * 2**floor(log2(1 / i))
2**(bits + 1) < (i**2 + i) / i
2**(bits + 1) < i + 1

您可能相差一个左右，但您应该看到在i = 2**(bits + 1) - 1 附近，差异不再是可解决的。事实上，对于float16 的 11 位尾数，我们看到：

>>> np.float16(1 / (2**12 - 1)) - np.float16(1 / (2**12))
0.0

这里的实际数字要少一点（请记住我们拿走floor 的近似值）：

>>> np.float16(1 / (2**12 - 5)) - np.float16(1 / (2**12 - 4))
0.0
>>> np.float16(1 / (2**12 - 6)) - np.float16(1 / (2**12 - 5))
2.4e-07

正如您在 cmets 中指出的，i 是

>>> 2**12 - 6
4090

您可以以类似的方式计算所有其他浮点类型的精确值。但这确实留给读者作为练习。