深度长文ArticleSep 13, 20261 分钟 min read

「极限不存在」是假的：equation-solver 与 matrix 的七处数值失真

给 equation-solver 与 matrix 搭了一台探针，跑了七十余组数值对照。最危险的不是算法是阈值：一个 Richardson 外推权重把「极限不存在」印在了 x² 上，两处 1e-12 的绝对奇异阈值把 det = 1.008e+40 的矩阵判成不可逆。

我把 x^2 在 x = 1 处的极限提交给 equation-solver，它回答「极限不存在」。

我盯着这个答案看了好一会儿，第一反应是自己哪里写错了——x² 在 1 处的极限怎么会不存在。于是把输入框里的函数又打了一遍。不是打错了。工具真的认为 lim(x → 1) x² 不存在。

这篇文章想说的第一件事就是：数值工具的错误不会以报错的形式出现。它会以一个看起来完全合理的答案出现，合理到你会怀疑自己。

1. 探针怎么搭

src/tools/calculators.ts 有 4057 行，但整个模块只有一个导出：

export const CALCULATOR_TOOLS: ToolEntry[];

gaussSolve、lupDecompose、solveCubic、那四个极限求值器全是模块私有。想对它们灌任意输入，只能走构建期探针——把 config 从注册表里捞出来，直接调 compute()：

const eq = CALCULATOR_TOOLS.find((t) => t.slug === 'equation-solver')!;
// rows 里 label 是英文、labelZh 是产品文案——当证据只能引 label
const rows = eq.config!.compute(V({ type: 'limit', limExpr: 'x^2', limX0: '1', limDir: 'both' })).rows;

构建命令有个坑：rm -f 必须在 esbuild 之前。失败的 esbuild 会静默沿用上一次的成功产物，你以为在测新代码，其实在测旧的。

rm -f /tmp/probe.mjs && npx esbuild src/tools/.probe.ts --bundle --format=esm --outfile=/tmp/probe.mjs --log-level=error

七处失真和三组对照基准都是在这样一台探针上量出来的。探针和产物都留在 /tmp，不进仓库。

2. 「极限不存在」是假的

双重极限的误判，是这一轮最严重的发现。

type: 'limit' 且 limDir: 'both' 时，左右各跑一次单边逼近，各自外推到两个值，再拿差值卡一道 1e-7 的门槛：

const exists = limLeft !== null && limRight !== null && Math.abs(limLeft - limRight) < 1e-7;
// 差值过不了 1e-7，就判「极限不存在」

门槛的本意是对的——左右极限不相等，极限确实不存在。但它量到的不是数学，是算法自己留下的误差。

三组输入，三个假的「不存在」：

输入	Left	Right	真值
`x^2` @ 1	1.00000399997	0.999995999968	1
`x^2` @ 2	4.00000799997	3.99999199997	4
`(x^2-1)/(x-1)` @ 1	2.00000200003	1.99999800012	2

|L−R| 分别是 8.0000e-6、1.6000e-5、3.9999e-6，全超过 1e-7。

根因是一条四阶外推用错了步长比。 两个逼近器跑同一条步长序列 h = 1e-1 … 1e-6，相邻两步差 10 倍。但外推权是 (4, 1)/3——那是给 4 倍步长比设计的。比例写错了，O(h) 项从来没被消掉。

设 $v(h) = L + c\,h + \tfrac{1}{2}d\,h^2$ ，取 $h_1 = 10^{-5}$ 、 $h_2 = 10^{-6}$ ，代入：

\frac{4v_2 - v_1}{3} = L - 2\cdot10^{-6}\,c - 1.6\cdot10^{-11}\,d

左右两侧只差 $c$ 的符号，所以

|L_{\text{left}} - L_{\text{right}}| = 4\times10^{-6}\cdot|f'(x_0)|

用 x² 在 1 处验一下： $c = f'(1) = 2$ ， $d = 2$ ，右极限 $= 1 - 4\times10^{-6} - 3.2\times10^{-11} = 0.999995999968$ ，和工具印出来的 Right 逐位一致；左极限 $= 1.00000399997$ 。误差和公式对上了，说明外推器是在做减法时算错了比例，不是在随机漂移。

门槛于是变成了一道斜率闸门。 |L−R| < 1e-7 要求 |f′(x₀)| < 1/40。我用 s·x² 在 1 处扫了一组：

