卷 III · 改写CH 09深度 9/20

求根公式是错的

卷 II 教你诊断，卷 III 开始动手改。规则很严格：不换语言、不加精度、不改数据、不换库，只把式子换一种在数学上完全等价的写法。第一个例子挑了个所有人都熟的——那条从初中背到现在的求根公式。它在一大类输入上会给出零位有效数字，而修法是改一行。

韦达定理数学等价 ≠ 数值等价log1p / expm1
▷ 猜一下差多少

解一个二次方程：

x² + 1 000 000 000·x + 1 = 0        （a=1, b=10⁹, c=1）

它的两个根很好估：大根约 −10⁹，小根约 −10⁻⁹（因为两根之积等于 c/a = 1）。

用求根公式算小根：

const d = Math.sqrt(b * b - 4 * a * c);
const x2 = (-b + d) / (2 * a);

问：x2 打出来是什么？

A −1e-9。公式是对的，double 有 16 位，绰绰有余 B −7.450580596923828e-9。大概对个一位 C 0 D NaN 或者 Infinity。开方或者除法出问题了

顺便想想：如果把 b 换成 10⁵（小很多），还有几位是对的？

把它拆开看

b² = 1000000000000000000            （10¹⁸）
ulp(10¹⁸) = 128                     ← 这一带的刻度是 128

b² − 4ac = 10¹⁸ − 4
         = 1000000000000000000      ← ⟨那个 4 被整个抹掉了⟩
                                       因为 4 比半格（64）还小

√(b²−4ac) = 1000000000             （正好 10⁹，一分不差）
−b        = −1000000000

−b + √(…) = 0                       ← 两个 10⁹ 相减，一位有效数字都不剩

x₂ = 0 / 2 = 0

答案是 C：0。

注意这里发生了两次抹零，而且是两种不同的抹零：

  1. 判别式里的 −4ac 被吃掉了。这是第 2 章：在 10¹⁸ 那一带，刻度是 128，一个 4 连半格都不到。
  2. 然后 −b + √(…) 相消。这是第 3 章：两个几乎相等的大数相减，把仅剩的信息也消光了。

第一次抹零决定了结果必然是 0——因为 √(b²) 精确等于 b，相减必得 0。那个 4 就是全部答案所在的地方，而它在第一步就没了。

它是怎么一路坏下去的

把 b 从小到大扫一遍，能看见有效位数是怎么一位一位掉的：

b公式版 x₂改写版 x₂相对误差还剩几位
10³−0.0010000010000226212−0.0010000010000022.1 × 10⁻¹¹10.7
10⁵−1.0000003385357559e-5−1.0000000001000001e-53.4 × 10⁻⁷6.5
10⁷−9.96515154838562e-8−1.00000000000001e-73.5 × 10⁻³2.5
10⁸−7.450580596923828e-9−1e-80.2550.6
10⁹0−1e-91.0000

b = 10³ 就已经掉了五位。不需要极端输入，一千就够。

这张表还顺带解释了为什么这个 bug 特别难发现：在小的 b 上它工作得挺好，测试用例基本都会过；坏得最厉害的地方（b 很大、c 很小）恰恰是「刚性」问题、病态系统、极端参数——也就是你最需要它对的时候。

改一行

修法用的是韦达定理，初中同一节课教的东西：两根之积等于 c/a。

// 这么写：两个根都用公式
const d  = Math.sqrt(b * b - 4 * a * c);
const x1 = (-b - d) / (2 * a);     // 这个是安全的：两个同号的数相加，没有相消
const x2 = (-b + d) / (2 * a);     // ★ 这个是灾难：相消

// 改成：只算安全的那个，另一个用乘积推出来
const d  = Math.sqrt(b * b - 4 * a * c);
const x1 = (-b - Math.sign(b) * d) / (2 * a);   // 永远取「加大不减小」的那一侧
const x2 = c / (a * x1);                        // ★ 韦达：x₁·x₂ = c/a

关键在于 x1 那一行：−b 和 ±d 符号相同，所以是两个同号的数相加——加法从不相消，永远安全。拿到精确的 x1 之后，x2 用一次除法就出来了，同样不相消。

整条路径上没有一处「两个相近的数相减」。

