一个浮点数的倒数是否总是能精确还原？

Question

一个浮点数的倒数是否总是能精确还原？

floating-pointprecisionfloating-accuracyieee-754

9

对于IEEE-754算术，是否保证倒数的最后一位精度为0或1？从而，倒数的倒数是否存在保证的误差界限？

- Raymond Hettinger

2

这个问题与实现单精度乘以常数的问题有些相关，方法是通过双精度乘以倒数再舍入到单精度。这个问题也被Mark回答过：https://dev59.com/vXjZa4cB1Zd3GeqPdWSR - Pascal Cuoq

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Mark Dickinson · Accepted Answer

[假设使用IEEE 754二进制格式，采用四舍五入作为取整模式。]

由于倒数（计算为1/x）是基本的算术操作，1 可以被精确地表示，并且标准保证了算术运算的正确舍入，因此倒数结果的真实值与计算结果之间的差距在最后一位小数点的范围内不超过 0.5 个单位。（这适用于标准中指定的任何基本算术操作。）

一般而言，一个值x的倒数的倒数不一定等于x。下面是一个在IEEE 754 binary64格式下的快速示例：

>>> x = 1.8
>>> 1.0 / (1.0 / x)
1.7999999999999998
>>> x == 1.0 / (1.0 / x)
False

假设溢出和下溢已经避免，并且 x 是有限的、非零的、并且能够准确的表示，那么以下结论是成立的:

1.0 / (1.0 / x) 的值与 x 之间的差距不会超过最后一位的1。
如果 x 的有效数字(通常归一化到范围 [1.0,2.0))小于 sqrt(2)，那么倒数也是这样的：1.0 / (1.0 / x) == x。

证明概述: 不失一般性，我们可以假设 x 是正的，并将 x 缩放为二的幂以使它位于范围 [1.0,2.0)。显然，在 x 是一个二的幂的情况下，上述结果是正确的，所以我们假设它不是（在下面的第二步中这将是有用的）。下面的证明是针对 IEEE 754 二进制64位格式的精度给出的，但它直接适用于任何IEEE 754二进制格式。写 1/x 表示倒数的真实值，在舍入之前，让 y 是最接近的可表示的二进制64位浮点数到 1/x。则：

由于 y 是最接近 1/x 的浮点数，并且两者都在间隔为 2^-53 的区间 [0.5,1.0] 中，因此我们有 |y-1/x| <= 2^-54。实际上，我们可以做得更好:
实际上，我们上面给出的是一个严格的不等式: |y-1/x| < 2^-54。如果 |y-1/x| 恰好等于 2^-54，那么 1/x 将在任意精度的二进制浮点数中被精确表示（因为 y 和 2^-54 都是可表示的）。但是，只有对于二的幂，才有二进制浮点数 x 使得 1/x 在某些精度下被精确地表示，而我们已经排除了这种情况。
如果 x < sqrt(2)，则 1 / x > x / 2，因此（将两者四舍五入为可表示的浮点数），我们有 y >= x / 2，所以 x / y <= 2。
现在 x - 1/y = (y - 1/x) x/y，从对 |y - 1/x| 和 x/y 的限制中（仍假设 x < sqrt(2)），我们得到 |x - 1/y| < 2^-53。由此可知，x 是最接近 1/y 的可表示浮点数，1/y 四舍五入为 x，并且往返转换成功。这完成了第二部分的证明。
在一般情况下，x < 2，我们有 x / y < 4，因此 |x - 1/y| < 2^-52。这使得 1/y 距离 x 最多只有 1 ulp，这完成了第一部分的证明。
以下是对 sqrt(2) 阈值的演示：使用 Python，在范围为 [1.0，2.0) 的一百万个随机浮点数中，找出那些不能通过倒数回路的浮点数。所有小于 sqrt(2) 的样本都可以进行往返转换。
```
>>> import random
>>> samples = [random.uniform(1.0, 2.0) for _ in range(10**6)]
>>> bad = [x for x in samples if 1.0 / (1.0 / x) != x]
>>> len(bad)
171279
>>> min(bad)
1.4150519879892107
>>> import math
>>> math.sqrt(2)
1.4142135623730951
```
在评论中，@Kelly Bundy 发现 IEEE 754 binary64 值 x 最接近 √2 的值，使得 1.0 / (1.0 / x) != x，其值为 0x1.6a09e6964b6acp+0，约为 1.4142135731630487，与 √2 相差大约 48.6 百万个 ulps。

同时证明了误差最大一般不超过 1 ulp（对于二进制64位浮点数，在区间 [1.0, 2.0) 中，最后一位的单位是 2^-52）：
```
>>> samples = [random.uniform(1.0, 2.0) for _ in range(10**6)]
>>> max(abs(x - 1.0 / (1.0 / x)) for x in samples)
2.220446049250313e-16
>>> 2**-52
2.220446049250313e-16
```
这是一个使用IEEE 754二进制64格式的示例，说明避免下溢限制是必要的：
```
>>> x = 1.3e308
>>> x_roundtrip = 1.0 / (1.0 / x)
>>> x.hex()
'0x1.72409614c1e6ap+1023'
>>> x_roundtrip.hex()
'0x1.72409614c1e6cp+1023'
```
这里的x_roundtrip结果与原始结果相差两个单位的最后一位，因为1/x小于最小可表示浮点数，因此没有像x一样以相同的精度表示。
最后说明：由于IEEE 754-2008也涵盖了十进制浮点类型，因此我应该提到上述证明可以几乎毫无改动地延续到十进制情况，证明了对于小于sqrt(10)的有效数字的浮点数，会发生往返旅行，而对于一般的十进制浮点数（再次避免溢出和下溢），我们永远不会超过最后一位误差一个单位。然而，需要一些数论技巧来证明关键不等式|x-1/y| < 1/2 10^(1-p)总是严格的：最终需要证明量1+16 10^(2p-1)从未是平方数（这是真的，但在本网站范围内包括此处的证明可能有些超出范围）。
```
>>> from decimal import Decimal, getcontext
>>> import random
>>> getcontext().prec = 6
>>> samples = [+Decimal(random.uniform(1.0, 10.0)) for _ in range(10**6)]
>>> bad = [x for x in samples if 1 / (1 / x) != x]
>>> min(bad)
Decimal('3.16782')
```