| 斜率 $s$ （ $f' = 2s$ ） | |L−R| | 判决 | |—|—|—| | 0.01 | 8.000e-8 | 存在 | | 0.012 | 9.600e-8 | 存在 | | 0.0125 | 1.000e-7 | 不存在 | | 0.013 | 1.040e-7 | 不存在 | | 0.05 | 4.000e-7 | 不存在 |

$s = 0.0125$ 时 |L−R| 恰好等于 1e-7，门槛是严格的 <，判死。边界正好落在 |f′| = 1/40。

sin(x)/x 在 0 处、x² 在 0 处都能正常返回——不是因为算得准，是因为它们在该点的导数恰好是 0，误差项系数为零，侥幸躲过。这跟「工具会判断极限是否存在」是两回事。

3. 一个 token 的修法

外推权改成 (10, 1)/9，就是给 10 倍步长比用的那一组：

const extrap = (10 * v2 - v1) / 9;  // 原来是 (4 * v2 - v1) / 3

残差从 $L - 2\times10^{-6}c$ 降到

\frac{10v_2 - v_1}{9} = L - 5\cdot10^{-12}\,d

O(h) 项没了，只剩 O(h²)，而且这一项在左右两侧同号，|L−R| 直接归零。改完把前面三组重跑：x^2 @1、x^2 @2、(x^2-1)/(x-1) @1 的 Left/Right/Two-sided 全变成真值，1 就是 1、4 就是 4、2 就是 2；0.0125·x² 从「不存在」翻成「存在」。

真正发散的 (2x+1)/(x-2) 在 2 处仍然判「不存在」（−5499998.00046 对 +5500001.99922），说明门槛没被废掉，只是终于量到了对的东西。

改动一共两处——:1687 和 :1787，一个是 calcInf，一个是 evalDirectional，同一个外推器被复制了两遍。这也是它一直没被发现的原因之一：修一处，另一处还在错。

4. `1/x²` 报出 1330000000000

同一个外推器，另一个失败方向。

1/x^2 在 0 处：Two-sided 1330000000000，Left 1330000000000，Right 1330000000000。三个值一模一样，门槛 |L−R| < 1e-7 轻松通过，工具宣称这个极限存在。

1/x 在 0 处就没这么好运——左右分别是 −1300000 和 +1300000，符号一翻，门槛拦住，判「不存在」。x²+1/x 同样判不存在。

区别只在奇偶性。奇次幂的奇点两侧异号，偶次幂的同号同向，而那道门槛只测差值不测量级，所以它拦得住符号翻转，拦不住一起爆炸。

把 §3 的修法套上，这个值从 1330000000000 变成 1110000000000——两侧仍然相等，门槛仍然通过。外推比例修好了，发散检测没修。

1330000000000 也不是什么新数字。工具对 x^2 → +∞ 印同样的值，因为那是同一条 1/t² 序列喂给同一个坏外推器。发散检测器和那条收敛慢的提示（:1705）只接在 dir === '±inf' 的分支上；有限点 x₀ 的分支既没有发散检测，也没有提示。

5. `auto` 静默线性化

方程组的 auto 模式把二次项直接扔掉了。

x^2 + y^2 = 1, x + y = 1 的真解是两个点：(0, 1) 和 (1, 0)。工具返回：

Identified Problem Type: 2x2 Linear System
Standard Matrix Form: 1x + 1y = 1;  1x + 1y = 1
Coefficient Determinant det(A): 0
System Solution: Infinitely many solutions (dependent equations)

它把两条方程都截成了 x + y = 1——x² 和 y² 没了，然后拿两条相同的方程去解，得到「无穷多解」。从它的视角看，推理是自洽的。

对照组证明线性分支本身没坏：x+y=1, 2x+2y=2 → det 0，无穷多解；x+y=1, 2x+2y=3 → det 0，No solution (parallel inconsistent lines)；2x - y = 3, x + 4y = 5 → det 9，x = 17/9、y = 7/9，都对。

检测非线性只需要看一件事：方程里有没有第二个自变量，输入里有没有逗号、分号或换行分隔。现在缺的就是这一层判断——auto 看到多个方程就当线性组处理，不管里面是不是二次。

