4. 数值计算与浮点误差

在公式生成代码的场景里，最容易被忽视、也最容易造成“公式明明正确、结果却对不上”的原因，往往不是推导错误，而是数值计算本身。计算机用有限位二进制近似表示实数，任何一次加减乘除都可能引入微小偏差；当公式中存在相减抵消、大数吃小数、累加次数很多或矩阵接近奇异时，这些偏差会被放大，最终让校验失败。本节围绕浮点表示、舍入误差、数值稳定性、误差传播、条件数与病态问题，说明如何在生成代码时主动规避常见陷阱。

一、浮点数表示与舍入误差

以 IEEE 754 双精度（binary64）为例，一个数由符号位、11 位指数和 52 位尾数组成，规格化数的实际值为 (-1)^s × 1.f × 2^(e-1023)。尾数只有 52 位，意味着相邻可表示数之间的间隔（ulp）随数量级增大而增大。在 1 附近 ulp 约为 2^-52 ≈ 2.22e-16，而在 1e16 附近 ulp 已经达到 2，此时 1e16 + 1 在双精度下仍等于 1e16。

由此产生两类典型误差：一是表示误差，如 0.1 无法用二进制有限位精确表示，实际存储为略大于 0.1 的值；二是舍入误差，每次运算结果被舍入到最近的可表示数。单次误差很小，但误差会随运算累积和传播。代码校验时若直接用 == 比较两个浮点结果，几乎必然出现“应为 0.3、实为 0.30000000000000004”式的失败，正确做法是采用带容差的比较：

import math

def close(a, b, rtol=1e-9, atol=1e-12):
    # 相对+绝对混合容差，兼顾大数与接近零的情形
    return math.isclose(a, b, rel_tol=rtol, abs_tol=atol)

print(0.1 + 0.2 == 0.3)          # False
print(close(0.1 + 0.2, 0.3))     # True

容差不是随手取一个数。一般可按“结果量级 × 相对容差 + 绝对容差”设定，相对容差常取 1e-9 到 1e-12，绝对容差用于结果接近 0 的场景。若公式涉及多步迭代，容差还应随步数适当放宽，否则会把正常累积误差误判为错误。

二、数值稳定性与误差传播

数值稳定性指算法对输入扰动和舍入误差是否敏感。同一数学公式可以有多种等价写法，但计算稳定性差别很大。经典例子是二次方程求根：直接套用 (-b ± √(b²-4ac)) / (2a)，当 b² ≫ 4ac 时，-b + √(b²-4ac) 是两个相近数相减，发生“灾难性抵消”，有效数字大量丢失。稳定写法是先算不易抵消的那个根，再用韦达定理 x1·x2 = c/a 求另一个根：

import math

def roots(a, b, c):
    d = b*b - 4*a*c
    if d < 0:
        return None
    sq = math.sqrt(d)
    # 选取同号相加，避免相近数相减
    if b >= 0:
        x1 = (-b - sq) / (2*a)
    else:
        x1 = (-b + sq) / (2*a)
    x2 = c / (a * x1)   # 韦达定理，避免第二次抵消
    return x1, x2

误差传播分析关注“输入误差如何传到输出”。对 y = f(x1, x2, …)，一阶近似下输出误差约为各输入误差乘以偏导数的加权和。若某偏导数绝对值很大，说明该输入的小扰动会被显著放大。例如 y = 1/(x-1) 在 x 接近 1 时，x 的微小误差会让 y 剧烈变化。生成代码时，应识别这类敏感表达式，在文档或断言中标注输入的有效范围。

三、条件数与病态问题识别

条件数衡量问题本身对扰动的敏感程度，与算法无关。对线性方程组 Ax = b，条件数 cond(A) = ‖A‖·‖A⁻¹‖，条件数越大问题越病态。条件数为 10^k 时，输入相对误差约被放大 k 个数量级。判断是否病态，可计算条件数、观察矩阵行列式是否接近 0、检查主元是否过小。例如希尔伯特矩阵随阶数增加迅速病态，5 阶条件数已达约 4.8e5，10 阶超过 1e13，此时双精度求解结果基本不可信。

四、代码生成中规避常见数值陷阱

把上述分析落到生成代码的实践，可归纳为几条可操作的规则。第一，比较用容差，禁止对浮点结果直接用等号；断言中写明相对与绝对容差，并说明依据。第二，避免相近数相减，对形如 √(x+1)-√x 的表达式改写为 1/(√(x+1)+√x)，对 1-cos x 在小 x 时改用 2sin²(x/2)。第三，控制累加顺序，大量同号数求和时按绝对值从小到大累加，或使用 Kahan 补偿求和：

def kahan_sum(values):
    s = 0.0
    c = 0.0            # 补偿量
    for v in values:
        y = v - c
        t = s + y
        c = (t - s) - y
        s = t
    return s

第四，先缩放再计算，矩阵求解前做行归一化，避免主元量级悬殊；判断主元绝对值小于阈值时给出警告而非静默继续。第五，显式声明输入有效域，对 1/(x-1)、log(x)、√x 等在生成代码中加入范围检查，把“数学上成立、数值上不可靠”的输入挡在计算之前。第六，用高精度或符号计算做交叉校验，例如用 Python 的 decimal 或 fractions 对关键公式复算一遍，作为容差设定的参考基准。

综上，浮点误差无法消除，但可以被识别、量化和控制。生成代码时把“表示精度、运算稳定性、问题条件数、输入有效域”作为四个检查维度，配合容差比较与等价变形，才能让公式生成代码在真实数据上既算得对、也验得准。