浮点倒数总是往返吗？

Question

对于 IEEE-754 算法，是否可以保证倒数的最后一位精度为 0 或 1 个单位？据此，倒数的倒数是否有保证的误差范围？

Answer 1

[以下所有内容均采用固定的 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 缩放为 2 的幂，使其位于 [1.0, 2.0) 范围内。在 x 是 2 的幂的情况下，上面的结果显然是正确的，所以我们假设它不是（这将在下面的第二步中有用）。下面的证明是针对特定于 IEEE 754 binary64 格式的精度给出的，但它直接适用于任何 IEEE 754 二进制格式。

为倒数的true值写1 / x，在四舍五入之前，让y成为最近的（唯一的，事实证明）可表示的 binary64 浮点数到 1 / x。那么：

因为 y 是最接近 1 / x 的浮点数，并且 y 和 1/x 都在 binade [0.5, 1.0] 中，其中连续浮点数之间的间距恰好是 2^-53，我们有 |y - 1/x| <= 2^-54。其实我们可以做得更好一点：
我们上面实际上有一个严格的不等式：|y - 1/x| < 2^-54。如果 |y - 1/x| 正好等于 2^-54，那么 1/x 可以用任意精度的二进制浮点数精确表示（因为 y 和 2^-54 是）。但是，1/x 可以精确表示的唯一二进制浮点数 x 是 2 的幂，我们已经排除了这种情况。
if x < sqrt(2) then 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 并且往返成功。这就完成了第2部分的证明。
在一般情况下x < 2，我们有x / y < 4，所以|x - 1/y| < 2^-52。这使得 1/y 最多与 x 相差 1 ulp，从而完成了第 1 部分的证明。

下面是 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

并演示最大误差不超过1个ulp，一般来说（对于binary64格式，在binade [1.0, 2.0]中，最后一位的1个单位是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 binary64 格式的例子，表明避免下溢的限制是必要的：

>>> 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 也涵盖了十进制浮点类型，我应该提到上面的证明几乎逐字地延续到十进制的情况，确定有效数小于 [=73= 的浮点数]，往返发生，而对于一般的十进制浮点数（再次避免上溢和下溢），我们在最后一个地方永远不会超过一个单位。然而，需要一些数论技巧来证明关键不等式 |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')

浮点倒数总是往返吗？

Does a floating-point reciprocal always round-trip?

floating-point

precision

floating-accuracy

ieee-754