同一条输出链上还有几处格式化毛病：2x + -1y = 3（负号跑到系数中间）、1x³ + -6x² + 11x + -6 = 0，以及零系数不消——1x³ + 0x² + -3x + 0 = 0。

6. 三次方程的根按 acos 顺序返回

solveCubic 用三角法解三实根情形，三个根是 cos(φ/3)、cos((φ+2π)/3)、cos((φ+4π)/3)——按相位顺序返回，不排序。

x³ − 6x² + 11x − 6 = 0  →  Root x1: 3 | Root x2: 1 | Root x3: 2    （真解 1, 2, 3）
x³ − 7x² + 14x − 8 = 0  →  Root x1: 4 | Root x2: 1 | Root x3: 2    （真解 1, 2, 4）

第三组更说明问题：x³ − 3x = 0 的三个根是 0 和 ±√3，工具返回 1.73205080757 | −1.73205080757 | −3.67394e-16——精确为零的那个根印成 −3.67394e-16，还排在最后。

我用系数路径（a:1, b:0, c:-3, d:0）重跑一遍，结果逐位一致，所以问题在 solveCubic 本身，不在 auto 的解析器。

修法是一行 roots.sort((a, b) => a - b)，但那个 −3.67394e-16 还得配一层 |x| < eps → 0 的清洗。否则排序之后它仍然排在最前面——只是从「最后一个不干净的根」变成「第一个不干净的根」。

7. Citardauq：打印的 Δ 不是用的 Δ

这一组是全篇唯一「算对了」的，但值得写出来——它说明稳定形式为什么是必需品。

a=1, b=1e8, c=1：

Discriminant Δ = b² − 4ac: 1e+16
Root x₁: -100000000
Root x₂: -1e-8
Parabola Vertex (xv, yv): (-50000000, -2499999999999999)

两个根都是对的，包括那个极小的 −1e-8。说明代码走的是稳定形式 $q = -\tfrac{1}{2}(b + \operatorname{sign}(b)\sqrt{\Delta})$ 、 $x_1 = q/a$ 、 $x_2 = c/q$ ，而不是教科书上的 $x = (-b \pm \sqrt{\Delta})/(2a)$ 。

有意思的是打印和计算用的不是同一个 Δ：formatNumber 把 9999999999999996 显示成 1e+16，但算术全程用的是完整双精度值。1e16 - 4 === 1e16 是 false，所以这次减法恰好代表得下——b 再大一点，就代表不下了。

拿教科书形式做对照，同一个输入：

b	小根（教科书形式）	真值	相对误差
1e8	−7.450580596923828e-9	−1e-8	−25.49%
1e10	0	−1e-10	−100.00%
1e12	0	−1e-12	−100.00%

b 到 1e10 就完全归零了——1e20 − 4 在双精度里就是 1e20，开方减 b 之后剩下 0。这不是精度损失，是灾难性抵消把答案整个吃掉。

formatNumber 还有个连带的小坑：Number.isSafeInteger(−2499999999999999) 为 true，所以顶点 y 坐标全展开打印；Number.isSafeInteger(1e16) 为 false，Δ 就走 1e12 以上的指数记法。同一条输出里两种风格混着，看起来像两套格式规则。

a = 0 的退化路径倒是干净的：a:0, b:2, c:4 → Linear Root x: -2，正确。

8. 两套阈值，同一个矩阵，双向分歧

这一轮第二个最贵的发现。

多项式回归在配齐次方程组时，解的是范德蒙型的正规矩阵 $A_{ij} = \sum x^{i+j}$ 。我在探针里把这个矩阵构造出来，用同样的数字分别喂给 polynomial-regression 和 matrix 工具。

第一组，x ∈ {1000, 2000, 3000, 4000, 5000}，3 次拟合：maxNorm = 2.0515e+22，gaussEps = max(1e-12, maxNorm·1e-13) = 2.0515e+9。

polynomial-regression：Result: — (singular system: the x values need more spread)
matrix 工具：Determinant det(A): 1.008e+40、Trace tr(A): 2.0515e+22，矩阵求逆正常完成

第二组，x ∈ {0, 0.02, 0.04, 0.06, 0.08, 0.1}，3 次拟合：maxNorm = 6.000e+0，gaussEps = 1.000e-12。

