numpy浮点数的精度限制?

Question

给定序列

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

如何计算在哪个点上我将无法以精度E区分两个连续元素1/I和1/i + 1如果我使用numpy.float16？ 例如，"i"是什么?

其他np-float呢?

最小的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中可表示的最大i是什么?floatXX ?

---

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

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耗尽了时间。因为我很便宜，不想买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

正如你在评论中提到的，i是

>>> 2**12 - 6
4090

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