解释区间 [0, 1] 中明显联系的舍入方向上的惊人奇偶性

Mar*_*son 8 python floating-point rounding ieee-754

考虑和0.xx5之间形式的浮点数集合:0.01.0[0.005, 0.015, 0.025, 0.035, ..., 0.985, 0.995]

我可以在 Python 中轻松列出所有 100 个这样的数字:

>>> values = [n/1000 for n in range(5, 1000, 10)]
Run Code Online (Sandbox Code Playgroud)

让我们看看前几个和最后几个值来检查我们没有犯任何错误:

>>> values[:8]
[0.005, 0.015, 0.025, 0.035, 0.045, 0.055, 0.065, 0.075]
>>> values[-8:]
[0.925, 0.935, 0.945, 0.955, 0.965, 0.975, 0.985, 0.995]
Run Code Online (Sandbox Code Playgroud)

现在我想将这些数字中的每一个四舍五入到点后两位小数。一些数字将被四舍五入;有些会被四舍五入。我有兴趣准确计算有多少舍入。我也可以在 Python 中轻松计算:

>>> sum(round(value, 2) > value for value in values)
50
Run Code Online (Sandbox Code Playgroud)

所以事实证明,这 100 个数字中正好有一半被四舍五入了。

如果你不知道 Python 在底层使用二进制浮点,这个结果就不足为奇了。毕竟,Python 的文档明确指出该round函数使用舍入到偶数(又名银行家舍入)作为其舍入模式,因此您希望这些值交替向上舍入和向下舍入。

但是Python引擎盖下使用二进制浮点,这意味着有例外的少数(即0.1250.3750.6250.875),这些值是准确的关系,而仅仅是非常好的二进制近似那些关系。毫不奇怪,对舍入结果的仔细检查表明这些值不会交替上下舍入。相反,每个值向上或向下舍入取决于二进制近似值发生在十进制值的哪一侧。所以没有先验的理由期望值的一半向上取整,一半取下。这让我们得到的结果正好是 50 有点令人惊讶。

但也许我们只是走运了?毕竟,如果你折腾一个公平的硬币100次,恰好有50头是不是不寻常的结果:它会以8%左右的概率发生。但事实证明,这种模式在小数位数较多的情况下仍然存在。以下是四舍五入到小数点后 6 位的类似示例:

>>> values = [n/10**7 for n in range(5, 10**7, 10)]
>>> sum(round(value, 6) > value for value in values)
500000
Run Code Online (Sandbox Code Playgroud)

这里再次将明显的联系四舍五入到点后的小数点后 8 位:

>>> values = [n/10**9 for n in range(5, 10**9, 10)]
>>> sum(round(value, 8) > value for value in values)
50000000
Run Code Online (Sandbox Code Playgroud)

所以问题是:为什么正好有一半的案例四舍五入?或者换句话说,为什么在这些小数关系的所有二进制近似值中,大于真实值的近似值数量与较小的近似值数量完全匹配?(可以很容易地证明,对于精确的情况,我们有相同数量的向上舍入和向下舍入,因此我们可以忽略这些情况。)

笔记

  1. 我假设 Python 3。
  2. 在典型的台式机或笔记本电脑上,Python 的浮点数将使用 IEEE 754 binary64(“双精度”)浮点格式,整数的真正除法和round函数都将是正确的舍入运算,使用舍入关系到- 甚至舍入模式。虽然语言本身不能保证这一切,但这种行为非常普遍,我们假设在这个问题中使用了这样一台典型的机器。
  3. 这个问题的灵感来自 Python 错误报告:https : //bugs.python.org/issue41198

Mar*_*son 4

事实证明,我们可以证明一些更强大的东西,这与十进制表示或小数舍入没有特别关系。这是更强有力的声明:

定理。选择一个正整数,并考虑由分数, , , ...,组成的n <= 2^1021长度序列。使用 IEEE 754 舍入方向将每个分数转换为最接近的 IEEE 754 二进制 64 浮点值。那么转换后的值大于原始分数的分数的数量将恰好等于转换后的值小于原始分数的分数的数量。n1/2n3/2n5/2n(2n-1)/2nroundTiesToEven

涉及浮点序列的原始观察结果[0.005, 0.015, ..., 0.995]是从上述语句的情况得出的n = 100:在 100 种情况中的 96 种情况中, 的结果round(value, 2)取决于舍入为 Binary64 格式时引入的错误的符号,并且根据上述语句,48 的这些情况将产生正误差,而 48 将产生负误差,因此 48 将向上舍入,48 将向下舍入。其余 4 种情况 ( 0.125, 0.375, 0.625, 0.875) 转换为binary64值不发生变化的格式,然后采用庄家舍入规则进行round四舍五入0.1250.625向下舍入和向上0.3750.875入。

