解释区间 [0, 1] 中表观关系在舍入方向上的惊人奇偶性

Explain a surprising parity in the rounding direction of apparent ties in the interval [0, 1]

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

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

>>> values = [n/1000 for n in range(5, 1000, 10)]

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

>>> 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]

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

>>> sum(round(value, 2) > value for value in values)
50

原来这 100 个数字正好有一半被四舍五入了。

如果您不知道 Python 在幕后使用二进制浮点数,这个结果就不足为奇了。毕竟,Python 的 documentation 清楚地表明 round 函数使用 round-ties-to-even (a.k.a. Banker's rounding) 作为其舍入模式,因此您希望这些值交替向上和向下舍入。

但是 Python 在后台使用二进制浮点数,这意味着除了少数例外(即 0.1250.3750.6250.875),这些值 不是 精确关系,而只是对这些关系的非常好的二元近似。毫不奇怪,仔细检查舍入结果会发现这些值 而不是 交替上下舍入。相反,每个值向上或向下取整取决于二进制近似值恰好落在十进制值的哪一侧。所以没有 a priori 理由期望恰好一半的值向上舍入,恰好一半向下舍入。这让我们得到的结果正好是 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

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

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

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

备注

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

不是答案,只是想充实一下它的疑惑所在。它当然不是“随机的”,但要注意这还不够 ;-) 具体请看 2 位数的情况:

>>> from decimal import Decimal as D
>>> for i in range(5, 100, 10):
...     print('%2d' % i, D(i / 100))
    
 5 0.05000000000000000277555756156289135105907917022705078125
15 0.1499999999999999944488848768742172978818416595458984375
25 0.25
35 0.34999999999999997779553950749686919152736663818359375
45 0.450000000000000011102230246251565404236316680908203125
55 0.5500000000000000444089209850062616169452667236328125
65 0.65000000000000002220446049250313080847263336181640625
75 0.75
85 0.84999999999999997779553950749686919152736663818359375
95 0.9499999999999999555910790149937383830547332763671875

现在您可以将 i/100(100-i)/100 配对,它们的数学总和恰好为 1。所以这对,在上面,5 与 95,15 与 85,依此类推。 5 舍入的精确机器值,而 95 舍入的精确机器值,这是“预期的”:如果真实总和为 1,并且一个加数“向上舍入”,那么另一个肯定“向下舍入”。

但情况并非总是如此。 15 和 85 都向下舍入,25 和 75 是混合,35 和 65 是混合,但是 45 和 55 都向上舍入。

是什么在起作用使总的“上升”和“下降”情况完全平衡? Mark 表明它们对 10**310**710**9 有效,我也验证了指数 2、4、5、6、8、10 和 11 的精确平衡。

一个令人费解的线索

非常精致。如果我们不是除以 10**n,而是乘以它的倒数会怎么样。对比一下上面的:

>>> for i in range(5, 100, 10):
...     print('%2d' % i, D(i * (1 / 100)))

 5 0.05000000000000000277555756156289135105907917022705078125
15 0.1499999999999999944488848768742172978818416595458984375
25 0.25
35 0.350000000000000033306690738754696212708950042724609375
45 0.450000000000000011102230246251565404236316680908203125
55 0.5500000000000000444089209850062616169452667236328125
65 0.65000000000000002220446049250313080847263336181640625
75 0.75
85 0.84999999999999997779553950749686919152736663818359375
95 0.95000000000000006661338147750939242541790008544921875

现在有 7 个(而不是 5 个)案例汇总。

对于10**3,64(而不是50)四舍五入;对于 10**4,828(而不是 500),对于 10**5,9763(而不是 5000);等等。因此,在计算 i/10**n.

时,不超过一个舍入误差是 重要的事情

事实证明,可以证明更强大的东西,这与十进制表示或十进制舍入无关。这是更有力的说法:

