04 抹一次，和抹一千万次

▷ 猜一下差多少

一段最普通不过的循环：

let s = 0;
for (let i = 0; i < 10000000; i++) s += 0.1;   // 一千万次
console.log(s);                                 // 应该是 1000000

问：打出来的是什么？

A 正好 1000000。误差正负抵消，最后回到原位 B 1000000.0000000001 左右。差一两个 ulp C 999999.9998389754。差了大约万分之一 D 差好几个百分点

再猜第二个：如果把这一千万个数换成随机的正数（不再全是同一个 0.1），误差会变大还是变小？

先看结果

得到   = 999999.999838975374586880207061767578125    ← 这是精确值，不是显示截断
应为   = 1000000
绝对误差 = ⟨−0.0001610246254131198⟩
相对误差 = 1.61 × 10⁻¹⁰
差了     = 1 383 191 个 ulp        （ulp(1e6) = 1.1641532182693481e−10）

一百三十八万格。而每一次加法最多只抹掉半格。

为什么会这样？因为累加器一直在长大，而它长大的同时，脚下的刻度也在变粗。

循环刚开始时，s 只有几十，ulp 大约是 10⁻¹⁴，抹一次几乎不值一提。循环到一半，s 已经到了 50 万，ulp 变成 5.8×10⁻¹¹——抹一次的量比开头大了三个数量级。越加越粗，越粗越抹。

√n 还是 n：全看方向

现在回答第二个问题，也是这一章真正的判断。

每一次加法抹掉一点。这些「一点」加起来，长成什么样，取决于它们的方向是不是一致：

情形	抹零的方向	n 次之后的误差	例子
同一个数加很多遍	完全一致——每次都朝同一边偏	按 n 长	爱国者的时钟：每 tick 都少 9.54×10⁻⁸ 秒
各不相同的数	随机——一半向上一半向下	按 √n 长	加一百万个不同的测量值
最坏情况（教科书界）	假设全部同向	按 n 长	误差分析里的上界，实际很少达到

实测一下随机那一档。下面这张表的每一行，是把 n 个各不相同的数（exp(U(−1,1))，故意选了个二进制上不友好的分布）顺序累加，和精确和比对，跑 40 组取均方根：

n	实测误差（ulp 格数，均方根）	√n	最坏上界 n
100	2.2	10	100
1 000	4.7	32	1 000
10 000	17.3	100	10 000
100 000	59.9	316	100 000
1 000 000	147.3	1 000	1 000 000

n 涨一万倍（100 → 10⁶），误差只涨了 67 倍——比 √n 预测的 100 倍还略好一点，离最坏情况的 10 000 倍差了两个数量级。随机游走这件事，是浮点数最大的一份运气。

而 s += 0.1 那个循环拿不到这份运气：每一次加的都是同一个数，抹零的方向也就每一次都一样。爱国者导弹的三分之一秒，就是这么攒出来的。

◆ 主线

看到一个循环里的累加，先问两个问题：

累加器会长到多大？它长得越大，脚下的刻度越粗，每一步抹掉的越多。
加的东西是不是每次都一样？一样 ⇒ 误差按 n 长；各不相同 ⇒ 按 √n 长。

两个问题都答「糟」的循环，就是这本书里最经典的地雷：规则性 + 长时间运行。而这恰恰是嵌入式、控制、计时、动画、模拟里最常见的那种循环。

顺序也算数：一个 726 比 32 的例子

同一堆数，只是换个顺序加，误差可以差二十倍。看调和级数 1 + 1/2 + 1/3 + … + 1/10⁷：

顺序	结果	差几个 ulp
从大到小（1, 1/2, 1/3, …）	16.695311365857272	726
从小到大（…, 1/3, 1/2, 1）	16.695311365859965	32
精确（正确舍入）	16.695311365859851	0

道理和上面那张图是同一个：从小到大加，累加器在很长一段时间里都还很小，脚下的刻度还很细；从大到小加，累加器一上来就冲到最大，之后每一步都踩在最粗的格子上。

「先加小的」因此成了一条流传很广的建议。下一章之后你会发现它不太靠得住——第 8 章会给一份数据，在那上面「先排序再加」赢了，但赢的理由很可耻；而在另一些数据上它会输给一个既不排序也不补偿的办法。

一个反例：有时候一次都不抹

顺手记一条，因为它能帮你判断「什么时候可以完全放心」。

把一百万个「分母是 2 的幂」的数（比如 k / 2³¹ 这种）加起来，实测误差是 0 个 ulp——一次都没抹。

原因不神秘：这些数全都精确地落在格子上，它们的和也落在格子上（只要总和不超过 2⁵³ 个最小单位）。没有东西需要被抹。

这条反例的实用价值是：如果你的量本来就是「某个固定小单位的整数倍」——分、毫秒、像素、格子数——那么把它当整数处理，累加就是精确的。这是第 17 章「钱不能用 double」的理论依据。

⌨ 自己跑一遍

把「误差按 √n 长」这件事亲手量一遍。关键是要有一个精确的参照——Python 的 math.fsum 正好就是（它给出正确舍入的和，等价于用无限精度算完再存一次）：

import math, random, struct
def ulps(a, b):                      # 两个 double 之间隔了几格
    o = lambda x: (lambda u: u if u >= 0 else -0x8000000000000000 - u)(
        struct.unpack('<q', struct.pack('<d', x))[0])
    return abs(o(a) - o(b))