符号。在这里和下面,我使用伪数学符号,而不是 Python 符号:^表示求幂而不是按位异或,/表示精确除法,而不是浮点除法。

例子

认为n = 11。然后我们考虑序列1/22, 3/22, ..., 21/22。以十进制表示的精确值有一个很好的简单的循环形式:

 1/22 = 0.04545454545454545...
 3/22 = 0.13636363636363636...
 5/22 = 0.22727272727272727...
 7/22 = 0.31818181818181818...
 9/22 = 0.40909090909090909...
11/22 = 0.50000000000000000...
13/22 = 0.59090909090909090...
15/22 = 0.68181818181818181...
17/22 = 0.77272727272727272...
19/22 = 0.86363636363636363...
21/22 = 0.95454545454545454...
Run Code Online (Sandbox Code Playgroud)

最接近的可精确表示的 IEEE 754 二进制 64 浮点值是:

 1/22 -> 0.04545454545454545580707161889222334139049053192138671875
 3/22 -> 0.13636363636363635354342704886221326887607574462890625
 5/22 -> 0.2272727272727272651575702866466599516570568084716796875
 7/22 -> 0.318181818181818176771713524431106634438037872314453125
 9/22 -> 0.409090909090909116141432377844466827809810638427734375
11/22 -> 0.5
13/22 -> 0.59090909090909093936971885341336019337177276611328125
15/22 -> 0.68181818181818176771713524431106634438037872314453125
17/22 -> 0.7727272727272727070868540977244265377521514892578125
19/22 -> 0.86363636363636364645657295113778673112392425537109375
21/22 -> 0.954545454545454585826291804551146924495697021484375
Run Code Online (Sandbox Code Playgroud)

通过直接检查我们发现,转换为float时,1/22、9/22、13/22、19/22和21/22向上取整,而3/22、5/22、7/22、15/22 17/22 向下舍入。(11/22 已经可以精确表示,因此不会发生舍入。)因此 11 个值中的 5 个被向上舍入,5 个被向下舍入。声称无论 的值如何,都会出现这种完美的平衡n

计算实验

对于那些可能更相信数值实验而不是形式证明的人,这里有一些代码(Python 语言)。

首先,让我们使用 Pythonfractions模块编写一个函数来创建我们感兴趣的序列:

from fractions import Fraction

def sequence(n):
    """ [1/2n, 3/2n, ..., (2n-1)/2n] """
    return [Fraction(2*i+1, 2*n) for i in range(n)]
Run Code Online (Sandbox Code Playgroud)

接下来,这是一个计算给定分数 的“舍入方向”的函数f,我们将其定义为1如果最接近的浮点数f大于f-1如果它较小,并且0如果它相等(即,如果f事实证明可以完全表示)采用 IEEE 754 二进制 64 格式)。Fraction请注意,在典型的使用 IEEE 754 的机器上,从到 的转换float是正确舍入的,并且 a和 aroundTiesToEven之间的顺序比较是使用所涉及数字的精确值来计算的。Fractionfloat

def rounding_direction(f):
    """ 1 if float(f) > f, -1 if float(f) < f, 0 otherwise """
    x = float(f)
    if x > f:
        return 1
    elif x < f:
        return -1
    else:
        return 0
Run Code Online (Sandbox Code Playgroud)

现在要计算给定序列的各种舍入方向,最简单的方法是使用collections.Counter

from collections import Counter

def round_direction_counts(n):
    """ Count of rounding directions for sequence(n). """
    return Counter(rounding_direction(value)
                   for value in sequence(n))
Run Code Online (Sandbox Code Playgroud)

现在我们可以输入任何我们想要观察的整数,以确保 的计数1始终与 的计数相匹配-1。这里有一些例子,从开始n = 100整个事情的例子开始:

>>> round_direction_counts(100)
Counter({1: 48, -1: 48, 0: 4})
>>> round_direction_counts(237)
Counter({-1: 118, 1: 118, 0: 1})
>>> round_direction_counts(24)
Counter({-1: 8, 0: 8, 1: 8})
>>> round_direction_counts(11523)
Counter({1: 5761, -1: 5761, 0: 1})
Run Code Online (Sandbox Code Playgroud)

上面的代码未经优化且相当慢,但我用它来运行测试n = 50000并检查每种情况下的计数是否平衡。

另外,这是一种可视化小舍入的简单方法n:它生成一个字符串,其中包含+向上舍入的情况、-向下舍入的情况以及.可精确表示的情况。+所以我们的定理说每个签名的字符数与字符数相同-