Theorem. Choose a positive integer n <= 2^1021, and consider the sequence of length n consisting of the fractions 1/2n, 3/2n, 5/2n, ..., (2n-1)/2n. Convert each fraction to the nearest IEEE 754 binary64 floating-point value, using the IEEE 754 roundTiesToEven rounding direction. Then the number of fractions for which the converted value is larger than the original fraction will exactly equal the number of fractions for which the converted value is smaller than the original fraction.

涉及浮点序列 [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/223/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...

最接近的可精确表示的 IEEE 754 binary64 浮点值是:

 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

而且我们通过直接检查发现,当转换为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 中)。

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

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)]

接下来,这里有一个函数来计算给定分数 f 的“舍入方向”,如果最接近 f 的浮点数更大,我们将其定义为 1f-1(如果它更小)和 0(如果它相等)(即,如果 f 结果完全可以用 IEEE 754 binary64 格式表示)。请注意,在使用典型 IEEE 754 的机器上,从 Fractionfloat 的转换在 roundTiesToEven 下正确舍入,并且 Fraction 和 [= 之间的顺序比较51=] 是使用所涉及数字的精确值计算的。

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

现在要计算给定序列的各种舍入方向,最简单的方法是使用 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))

现在我们可以输入任何我们喜欢的整数来观察 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})

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

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

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

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

>>> signature(10)
'+-.-+++.--'
>>> signature(11)
'+---+.+--++'
>>> signature(23)
'---+++-+-+-.-++--++--++'
>>> signature(59)
'-+-+++--+--+-+++---++---+++--.-+-+--+-+--+-+-++-+-++-+-++-+'
>>> signature(50)
'+-++-++-++-+.+--+--+--+--+++---+++---.+++---+++---'

陈述证明

我给出的原始证明过于复杂。根据 Tim Peters 的建议,我意识到有一个更简单的方法。如果您真的感兴趣,您可以在编辑历史中找到旧的。

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

Observation 1. For any (non-tiny, non-huge) positive fraction x, x rounds "the same way" as 2x.

如果 y 是最接近 x 的 binary64 浮点数,那么 2y 是最接近 2x 的 binary64 浮点数。因此,如果 x 向上舍入,2x 也会向上舍入,如果 x 向下舍入,2x 也会向下舍入。如果 x 是可精确表示的,那么 2x.

也是

小字:“非微小,非巨大”应解释为意味着我们避免了 IEEE 754 binary64 指数范围的极端情况。严格来说,上面的说法适用于区间[-2^1022, 2^1023)内的所有x。在该范围的顶端有一个涉及无穷大的极端情况需要注意:如果 x 舍入为 2^1023,则 2x 舍入为 inf,因此声明仍然适用于那个极端情况。

观察结果 1 意味着(再次假设避免了下溢和上溢),我们可以将任何分数 x 缩放为 2 的任意幂,而不会影响它在转换为 binary64 时舍入的方向。

Observation 2. If x is a fraction in the closed interval [1, 2], then 3 - x rounds the opposite way to x.

这是因为如果 y 是最接近 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.

也是可表示的

Observation 3. The sequence 1/2n, 3/2n, 5/2n, ..., (2n-1)/2n of fractions is equal to the sequence n/n, (n+1)/n, (n+2)/n, ..., (2n-1)/n, up to scaling by powers of two and reordering.

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

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

现在 n/n 恰好是一,因此可以精确表示。在 n 是偶数的情况下, 3/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的大小限制是为了保证我们的序列元素none在次正规范围内,这样Observation 1就可以使用了。最小的正 binary64 正常值是 2^-1022,所以我们的证明适用于所有 n <= 2^1021.

为了具体起见,我将通过 Mark 的解释(正如我 mod 在评论中证实的那样)来解释在我发布详尽结果的 2 位数案例中看到的所有内容。

那里我们正在寻找 i / 100 以获得 i in range(5, 100, 10),它正在寻找 (10*i + 5) / 100 以获得 i in range(10),这是相同的(将分子和分母除以 5 ) 就像在 (2*i + 1) / 20 中查看 i in range(10).

