4. 浮点加减法的误差传播

浮点加减法看起来只是“把两个数加起来”，但在 IEEE 754 双精度（binary64）下，一次加法要经过对阶、尾数相加、规格化、舍入四个阶段，每一步都可能引入误差。理解误差从哪里来，比记住“浮点不精确”这句结论重要得多。本节围绕三个关键机制展开：对齐移位造成的低位丢失、尾数截断/舍入造成的不可逆信息损失，以及灾难性抵消造成的有效位急剧衰减，并用可运行的代码演示“大数吃小数”和“抵消放大误差”两类典型现象。

一、对齐移位：指数差决定谁被“挤出去”

设两个双精度数 x、y，尾数各 53 位（含隐含位）。相加前必须让两者阶码一致，通常把阶码较小的数右移，使其指数等于较大者。右移的位数 d = |Ex − Ey|。被右移的数每右移一位，最低位就移出有效尾数范围，只能进入保护位、舍入位和粘滞位（guard/round/sticky）。当 d 超过约 53 时，较小数的全部有效位都被挤出，它对结果的贡献只剩“是否触发舍入”这一比特信息。

这就是“大数吃小数”的根源：不是加法算错了，而是小数的信息在对齐阶段就已经离开了可表示范围。下面用 Python 演示，所有输出均为示例环境下的合理结果，读者可在本机复现。

import math

a = 1.0e16          # 约 2^53.15
b = 1.0
print(a + b == a)   # True：b 被完全吃掉
print((a + b) - a)  # 0.0，b 的信息彻底丢失

# 连续累加：小数被反复吞掉
s = 0.0
for _ in range(10_000_000):
    s += 1.0e-8
print(s)            # 远小于 0.1，误差随规模累积

# 用 math.fsum 做精确求和作为对照
vals = [1.0e-8] * 10_000_000
print(math.fsum(vals))  # 接近 0.1

注意 a + b == a 为 True 并不代表程序有 bug，而是 binary64 在 1e16 附近的 ulp（单位末位）约为 2，1.0 小于半个 ulp，舍入后结果回到 a。判断“是否被吃掉”的实用标准是：|y| < 0.5·ulp(x) 时，x + y 在默认舍入模式下等于 x。

二、尾数截断与舍入：误差的第二次注入

对齐之后，两个 53 位尾数相加可能产生 54 位甚至更多位的结果，必须舍入回 53 位。IEEE 754 默认使用“就近舍入、偶数优先”（round-to-nearest, ties-to-even），它保证单次运算的相对误差不超过 0.5 ulp，即机器 epsilon 的一半（binary64 的 eps ≈ 2.22e-16）。但“单次误差有界”不等于“多次运算误差有界”：误差会随运算次数累积，也可能被后续的减法放大。

一个常被忽视的细节是双重舍入：在 x87 80 位扩展精度或某些融合乘加（FMA）场景中，中间结果先舍入到扩展精度、再舍入到目标精度，可能比单次舍入多出 1 ulp 的偏差。工程上应尽量让中间结果保持在目标精度，或使用 FMA 减少一次舍入。

三、灾难性抵消：减法把相对误差放大

当两个相近的数相减时，结果的绝对值很小，而参与运算的两个数各自的绝对误差并未减小，于是相对误差被放大。设 x、y 相近，各自有绝对误差 ε，则 x − y 的相对误差约为 ε/|x − y|。|x − y| 越小，放大越严重。这就是灾难性抵消（catastrophic cancellation）。

# 经典示例：用公式 (1 - cos x) / x^2 计算 x 很小时的函数值
import math

def f_naive(x):
    return (1.0 - math.cos(x)) / (x * x)

def f_stable(x):
    # 1 - cos x = 2 sin^2(x/2)，避免相近数相减
    s = math.sin(x / 2.0)
    return 2.0 * s * s / (x * x)

for x in (1e-2, 1e-4, 1e-6, 1e-8):
    print(f"{x:.0e}  naive={f_naive(x):.12f}  stable={f_stable(x):.12f}")

示例输出中，x = 1e-8 时 naive 版本会得到 0.0 或明显偏离 0.5 的值，而 stable 版本稳定趋近 0.5。原因很直接：cos(1e-8) 与 1.0 的差约为 5e-17，已经接近甚至低于 1.0 附近的 ulp（约 2.2e-16），相减后有效位几乎耗尽。改写为 2sin²(x/2) 后，不再有相近数相减，误差被控制在正常舍入水平。

四、误差传播的完整链路

把上述机制串起来，一次浮点加法从输入到输出的误差路径可以用下图表示。图中区分了“对齐移位导致的低位丢失”“舍入注入的 0.5 ulp 误差”和“抵消放大相对误差”三个环节，并标出可干预的位置。

五、工程上的应对策略

针对上述三类误差，工程实践中有几条成熟做法。第一，控制对齐移位的规模：累加大量量级差异大的数时，先按绝对值从小到大排序再累加，可显著减少“大数吃小数”。第二，使用补偿求和：Kahan 求和用一个补偿变量记录每次舍入丢失的低位，把误差重新注入后续累加，能把求和误差从 O(n·eps) 降到接近 O(eps)。第三，改写公式规避抵消：如上面的 1−cos x 改为 2sin²(x/2)，或用 math.expm1、math.log1p 等专门函数处理小量。第四，用高精度做参照：Python 的 decimal、fractions 或 math.fsum 可作为验证基准，帮助判断误差是否来自算法而非实现。

def kahan_sum(values):
    total = 0.0
    comp = 0.0          # 补偿项，保存被舍入丢失的低位
    for v in values:
        y = v - comp
        t = total + y
        comp = (t - total) - y
        total = t
    return total

vals = [1.0e-8] * 10_000_000
print(kahan_sum(vals))   # 明显优于朴素累加，接近 0.1

需要强调的是，这些手段解决的是误差传播问题，而不是让浮点变成精确算术。判断一段数值代码是否可靠，应关注三点：输入的有效位数量、中间是否出现相近数相减、累加规模与量级差异。把这三项作为代码审查的固定检查项，比事后调试更有效。