验算一下（把根代回原方程，看残差）：

公式版  x₂ = 0
        代回 x² + 10⁹x + 1  ⇒  1          ✗ 相对残差 100%

改写版  x₂ = −1e-9
        代回 x² + 10⁹x + 1  ⇒  0          ✓ 精确
◆ 主线

数学上等价，数值上不等价。

这是卷 III 的全部内容，也是这本书最实用的一句话。代数恒等式告诉你「这两个式子的值相同」，它没有告诉你「这两条计算路径的条件数相同」。

找改写机会的方法很机械，三步：

  1. 把式子里所有的减法圈出来。（包括伪装成加法的：加上一个负数）
  2. 对每一处问：这两个操作数会不会几乎相等？会 ⇒ 这里会相消。
  3. 用代数把这处减法搬走：有理化、恒等式、因式分解、换用同一族里没有减法的那个表达式。

三步做完通常就结束了。它不花运行时开销，只花一次思考。

标准库里那些「看起来多余」的函数

知道这一招之后，标准库里有几个函数会突然变得有意义。它们存在的唯一理由就是这一章：

看起来多余的函数它替掉的写法为什么
log1p(x)log(1 + x)x 很小时，1 + x 在存进 double 时就把 x 抹没了
expm1(x)exp(x) - 1x 很小时，exp(x) 约等于 1，减 1 是灾难性相消
hypot(x, y)sqrt(x*x + y*y)防的是另一种事：中间的平方会溢出／下溢
fma(a, b, c)a * b + c只舍入一次而不是两次。第 19 章会看到它带来的麻烦
logaddexp(a, b)log(exp(a) + exp(b))下一章的主角

实测 log1p：

xMath.log(1 + x)Math.log1p(x)
10⁻⁸9.999999889225291e-99.999999950000001e-9
10⁻¹²1.000088900581841e-129.999999999995e-13
10⁻¹⁶01e-16
10⁻¹⁷01e-17

x = 10⁻¹² 那一行值得看：log(1+x) 给出 1.000088900581841e-12，看起来非常像对的——量级对，前四位「1.000」也对。实际上从第五位起全是噪声。这类「看起来很对的错答案」，比返回 0 或 NaN 危险得多。

✎ 术语正名

数值等价这个词其实不该存在——两个数学上等价的表达式，在浮点里几乎从来不等价。

更准确的说法是：每个表达式对应一条计算路径，每条路径有自己的条件数。(a+b)² 和 a²+2ab+b² 是同一个函数的两条不同路径；(-b+√d)/2a 和 c/(a·x₁) 也是。

所以「优化表达式」这件事在数值计算里有两个互不相干的目标：更快（编译器干的）和更准（只有你能干）。而且它们经常冲突——第 19 章会看到，编译器为了更快做的表达式重排，正是把「更准」那一半破坏掉的东西。

⌨ 自己跑一遍

把这张表自己打一遍，尤其要看 b 从 10³ 到 10⁹ 那个滑坡：

import math
print(f"{'b':>8} {'公式版':>26} {'改写版':>26} {'相对误差':>11} {'剩几位'}")
for k in (3, 5, 7, 8, 9):
    b, a, c = 10.0 ** k, 1.0, 1.0
    d  = math.sqrt(b * b - 4 * a * c)
    bad  = (-b + d) / (2 * a)                 # 公式版：相消
    x1   = (-b - d) / (2 * a)                 # 安全的那个根
    good = c / (a * x1)                       # 韦达定理
    rel  = abs(bad - good) / abs(good)
    left = 16.0 if rel == 0 else max(0.0, -math.log10(rel))
    print(f'1e{k:<6} {bad:>26.17g} {good:>26.17g} {rel:>11.1e} {left:>6.1f}')

# 顺手看看判别式里那个 4 是怎么没的
b = 1e9
print('b*b       =', repr(b * b))
print('b*b - 4   =', repr(b * b - 4))          # 一模一样
print('ulp(1e18) =', math.ulp(1e18))           # 128

python3 quad.py

在线：python.org/shell。JS 版把 math. 换成 Math.，结果逐位相同（math.ulp 对应第 2 章那个位运算版本）。

▸ 在现实里

