公式在 IEEE 754 中失败的概率

Question

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

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

计算为 False。

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

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

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

鉴于下面的答案，我认为原始问题的以下部分是遥不可及的，至少现在是这样。

Is there is a tool that computes the failure probability for a given formula without running a simulation.

Formulas can be assumed to be simple, e.g. involving the use of parentheses, addition, subtraction, and possibly multiplication and division.

（以下可能是 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)
})

情节看起来像这样：

特别是，当 n 是 2 的幂并且似乎具有分形图案时，失败概率下降到 0。看起来每个 "dip" 的失败概率都与之前的 "peak" 相同。任何关于为什么会发生这种情况的解释都会很棒！

Answer 1

当然可以手动评估这些东西，但我所知道的唯一方法是乏味的，并且涉及大量的个案枚举。

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

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

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

a + b < 1 和b + c < 1 的情况。在这种情况下，a + b 和b + c 都计算没有错误，并且(a + b) + c 和 a + (b + c) 因此都给出最接近精确结果的浮点数，像往常一样四舍五入到 even 。所以在这种情况下，一致的概率是1.
a + b < 1 和 b + c >= 1 的情况。这里 (a + b) + c 将是真实总和的正确舍入值，但 a + (b + c) 可能不是。我们可以根据 a、b 和 c 的最低有效位的奇偶性进一步划分为子情况。让我们滥用术语并称 a "odd" 如果它的形式为 n/2**53 且 n 为奇数，如果它的形式为 n/2**53 且 n 为偶数则称其为 "even"，对于 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 的概率达成一致。
a + b >= 1 且b + c < 1 的情况。这与之前交换a 和c 的角色后的情况相同；一致的概率又是 3/4.
a + b >= 1 且 b + c >= 1，但 a + b + c < 2。同样，可以拆分 a、b 和 c 的奇偶校验并查看依次产生8个案例。对于 even-even-even 和 odd-odd-odd 的情况，我们总是能达成一致。对于奇-偶-奇的情况，一致的概率结果为 3/4（通过进一步的子分析）。对于所有其他情况，它是 1/2。对于这种情况，将它们放在一起得到 21/32 的总概率。
情况 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 位。

当然，上面遗漏了很多细节，但我希望它能让您了解如何做到这一点。

公式在 IEEE 754 中失败的概率

Probability that a formula fails in IEEE 754

python

numpy

floating-point

probability

ieee-754