减法不产生误差
数值计算里最有名的坑叫「灾难性相消」，它的名字造成了一个持久的误解：好像减法是一种危险操作，会把误差搞出来。事实正好相反——两个相近的 double 相减，结果精确到一个 ulp 都不差，这是可以证明的定理。可就是这个最诚实的运算，能在一步之内让八位有效数字消失。
要算 √(x+1) − √x，在 x = 100 000 000（一亿）处。这是个再普通不过的表达式，中学生都会。照着写：
Math.sqrt(x + 1) - Math.sqrt(x)
问：这样算出来的答案，有几位有效数字是对的？（double 一共能给你 15 到 16 位）
顺便想一下：如果把 x 换成一万亿（10¹²），情况会变好还是变坏？
先证明减法是无辜的
有一条很短的定理，叫 Sterbenz 引理：
如果两个浮点数 a 和 b 满足 b/2 ≤ a ≤ 2b（也就是它们相差不到一倍），那么 a − b 可以被精确表示，减法的结果一个 ulp 都不会被抹。
为什么？想想上一章那把尺子。a 和 b 相差不到一倍，说明它们最多跨了一个段，两者所在位置的刻度间距最多差一倍。它们的差是刻度间距的整数倍，而这个差比 a 本身小，落在更靠近原点、刻度更密的地方。既然差值是粗刻度的整数倍，它当然也落在细刻度上——结果一定在格子上，不需要抹。
所以：两个相近的数相减，是浮点运算里最精确的一种操作。这不是修辞，是定理。
那八位有效数字是怎么没的？
误差不是被创造的，是被提拔的
把 x = 1e8 那一步拆开看。√(x+1) 和 √x 各自都要先存成 double，各自都被抹过一次：
√(100000001) = 10000.0000500000⟨00555883…⟩ ← 存进 double，误差约 1e−12
√(100000000) = 10000 ← 正好，精确
相减（Sterbenz：这一步精确，一格都不抹）：
= 0.0000500000⟨00555883…⟩
正确答案是 = 0.00004999999987⟨5000…⟩
↑ 从这一位开始就不对了
看清楚发生了什么。√(x+1) 那个数的绝对误差大约是 10⁻¹²——相对它自己（一万）来说是 10⁻¹⁶，非常小，完全合格。
然后减法把前面那四位「10000.」全部消掉了。剩下的结果只有 5×10⁻⁵ 那么大，可是那个 10⁻¹² 的绝对误差原封不动地留了下来。它相对于新结果的比例，一下子从 10⁻¹⁶ 变成了 2×10⁻⁸。
相减前：绝对误差 1e−12，相对误差 1e−16 ← 合格 相减后：绝对误差 1e−12（一点没变），相对误差 2e−8 ← 放大了一亿倍 # 减法一个新误差都没造。它只是把结果缩小了一亿倍， # 于是原本藏在小数点后面很远的那个误差，被顶到了前排。
这就是全书第一个「放大」的例子，也是整个卷 II 的引子。
丢几位可以算出来
这件事有一个很好记的定量版本：
计算 a − b 时，如果 b/a 接近 1，你丢掉的有效位数约等于 −log₂|1 − b/a|。
换成十进制的说法：a 和 b 前面有几位相同，就丢几位。相同得越多，丢得越多。
a = 1.0000000012345678
b = 1.0000000000000000
└────┬───┘
前 10 位一模一样
a − b = 1.2345677813385691e−9
└───┬───┘
只剩 8 位可信，后面全是 a 存进来时被抹掉的那点东西
回到开头的题：√(1e8+1) 和 √(1e8) 前面有 8 位相同（10000.00），所以丢 8 位。答案是 C。
把 x 换成 10¹²？那两个开方值前面相同的位数更多——实测只剩 5.1 位可信。x 越大越糟，这一点也很反直觉：大多数人会觉得数越大越「稳当」。
另一个经典：(1 − cos x) / x²
这个式子在物理和图形里到处都是（它的极限是 1/2）。看它在 x 变小时的表现：
| x | 直接照写 (1-Math.cos(x))/(x*x) | 改写后 |
|---|---|---|
| 10⁻³ | 0.49999995832550326 | 0.4999999583333347 |
| 10⁻⁵ | 0.5000000413701854 | 0.4999999999958333 |
| 10⁻⁷ | 0.4996003610813205 | 0.4999999999999996 |
| 10⁻⁸ | 0 | 0.5 |
| 10⁻⁹ | 0 | 0.5 |
x 小到 10⁻⁸ 时，cos x 已经是 0.99999999999999999…，它存进 double 之后就是 1。1 − 1 = 0。整个式子返回 0，而正确答案是 0.5。
注意这里的可怕之处：没有报错，没有 NaN，没有 Infinity。它返回了一个格式完全正确、类型完全正确、看起来无比正常的 0。这类失败是这本书里最值得警惕的一种，它不会在测试里跳出来。
改写的方法是三角恒等式 1 − cos x = 2sin²(x/2)：
// 这么写：两个几乎相等的数相减 (1 - Math.cos(x)) / (x * x) // 改成：没有减法了 2 * Math.sin(x / 2) ** 2 / (x * x)
右边那个式子里，sin(x/2) 本身就是个小数，直接算就好，不需要从 1 里减出来。同一个数学量，换一种代数写法，减法就消失了。这是卷 III 的全部内容。
灾难性相消（catastrophic cancellation）这个词里的「灾难性」，指的不是减法这个动作，而是后果——你从此再也不知道自己手上还剩几位是真的。
与之相对的叫良性相消（benign cancellation）：如果参与相减的两个数本身就是精确的（比如都是小整数，或者都是刚读进来的原始数据），那么相消之后什么也不会发生，结果照样精确。
所以判据不是「有没有相减」，而是：被减的这两个数，自己身上带着误差吗？带，就危险；不带，就没事。2.0 - 1.0 永远精确，Math.sqrt(2)**2 - 2 就不是。
相消最好的体感来自「把丢掉的位数打出来」。下面这段会打出一张表，每行是一个 x，最后一列是还剩几位可信：
import math
print(f"{'x':>8} {'直接相减':>24} {'改写后':>24} {'剩几位'}")
for k in range(2, 15):
x = 10.0 ** k
naive = math.sqrt(x + 1) - math.sqrt(x)
good = 1 / (math.sqrt(x + 1) + math.sqrt(x))
rel = abs(naive - good) / good
left = 16 if rel == 0 else max(0.0, -math.log10(rel))
print(f"1e{k:<6} {naive:>24.17g} {good:>24.17g} {left:.1f}")
会打出（节选）：
1e2 0.049875621120889946 0.049875621120890272 14.2 1e6 0.00049999987504634191 0.00049999987500006253 10.0 1e8 5.0000000555883162e-05 4.9999999874999996e-05 ★ 7.9 1e12 5.0000380724668503e-07 4.9999999999987493e-07 ★ 5.1 1e14 5.029141902923584e-08 4.9999999999999872e-08 ★ 2.2 1e16 0 5.0000000000000001e-09 ★ 0.0
每往右走一个数量级，就掉一位左右。到 10¹⁶ 直接归零——那时 √(x+1) 和 √x 存进 double 之后已经是同一个数了。
python3 -c "$(pbpaste)" ／ 或者存成 cancel.py 再 python3 cancel.py
在线跑：python.org/shell。JS 版把 math. 换成 Math. 就行，结果逐位相同。
GPS。卫星到接收机的距离是两万公里量级，而定位精度要做到米级——直接相减两个两万公里的数，米级的信息落在第 8 位上。所以差分 GPS 和 RTK 的做法从来不是「算得更准」，而是先减掉一个共同的基准，让剩下的量本身就是小数。这和上面 1 − cos x 的改写是同一个动作。
财务的损益。「本期利润 = 本期收入 − 本期成本」，当利润率只有 1% 时，你在用两个 100 倍大的数去算一个小数——丢两位。第 8 章会把这件事量化到一分钱。
物理模拟里的能量。检查「能量守恒」时算的是 E_now − E_start，两个几乎相等的大数。所以这类检查一定要用相对量，而且要清楚自己的判据本身能有几位可信。
回归模型的显著性。统计里 Var = E[x²] − (E[x])² 是同一个陷阱，第 10 章会看到它能把方差算成负数。
「精度不够就上更高精度。用 long double／Decimal／128 位，相消问题就解决了。」
高精度推迟问题，但不解决问题。相消丢掉的位数只跟「两个数前面有多少位相同」有关，跟你有多少位无关。从 53 位换到 113 位，你还是丢 8 位——只不过是从 113 位里丢，剩得多一点。
而且这个办法有个上限：如果两个量在数学上就相等（比如 x → 0 时的 1 − cos x），任何有限精度都会在某个 x 处归零，只是那个 x 更小一点。
真正解决它的是换一种写法——把减法从式子里拿掉。这不花一分钱，而且效果是无限的：改写后的 2sin²(x/2)/x² 在任意小的 x 上都给 0.5。
正确答案是 C：只对 8 位。直接算给出 0.00005000000055588316，正确答案是 0.000049999999874999996，相对误差 1.362 × 10⁻⁸。x 换成 10¹² 会更糟——丢 12 位。
x = 1e8 时太悲观了，还剩 8 位。但如果 x 到 10¹⁶，D 就成了正确答案：那时 √(x+1) 和 √x 存进 double 后完全相同，结果精确地等于 0。猜 D 的直觉方向是对的，只是提前了八个数量级。
相消的标准解法是有理化——分子分母同乘共轭式，把减法搬到分母上（分母上是相加，不会相消）：
Math.sqrt(x+1) - Math.sqrt(x) 在 x=1e8 处只剩 8 位
1 / (Math.sqrt(x+1) + Math.sqrt(x)) 在任何 x 上都有 16 位
推导只有一行：(√(x+1) − √x) × (√(x+1) + √x) = (x+1) − x = 1，所以原式等于 1 / (√(x+1) + √x)。数学上完全等价，数值上差了八个数量级。
同一招的其他形态：1 − cos x → 2sin²(x/2)；log(1+x) → log1p(x)；eˣ − 1 → expm1(x)。标准库里之所以有 log1p 和 expm1 这两个看起来多余的函数，原因就是这一章。
这一章的一句话
减法是浮点里最诚实的运算，它一个误差都不制造；它只是把结果缩小，于是把早就存在的误差顶到了前排——所以看到「两个相近的数相减」，要问的不是这一步准不准，而是这两个数各自是从哪里来的。
下一章看另一种放大方式：不靠相消，纯靠次数。累加一千万个 0.1，得到的不是 1 000 000，而是 999 999.999838975374586880207061767578125——差了 1 383 191 个 ulp。但这里有个反直觉的地方：这个误差不是按次数线性长的，它按 √n 长。而爱国者导弹那个误差偏偏就是线性长的。区别在哪？