众所周知，浮点数是有一定误差的。经过多次计算，误差会累积。那么我们能不能从理论上分析误差累积的情况呢？

当然大多数情况误差都是可控的，只有比较极端的情况才会注意到误差的存在。

1. 浮点数

浮点数不介绍了。为了方便，浮点数的信息如下表所示：

类型	符号位	指数宽度	精度
float	1 bit	8 bit	23 bit
double	1 bit	11 bit	52 bit

2. 加法

假设有两个浮点数 A, B，且都是正数。

• A 累积误差是 $e_A$（A 和准确的 A 距离最大值，也可以叫做绝对误差），A 舍入误差是 $p_A$（A 可以表示 $A \pm p_A$ 范围内的数）。
• B 累积误差是 $e_B$，B 舍入误差是 $p_B$。

如果这两个数做加法，那么首先 $e_A, e_B$ 可以直接相加。

如果 $A > B$ 且 B 的一些有效位可能已经超出 A 的精度了，这就会有舍入误差 $p_A$。

如果 A, B 有效位互相不超出对方的精度，那么 A + B 一定会进位，所以还是会有舍入误差 $p_A$。

所以理论上累积误差 $e_{A+B}$ 应该是 $e_A+e_B+\max(p_A, p_B)$。

---

因为舍入误差取决于这个浮点数的大小，并不是定值，不好计算。我们用一个常数 u 来表示舍入误差和浮点数的比值。（这么做会稍微不精确一些）

$e_{A+B}=e_A+e_B+\max(A, B)\cdot u$

然后方便计算，把 max 改成加法。（这么做又会不精确一些）

$e_{A+B}=e_A+e_B+(A+B)\cdot u$

这个公式就和论文 [1] 基本一致了。

其他运算基本和加法一样的分析思路。

3. 例：求和

一个最简单的例子就是对 n 个数求和。（其实论文 [1] 提到 5 个方法，但是我没时间完整地读论文）

3.1. 朴素做法

最朴素的做法就是顺序地加起来，代码如下：

float sum(const std::vector<float> &a) {
    float res = 0;
    for (float i : a) {
        res += i;
    }
    return res;
}

这个做法会产生很大的误差。

前 2 个数相加的舍入误差是 $(a_0+a_1)\cdot u$，然后用 $a_0+a_1$ 和第三个数相加，舍入误差是 $(a_0+a_1+a_2)\cdot u$，以此类推。

（注意：用 $a_0+a_1$ 和第三个数相加，更准确的应该是用 $a_0+a_1+(a_0+a_1)\cdot u$ 和第三个数相加，舍入误差是 $(a_0+a_1+a_2)\cdot u + (a_0+a_1)\cdot u^2$。因为 $u^2$ 有点小，就忽略了）

最后把所有舍入误差加起来得到：

$\left[(n-1)a_0+(n-1)a_1+(n-2)a_2+(n-3)a_3+\ldots+a_{n-1}\right]\cdot u$

如果 $a_i$ 大小都约等于 1，那么误差约等于：（因为这是个特殊情况，具体情况还得具体分析）

$\left[\dfrac{n\cdot(1+n)}{2}-1\right]\cdot u$

可以看到误差增长是 $O(n^2)$ 的。

3.2. Pairwise 求和

一个较优的方案是 Pairwise 求和 (opens new window)，它通过分治的方法，计算左半部分的和，以及右半部分的和，将这两个值加起来作为答案。代码如下：

float sum(float *l, float *r) {
    if (l == r) { return 0; }
    if (l + 1 == r) { return *l; }
    float *m = l + (r - l) / 2;
    return sum(l, m) + sum(m, r);
}

分析一下这个算法，可以发现其实同一递归深度的舍入误差加起来，都是 $(a_0+a_1+...+a_{n-1})\cdot u$，而递归深度是 $\log n$，所以总的误差是：

$(a_0+a_1+...+a_{n-1})\cdot u\log n$

如果 $a_i$ 大小都约等于 1，那么误差约等于：

$un\log n$

误差增长是 $O(n\log n)$，是有一些改进的。

3.3. Kahan 求和

Kahan 求和 (opens new window) 是一种补偿的求和算法。代码如下：

float sum(const std::vector<float> &a) {
    float res = 0;
    float c = 0;
    for (auto y : a) {
        y -= c;
        float t = res + y;
        c = (t - res) - y;
        res = t;
    }
    return res;
}

个人理解相当于有了两个浮点数的精度，所以累积误差应该是朴素做法里的 u 换成 $u^2$。所以这个算法只能扩展 double 求和的精度，但是 double 的精度太高了，想不到合适的使用场景（ACM 卡精度题除外）。

4. 参考

[1] Higham N J. The accuracy of floating point summation[J]. SIAM Journal on Scientific Computing, 1993, 14(4): 783-799.