for n in (100, 1000, 10_000, 100_000, 1_000_000):
    random.seed(42)
    a = [math.exp(random.uniform(-1, 1)) for _ in range(n)]
    naive = 0.0
    for v in a: naive += v
    print(f"n={n:>9}  差 {ulps(naive, math.fsum(a)):>5} 格   √n = {n ** 0.5:.0f}")

s = 0.0
for _ in range(10_000_000): s += 0.1
print('累加 0.1 一千万次 =', repr(s), '  差', ulps(s, 1e6), '格')

最后一行会打出 999999.9998389754 差 1383191 格。前面那几行的具体数字会随种子变，但增长速度一定贴着 √n，不会贴着 n。

python3 sqrtn.py # 一千万次循环，纯 Python 大概要跑一两秒

JS 版本几乎一样，把 math.fsum 换成第 8 章会给的精确求和即可。要在线跑用 python.org/shell（一千万次可能会超时，把它改成一百万次，误差约 138 319 格）。

▸ 在现实里

温哥华证券交易所，1982–1983。这是累积抹零最有名的一桩公案。该所 1982 年 1 月推出一个新股指，起点定为 1000.000。指数每成交一笔就重算一次，而重算的结果被截断到小数点后三位，不是四舍五入——截断意味着方向永远向下，也就是上面那张表里最糟的一档。

每次截断平均少 0.0005 点，每天重算约 3000 次，一天就少 1.5 点；二十二个月大约 484 个交易日，累计约 700 点。1983 年 11 月，指数读数是 524.811，而按正确舍入重算的真值是 1098.892——差 574 点，和这个量级的估算对得上。修正当天，指数「一夜上涨」574 点。

动画与游戏循环。t += dt 是最常见的写法，也正是爱国者那个写法。跑几个小时之后，游戏内时间和真实时间会缓慢分家。稳妥的做法是 t = frameCount * dt，或者干脆用整数微秒计时。

训练里的滑动平均。动量、EMA、Adam 的一阶二阶矩，本质都是长时间累加。它们通常没事，因为 m = β·m + (1−β)·g 是指数遗忘的——老误差会被 β 一遍遍衰减掉，不会无限攒。这是一条很有用的通则：带遗忘的累加不积累误差，不带遗忘的才积累。

✗ 这个直觉是错的

「误差有正有负，加起来会互相抵消，所以长期跑不会有问题。」

「有正有负」这个前提要看情况。每次加同一个数时，抹零的方向是被这个数的二进制展开钉死的，完全不随机——爱国者每一次都向下，温哥华每一次都向下。

就算方向真的随机，也不是「抵消」，是随机游走：n 步之后离原点的期望距离是 √n 步，不是 0 步。一千万次之后，误差是单次误差的三千倍，不是零。

正确的问法：我这个循环里，误差是按 n 长、按 √n 长，还是根本不长？三个答案对应三种完全不同的应对——按 n 长要改写法（换成乘一次），按 √n 长要看 n 到底多大（一百万次通常无所谓），根本不长的就可以放着不管。

◇ 对账

正确答案是 C：999999.999838975374586880207061767578125，差 1 383 191 个 ulp，相对误差 1.61 × 10⁻¹⁰。换成随机的正数，误差会小得多（同样一千万次，只有几百个 ulp）。

A 「误差正负抵消」——这是最流行的那条错误直觉，上面那一栏专门在讲它。0.1 存进 double 之后是偏大的（大 5.55×10⁻¹⁸），可累加的结果却偏小：因为主导误差来自「累加器长大之后每一步的抹零」，不是来自 0.1 本身。这两件事的方向不一定一致，但都不是随机的。 B 「差一两个 ulp」——这是把一千万次运算当成了一次。单次抹零确实不超过半个 ulp，但这里有一千万次，而且它们同向。差了六个数量级。 D 「差好几个百分点」——过头了。相对误差是 1.6×10⁻¹⁰，也就是百亿分之一点六。这个数在大多数场合完全无害。这一章要讲的不是「double 不能用」，而是它在什么条件下会从 10⁻¹⁶ 涨到 10⁻¹⁰，以及那个条件长什么样。同样的循环再跑一万倍长，就该报警了。

◈ 换一种写法

累积误差有三种改法，从便宜到贵：

for (…) s += x; —— 一路累加。误差按 √n 或 n 长。 ① 能不累加就不累加。s = n * x（每次加同一个数时）、t = frame * dt。误差降到一次。零成本。 ② 换成分治。把数组对半分，各自求和再相加。误差从 O(√n) 降到 O(√log n)，代码五行，速度不变——这就是 numpy 默认的做法。第 8 章展开。 ③ 补偿求和。用一个额外变量记住每一步被抹掉的量，下一步补回去。误差降到几乎为零，代价是慢两三倍。第 8 章会看到它也有失手的时候。

顺序也是免费的杠杆：从小到大加，上面那个调和级数从 726 格降到 32 格，一行 sort 的事。但这条建议有前提，第 8 章会给出它失效的场合。

这一章的一句话

累加的误差不看数据量，看两件事——累加器会长到多大（决定每一步抹多少），以及抹零的方向是否一致（决定它按 √n 长还是按 n 长）。

卷 I 到这里就结束了。你已经有了这本书的全部零件：一排不均匀的刻度、一次抹零、两种放大方式（相消和累积）。

卷 II 要做的是把「放大」这件事拆成两半——而这一拆，是整本书最值钱的一个动作。下一章先看第一半：有些问题天生就对输入敏感，跟你写什么代码毫无关系。有一个多项式，把它二十个系数里的一个改动 0.00000012，二十个根里会有十个当场离开实数轴，跑到复平面上去。