光线追踪里的射线—球相交就是解二次方程，而且恰恰落在最坏的参数区：射线原点离球很远（b 很大），球很小（c 相对很小）。用课本公式实现的求交函数，在远处会出现物体表面「闪烁」或者干脆穿透——这是渲染器里一类经典 bug，修法就是这一章这一行。Ray tracing 的经典教材（《Physically Based Rendering》）在讲二次方程时会专门用一节讲这件事。

金融里的隐含波动率、化学里的平衡浓度、弹道计算——凡是要解二次方程且系数量级悬殊的地方，都是同一个坑。

更普遍的一条：任何「小修正项」都在危险区。相对论修正、微扰展开、增量更新、差分——它们的共同结构是「一个大主项 + 一个小修正」，而你关心的往往正是那个小修正。只要它是从两个大数的差里出来的，它就已经没了。正确的做法永远是：直接为那个小量写一个式子，别让它作为差出现。

✗ 这个直觉是错的

「求根公式是数学定理，不可能有问题。有问题也是浮点数的锅，跟公式无关。」

公式当然没错——它在实数上永远成立。但你写下的不是公式，是一条计算路径。同一个定理有很多条路径，它们的数值表现天差地别。

这条直觉的危险之处在于它会让人停止寻找。「数学是对的，所以我这边没问题」——于是就去怪库、怪语言、怪硬件，或者干脆加个 if (Math.abs(x) < 1e-9) x = 0 把症状盖掉。

正确的反射是：看到答案不对，先在纸上把这个式子重写一遍，找那处减法。数值计算里绝大多数「玄学问题」，最后都是一处没被注意到的相消。

◇ 对账

正确答案是 C：x2 打出来是 0。b = 10⁵ 时还剩 6.5 位——已经丢掉了九位多。

A 「double 有 16 位，绰绰有余」——16 位是存储的位数，不是结果的位数。这一路上丢了两次：b²−4ac 里的 4 被抹掉（丢掉全部信息），−b+d 相消（把剩下的清零）。「位数够」和「算得出来」是两件事，这是整本书反复在说的话。 B 「−7.45e-9，对个一位」——这是 b = 10⁸ 时的答案，不是 10⁹ 时的。猜 B 说明直觉里已经有「会掉精度」这件事了，只是低估了它掉得多快：b 每涨十倍，就再掉一位；到 10⁹ 正好清零。这个「刚好归零」不是巧合——当 4ac < ulp(b²)/2 时，判别式必然精确等于 b²，结果必然精确等于 0。 D 「NaN 或 Infinity」——这个答案希望有个显式的失败信号。可惜没有。浮点数在这里表现得完美：开方合法，除法合法，结果是一个合法的有限数 0。没有异常、没有警告、没有 NaN。这本书里最危险的失败，全都长这样——安静、合法、看起来正常。这也是为什么必须靠事前的分析，而不能靠事后的报错。
◈ 换一种写法

这一章的改写模板，值得抄在手边：

x = (-b + Math.sqrt(b*b - 4*a*c)) / (2*a) const q = -(b + Math.sign(b) * Math.sqrt(b*b - 4*a*c)) / 2;
const x1 = q / a, x2 = c / q;　　—— 这是《Numerical Recipes》给的标准写法，两个根都安全

同一族的其他四条，一起记：

  • log(1+x) → log1p(x)；exp(x)-1 → expm1(x)
  • √(x+1) − √x → 1/(√(x+1) + √x)（有理化，第 3 章）
  • 1 − cos x → 2sin²(x/2)（三角恒等式）
  • (a−b)(a+b) ← 用它替 a² − b²，因为平方会先把两个数拉开量级

如果你只想记一条：任何时候你关心的是一个「小量」，就直接为那个小量写式子，绝不要让它以「两个大数之差」的形式出现。

这一章的一句话

你写下的不是数学公式，是一条计算路径；代数恒等式保证两条路径的终点相同，但完全不保证它们的路况相同。

下一章是同一招在统计上的应用，而且结果更荒唐：用教科书上那条一遍法的方差公式，去算四个数 10¹⁰+1, 10¹⁰+2, 10¹⁰+3, 10¹⁰+4 的方差。正确答案是 1.25。它算出来是 −16384。方差是一堆平方的平均，它在数学上不可能是负数。