“重新缩放技巧”包括将每个分子向左移动,直到它成为 >= 10。转换为二进制浮点数时,这与四舍五入无关紧要! 2 的幂的因数只影响指数,而不影响有效位(假设我们保持在正常范围内)。通过平移,我们将所有分子调整到 range(10, 20),因此当除以 20 时,我们得到半开范围 [0.5, 1.0) 中的有效分数,它们都具有相同的 2 次幂指数.

唯一的 k 使得 2**52 <= 10/20 * 2**k = 1/2 * 2**k < 2**53k=53(因此商的整数部分具有 53 位精度的 IEEE-754 双精度),所以我们'重新查看 i * 2**53 / 20 形式的转换比率 i in range(10, 20).

现在对于任何 n,将 n 表示为 2**t * o 其中 o 是奇数:

i * 2**k = j * 2**k (mod 2*n)当且仅当

i * 2**k = j * 2**k (mod 2**(t+1) * o)当且仅当(假设k >= t+1

i * 2**(k-t-1) = j * 2**(k-t-1) (mod o) iff(o 是奇数,所以与 2**(k-t-1) 互质)

i = j (mod o)

range(n, 2*n)n 个连续整数,因此 o 个元素的每个子切片 mod o 包含每个残基 class mod o 恰好一次,每个残基 class modulo orange(n, 2*n) 中恰好出现 2**t 次。最后一点在这里最重要,因为重新缩放技巧给我们留下了 range(n, 2*n).

的排列

我们正在使用 n = 10 = 2**1 * 5i * 2**53 / 20 = i * 2**51 / 5。在

q, r = divmod(i * 2**51, 5)

q为53位有效数,r为余数。如果余数为0,则q是精确的;如果余数为 1 或 2,则 q 稍微太小(“向下舍入”),如果余数为 3 或 4,则硬件将通过向 q 加 1 来“向上舍入”。但是这里我们不关心q,我们只想知道会发生哪个舍入动作,所以r才是我们关心的。

现在pow(2, 51, 5) = 3,所以,modulo 5,乘以2**51等于乘以3。取range(1, 20, 2)中的奇数并做重新缩放技巧,将所有内容压缩为 range(10, 20),然后乘以 2**51(与 3 相同),并找到余数 mod 5:

1  -> 16, * 3 % 5 = 3 up
3  -> 12, * 3 % 5 = 1 down
5  -> 10, * 3 % 5 = 0 exact
7  -> 14, * 3 % 5 = 2 down
9  -> 18, * 3 % 5 = 4 up
11 -> 11, * 3 % 5 = 3 up
13 -> 13, * 3 % 5 = 4 up
15 -> 15, * 3 % 5 = 0 exact
17 -> 17, * 3 % 5 = 1 down
19 -> 19, * 3 % 5 = 2 down

与之前发布的详尽结果相符。

不是答案,而是进一步评论。

我的假设是:

  • 原始 n/1000 的结果将四舍五入为小于或大于精确的小数值,方法是计算一个额外的精度位,然后使用 0 或 1在该额外位中确定是向上舍入还是向下舍入(相当于银行家舍入的二进制)

  • round 以某种方式将该值与精确的分数值进行比较,或者至少表现得好像它正在这样做(例如,在使用 more内部精度位,至少对于乘法)

  • 根据 精确 分数的一半可以向上取整而另一半向下取整的问题来信任它

如果是这样的话,那么这个问题就相当于说:

  • 如果你把分数写成二进制数,有多少分数在第 i 位(其中 i第 'th 位对应于 存储的最后一位之后的位置,根据我的假设,它将用于决定以哪种方式舍入数字)

考虑到这一点,这里有一些代码可以计算任意精度二进制数,然后对这些二进制数的第 i 位求和(对于非精确情况)和加上一半的非精确案例。

def get_binimal(x, y, places=100,
                normalise=True):
    """
    returns a 2-tuple containing: 
        - x/y as a binimal, e.g. for 
            x=3, y=4 it would be 110000000...
        - whether it is an exact fraction (in that example, True)

    if normalise=True then give fractional part of binimal that starts
    with 1. (i.e. IEEE mantissa)
    """
    if x > y:
        raise ValueError("x > y not supported")
    frac = ""
    val = x
    exact = False
    seen_one = False
    if normalise:
        places += 1  # allow for value which is always 1 (remove later)
    while len(frac) < places:
        val *= 2
        if val >= y:
            frac += "1"
            val -= y
            seen_one = True
            if val == 0:
                exact = True
        else:
            if seen_one or not normalise:
                frac += "0"
    if normalise:
        frac = frac[1:]  # discard the initial 1
    return (frac, exact)


places = 100

n_exact = 0
n = 100
divisor = n * 10
binimals = []
for x in range(5, divisor, 10):
    binimal, exact = get_binimal(x, divisor, places, True)
    print(binimal, exact, x, n)
    if exact:
        n_exact += 1
    else:
        binimals.append(binimal)
        
for i in range(places):
    print(i, n_exact // 2 + sum((b[i] == "1") for b in binimals))

运行 这个程序给出了例子:

0 50
1 50
2 50
3 50
4 50
5 50
6 50
7 50
8 50
... etc ...

一些观察结果,即:

  • 已确认(根据显示的结果加上对 n 的其他值的实验)这给出了与问题中观察到的相同的计数(即 n/2),所以上述假设似乎有效。

  • i 的值无关紧要,即 IEEE 64 位浮点数中的 53 个尾数位没有什么特别之处——任何其他长度都会给出相同的值。

  • 数字是否归一化无关紧要。请参阅我的 get_binimal 函数的 normalise 参数);如果将其设置为 True,则返回值类似于规范化的 IEEE 尾数,但计数不受影响。

显然二项展开将由重复序列组成,而 i 无关紧要的事实表明序列必须以这样的方式对齐,即 i =63=]i的数字总是相同的,因为重复序列的每个对齐都有相同的数字。

以 n=100 为例,并显示每个扩展的最后 20 位的计数(即位 80-99,因为我们要求 100 个位置)使用:

counts = collections.Counter([b[-20:] for b in binimals])
pprint.pprint(counts.items())

给出类似下面的内容,尽管在这里我手动编辑了顺序以便更清楚地显示重复序列:

[('00001010001111010111', 4),
 ('00010100011110101110', 4),
 ('00101000111101011100', 4),
 ('01010001111010111000', 4),
 ('10100011110101110000', 4),
 ('01000111101011100001', 4),
 ('10001111010111000010', 4),
 ('00011110101110000101', 4),
 ('00111101011100001010', 4),
 ('01111010111000010100', 4),
 ('11110101110000101000', 4),
 ('11101011100001010001', 4),
 ('11010111000010100011', 4),
 ('10101110000101000111', 4),
 ('01011100001010001111', 4),
 ('10111000010100011110', 4),
 ('01110000101000111101', 4),
 ('11100001010001111010', 4),
 ('11000010100011110101', 4),
 ('10000101000111101011', 4),

 ('00110011001100110011', 4),
 ('01100110011001100110', 4),
 ('11001100110011001100', 4),
 ('10011001100110011001', 4)]

有:

  • 20 位重复序列的 80 (=4 * 20) 个视图
  • 16 (=4 * 4) 个视图的 4 位重复序列对应于除以 5(例如 0.025 十进制 = (1/5) * 2^-3)
  • 4 个精确分数(未显示),例如小数 0.375 (= 3 * 2^-3)

正如我所说,这并未声称是完整答案

真正有趣的事情是这个结果似乎并没有被标准化数字所破坏。丢弃前导零肯定会改变各个分数的重复序列的对齐方式(通过改变位数来移动序列,具体取决于忽略了多少前导零),但这样做的方式是每次对齐的总计数被保存下来。我发现这可能是结果中最令人好奇的部分。

还有一件奇怪的事情 - 20 位重复序列由一个 10 位序列和其补码组成,例如以下两个相等数量的对齐将在每个位位置给出相同的总数:

10111000010100011110
01000111101011100001

4 位重复序列也类似。但结果似乎并不取决于此 - 相反,所有 20 个(和所有 4 个)比对的数量相等。