公式在IEEE 754中失败的可能性

Question

在我的电脑上，我可以检查一下

(0.1 + 0.2) + 0.3 == 0.1 + (0.2 + 0.3)

评估为False。

更一般而言，我可以使用以下模拟估算公式在统一且独立地选择时(a + b) + c == a + (b + c)大致失败17%的时间：a,b,c[0,1]

import numpy as np
import numexpr

np.random.seed(0)
formula = '(a + b) + c == a + (b + c)'


def failure_expectation(formula=formula, N=10**6):
    a, b, c = np.random.rand(3, N)
    return 1.0 - numexpr.evaluate(formula).mean()
# e.g. 0.171744

我想知道是否有可能手动得出这种概率，例如使用浮点标准中的定义以及对均匀分布的一些假设。

---

给定以下答案，我认为至少到现在为止，原始问题的以下部分是遥不可及的。

有没有一种工具可以在不运行模拟的情况下计算给定公式的失效概率。

可以假定公式很简单，例如，使用括号，加法，减法以及可能的乘除运算。

---

（以下内容可能是numpy随机数生成的产物，但似乎仍然很有趣。）

奖励问题基于NPE的观察。我们可以使用以下代码为一系列范围内的均匀分布生成故障概率[[-n,n] for n in range(100)]：

import pandas as pd

def failures_in_symmetric_interval(n):
    a, b, c = (np.random.rand(3, 10**4) - 0.5) * n
    return 1.0 - numexpr.evaluate(formula).mean()

s = pd.Series({
    n: failures_in_symmetric_interval(n)
    for n in range(100)
})

情节看起来像这样：

特别是，故障概率下降到0when n是的幂，2并且似乎具有分形模式。看起来每个“下降”的失败概率都等于先前的某个“峰值”。任何说明为何会发生这种情况的方法都很好！

Answer 1

手工评估这些东西绝对是可能的，但我所知道的唯一方法很乏味，并且涉及大量的具体情况枚举。

例如，对于确定 概率的具体示例(a + b) + c == a + (b + c)，该概率为 53/64，在机器 epsilon 的几倍之内。因此，不匹配的概率为 11/64，即 17.19% 左右，这与您从模拟中观察到的结果一致。

首先，请注意，在这种特殊情况下有一个主要的简化因素，那就是 Python 和 NumPy 的“uniform-on-[0, 1]”随机数始终采用 中n/2**53某个整数 n 的形式range(2**53)，并且在约束范围内对于底层 Mersenne Twister PRNG，每个这样的数字出现的可能性是相等的。2**62由于范围内大约有IEEE 754 二进制 64 位可表示值[0.0, 1.0]，这意味着这些 IEEE 754 值中的绝大多数不是由random.random()( 或np.random.rand()) 生成的。这一事实极大地简化了分析，但也意味着它有点作弊。

这是一个不完整的草图，只是为了让您了解所涉及的内容。为了计算 53/64 的值，我必须分为五个不同的情况：

1. a + b < 1 且 b + c < 1 的情况。在这种情况下，a + b 和 b + c 的计算均无错误，因此 (a + b) + c 和 a + (b + c)两者都给出了最接近准确结果的浮动，像往常一样将平局舍入为偶数。所以在这种情况下，达成一致的概率为 1。

2. a + b < 1 且 b + c >= 1 的情况。此处 (a + b) + c 将是真实总和的正确舍入值，但 a + (b + c) 可能不是。我们可以根据 a、b 和 c 的最低有效位的奇偶校验进一步划分子情况。让我们滥用术语，如果 n 为奇数形式，则称其为“奇数” n/2**53；如果 n 为偶数形式，则称其为“n/2**53偶数”，b 和 c 也类似。如果 b 和 c 具有相同的奇偶性（一半的情况下会发生），则 (b + c) 会被精确计算，并且 a + (b + c) 必须匹配 (a + b) + c。对于其他情况，每种情况达成一致的概率为 1/2；细节都非常相似，但例如，在 a 为奇数、b 为奇数、c 为偶数的情况下，(a + b) + c 被精确计算，而在计算 a + (b + c) 时，我们会产生两个舍入误差，每个量值均准确2**-53。如果这两个错误方向相反，它们就会抵消，我们就达成一致。如果没有，我们就不会。总体而言，在这种情况下达成一致的概率为 3/4。

3. a + b >= 1 且 b + c < 1 的情况。交换 a 和 c 的角色后，这与之前的情况相同；达成一致的概率又是 3/4。

4. a + b >= 1 且 b + c >= 1，但 a + b + c < 2。同样，我们可以拆分 a、b 和 c 的奇偶性，并依次查看所得的 8 种情况。对于偶-偶-偶和奇-奇-奇的情况，我们总是能达成一致。对于奇-偶-奇的情况，一致的概率是 3/4（通过进一步的子分析）。对于所有其他情况，该值为 1/2。对于这种情况，将这些放在一起得到的总概率为 21/32。

5. 情况 a + b + c >= 2。在这种情况下，由于我们将最终结果舍入为四倍的倍数2**-53，因此不仅需要查看 a、b 和 c 的奇偶校验，还需要查看最后两个有效位。我就不跟你讲那些血淋淋的细节了，但达成一致的概率是 13/16。

最后，我们可以将所有这些案例放在一起。为此，我们还需要知道三元组 (a, b, c) 在每种情况下落地的概率。a + b < 1 和 b + c < 1 的概率是由 0 <= a, b, c <= 1, a + b < 1, b + c < 1 描述的方金字塔的体积，其中是 1/3。其他四种情况的概率可以看出（通过一点立体几何，或者通过设置合适的积分）各为 1/6。

所以我们的总计是1/3 * 1 + 1/6 * 3/4 + 1/6 * 3/4 + 1/6 * 21/32 + 1/6 * 13/16，结果是53/64如所声称的那样。

最后一点：53/64 几乎肯定不是正确的答案 - 为了获得完全准确的答案，我们需要小心 a + b、b + c 或 a + b + c 的所有极端情况达到二进制边界（1.0 或 2.0）。当然可以改进上述方法来准确计算有多少2**109可能的三元组 (a, b, c) 满足 (a + b) + c == a + (b + c)，但不是在时机成熟之前让我去睡觉。但极端情况应占案例总数的 数量级1/2**53，因此我们对 53/64 的估计应至少精确到小数点后 15 位。

当然，上面遗漏了很多细节，但我希望它能提供一些关于如何做到这一点的想法。