Ray*_*ger 4 floating-point precision floating-accuracy ieee-754
对于IEEE-754算术,倒数的最后位置是否保证0或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).上面的结果在2的x精确幂的情况下显然是正确的,所以我们假设它不是(这将在下面的第二步中有用).下面的证明是针对IEEE 754 binary64格式的特定精度给出的,但它直接适应任何IEEE 754二进制格式.
在舍入之前写出倒数1 / x的真值,并且y将(最接近的)表示为最接近的可表示的binary64 float to 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是).但唯一的二进制彩车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并且往返成功.这样就完成了第2部分的证明.
在一般情况下x < 2,我们有x / y < 4,所以|x - 1/y| < 2^-52.这使得1/y距离最多1 ulp x,这完成了第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二进制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')