polynomial-regression：拟合成功，ŷ = 6.535748e-12x³ − 9.803622e-13x² + 1x − 1.516717e-16，R² = 1，RMSE = 1.74321e-16
matrix 工具：Determinant det(A): 0 + Inverse Matrix A⁻¹: Singular matrix (det ≈ 0, non-invertible)，但同一个矩阵的 Trace tr(A): 6.02215795296 照印无误

两个方向都出现过：同一个矩阵，一个说奇异一个说良态；一个说良态一个说奇异。

同一份文件里有三套阈值在打架：:61 的相对阈值 max(1e-12, maxNorm·1e-13)，和 :651、:674、:691 三处绝对阈值 1e-12。相对阈值随矩阵量级缩放，绝对阈值不动。量级一大，相对阈值把矩阵判死；量级一小，绝对阈值把矩阵判死。

第二组多项式回归的拟合结果本身也是垃圾——R² = 1、RMSE = 1.74e-16 看起来完美，但 6.535748e-12 是数值噪声，不是 3 次项系数。工具报告了成功，没有报告这个成功有多假。

9. 验过没问题的

不是所有东西都坏了。三个方法在对照基准上过关了，值得记下来，免得下一轮审查把它们也怀疑一遍。

经典 RK4，y′ = y + x、y(0) = 1、积到 x = 1，精确解 $2e - 2 = 3.43656365691809$ 。步数减半序列的误差比是 14.73 / 15.35 / 15.67 / 15.78 / 16.13——干净的 4 阶，一直顶到 N = 5 → 160。之后掉到 8.42、4.24，再往后翻转符号，从 N = 640 起误差钉在 +1.91e-12。所谓「验证到 1.9e-12」是准确的，但诚实的表述是：4 阶，直到撞上舍入地板。

反方向也对：从 (1, 2e−2) 往回积到 x = 0，N = 5000，返回 y(0) = 1，逐位精确。

五点中心差分（h = 1e-5）的七组对照，最大误差 2.1e-11，x² @ 1 是 0.0000e+0。这一支没问题。

自适应 Simpson：sin(x) 在 [0, π] 上误差恰好 0.0000e+0（值就是 2）；sqrt(x) 在 [0, 1] 上误差恶化到 −4.5007e-9，比其他几组差两个数量级——端点角奇点确实拖精度，但还在容差内。MAX_EVALS = 4096 是求值次数预算，不是递归深度上限；终止条件是 |delta| ≤ 15 × 1e-9 或预算耗尽，Richardson 校正取 delta / 15。

一个额外的坏消息：RK4 的步数字段会静默截断。odeSteps = 0 返回 3.43656365647，但同一行的说明印的是 h = 0.01, N = 100——用户要 0 步，得到 100 步。odeSteps = 999999 返回 3.43656365692，说明印 h = 0.0002, N = 5000——被截到上限。这个字段声明了 min: '1'、max: '5000'，但 compute() 里的守卫是静默 clamp，不是报错。两处截断都没提示。

10. 什么该修

按「错误答案会误导到什么程度」排：

:1687 和 :1787 的外推权——(4·v₂−v₁)/3 改成 (10·v₂−v₁)/9。一处 token，修掉一整类假的「极限不存在」。已验证。
有限点 x₀ 分支补发散检测——1/x² @ 0 现在报出 1330000000000 并宣称极限存在。外推比例修好也救不了它，需要独立的一层。
auto 模式判非线性——二次项被静默丢弃，x²+y²=1 被当成两条相同的直线。检测只需要「第二个自变量 + 分隔符」两个条件。
solveCubic 排序并清洗小根——roots.sort((a, b) => a - b) 加一层 |x| < eps → 0。
合并奇异判定——:61 的相对阈值和 :651、:674、:691 的绝对阈值统一口径。
RK4 步数守卫从静默 clamp 改成报错——要 0 步得到 100 步、要 999999 得到 5000，两处都没提示。

前五条是数值正确性，最后一条是交互诚实性。它们有一个共同点：没有一处会产生异常、崩溃或明显的失败信号。工具照常返回一行格式正常的结果，只是那个结果不对。

数值代码的类型检查器拦得住未定义变量，拦不住 Richardson 外推里的一个系数。这一轮的七处失真全都能通过 astro check，也不会让 E2E 变红——那些用例测的是交互和渲染，不测数值对不对。

矩阵计算器、方程求解器、多项式回归、微积分计算器。全部本地运算，数据不出浏览器。