def signature(n):
    """ String visualising rounding directions for given n. """
    return "".join(".+-"[rounding_direction(value)]
                   for value in sequence(n))
Run Code Online (Sandbox Code Playgroud)

还有一些例子,证明没有立即明显的模式:

>>> signature(10)
'+-.-+++.--'
>>> signature(11)
'+---+.+--++'
>>> signature(23)
'---+++-+-+-.-++--++--++'
>>> signature(59)
'-+-+++--+--+-+++---++---+++--.-+-+--+-+--+-+-++-+-++-+-++-+'
>>> signature(50)
'+-++-++-++-+.+--+--+--+--+++---+++---.+++---+++---'
Run Code Online (Sandbox Code Playgroud)

声明的证明

我给出的原始证明过于复杂。根据蒂姆·彼得斯的建议,我意识到还有一个更简单的方法。如果您确实感兴趣,可以在编辑历史记录中找到旧的。

证明基于三个简单的观察。其中两个是浮点事实;第三是数论观察。

观察 1.对于任何(非小、非大)正分数xx以“相同方式”舍入2x

如果y是最接近 的二进制64 浮点x,则2y是最接近 的二进制64 浮点2x。因此,如果x向上舍入,则 也进行舍入2x,如果x向下舍入,则 也进行舍入2x。如果x是完全可表示的,则 也是2x

小字:“非微小、非巨大”应解释为我们避免 IEEE 754 二进制 64 指数范围的极端情况。严格来说,上述语句适用于x区间 中的所有内容[-2^1022, 2^1023)。在该范围的顶端有一个涉及无穷大的极端情况需要小心:如果x舍入到2^1023,则2x舍入到inf,因此该语句在该极端情况下仍然成立。

观察 1 意味着(再次假设避免了下溢和溢出),我们可以x按任意 2 的幂缩放任何分数,而不会影响转换为二进制 64 时的舍入方向。

观察 2.如果x是闭区间中的分数[1, 2],则以3 - x相反的方式舍入为x

这是因为 ify是最接近 的浮点数x(这意味着y也必须在区间 中[1.0, 2.0]),那么由于 内浮点数的均匀间距[1, 2]3 - y也可以精确表示,并且是最接近 的浮点数3 - x。这甚至适用于“最接近”的 roundTiesToEven 定义下的关系,因为 的最后一位y是偶数当且仅当 的最后一位3 - y是。

因此,如果x向上舍入(即y大于x),则3 - y小于3 - x,因此3 - x向下舍入。类似地,如果x可以精确表示,则 也可以3 - x

观察 3.分数序列1/2n, 3/2n, 5/2n, ..., (2n-1)/2n等于序列n/n, (n+1)/n, (n+2)/n, ..., (2n-1)/n,直至按 2 的幂缩放并重新排序。

这只是一个更简单的语句的缩放版本,即1, 3, 5, ..., 2n-1整数序列等于序列n, n+1, ..., 2n-1,最多可按 2 的幂缩放并重新排序。该语句也许最容易从相反的方向看到:从序列 开始n, n+1, n+2, ...,2n-1,然后将每个整数除以其最大的二次幂除数。在每种情况下,剩下的都必须是小于 的奇数2n,并且很容易看出这样的奇数不能出现两次,因此通过计数,我们必须按1, 3, 5, ..., 2n - 1某种顺序获得 中的每个奇数。

有了这三个观察结果,我们就可以完成证明了。结合观察1和观察3,我们得到 的累积舍入方向(即向上舍入、向下舍入、保持不变的总数)与 的1/2n, 3/2n, ..., (2n-1)/2n累积舍入方向完全匹配n/n, (n+1)/n, ..., (2n-1)/n

现在n/n正好是一,所以是完全可以表示的。在偶数的情况下n3/2也出现在这个序列中,并且是完全可以表示的。其余的值可以相互配对,总计为3(n+1)/n与 配对(2n-1)/n、与(n+2)/n配对(2n-2)/n等等。现在根据观察 2,在每一对中,要么一个值向上舍入,一个值向下舍入,要么两个值都可以精确表示。

因此,该序列n/n, (n+1)/2n, ..., (2n-1)/n具有与向上舍入情况完全相同数量的向下舍入情况,因此原始序列1/2n, 3/2n, ..., (2n-1)/2n具有与向上舍入情况完全相同的数量的向下舍入情况。这样就完成了证明。

注意:原始语句中对 的大小的限制n是为了确保我们的序列元素都不位于低于正常范围内,以便可以使用观察 1。最小的正二进制 64 正常值为2^-1022,因此我们的证明适用于所有n <= 2^1021