修复前（251b3927 之前的 main）：
[506250000.0 506250001.0] | math variance # -32.0
[506250000.0 506250001.0] | math stddev # NaN
[1000000000 1000000001 1000000002] | math variance # 0.0
两个数相差 1，方差的真值是 0.25；三个连续整数，真值是 2/3。命令给出的却是一个负数、一个 NaN、一个 0。
方差是平方和，怎么会是负的？这正是有意思的地方：负号不是凭空冒出来的垃圾值，而是浮点减法在如实汇报它算出的东西——两个几乎相等的数相减，剩下的低位全是噪声，符号由舍入决定。math stddev 再对它开平方，噪声就变成了安静的 NaN。
这个修复 2026 年 10 月 10 日合入 nushell（#19108，squash 成 251b3927，里程碑 1.0）。值得说的不是"换了个更稳的公式"，而是选了哪一个：第一版改法修好了报告里的三条命令，却在常数序列上制造了新的错误；最后落地的版本，是在做两遍法之前先把原点平移到最小值。
结论先摆在这里：
- 出错的是单遍恒等式
S = Σx² − (Σx)²/n。数学上恒等，数值上却要把两个约 10¹⁷ 的数相减，才能得到 0.25；f64 在这个量级的 ulp 是 32，低位全是舍入噪声。结果可以是 0、任意正数，也可以是负数。 - 普通两遍法（先算均值，再算偏差平方和）修掉了消去，但均值本身要先求和：常数序列会被求和的舍入污染，
1e308会在求和时溢出。 - 先减去
min、在残差坐标里做两遍，两个问题一起消掉：常数序列的残差恰好全是 0；残差以range为上界，溢出的门槛从数据量级降到了range。但这不是免死金牌——平方和仍可能先溢出，即使真方差本身还能表示（第七节第 2 条）。
如果时间有限：第一至三节是主线，讲清问题、方案和取舍；第四、五节的完整实验与测试清单更像附录，可以回头再翻。
一、错不在推导，在消去
方差的定义式是两遍的：
S = Σ (xᵢ − x̄)², var = S / n （样本方差则除以 n−1）
要先知道均值，才能算偏差。为了只扫一遍数据，统计学教科书里给了一个恒等式：
S = Σxᵢ² − (Σxᵢ)²/n
推导没有问题，问题在计算：右边是两个几乎相等的巨大数相减。两项都约等于 n·a²（a 是数据的量级），要留下的 S 却只有 n·σ²；有效位数就在这一次减法里被吃掉了。
吃掉多少，取决于"数据离零有多远"相对于"数据有多散"。Chan、Golub 和 LeVeque 在 1983 年那篇综述里把这件事量化了：方差问题的条件数是
κ = ‖x‖₂ / √S ≥ 1
均值远大于标准差时 κ ≈ x̄/σ。教科书单遍算法的相对误差界是 N·κ²·u（u 是机器精度，f64 约 1.1×10⁻¹⁶）。也就是说，问题的敏感性按"均值/标准差"的平方增长。这是数据本身的条件数，换算法只能改变误差界里的常数和阶。
拿报告里的第一条命令走一遍：
5.0625×10⁸的平方是2.562890625×10¹⁷。这个量级的 f64 相邻可表示数相隔 32（ulp = 2⁵），一次平方最多带来 16 的舍入误差。- 两个平方相加得
5.1×10¹⁷，这里 ulp 已经是 64；Σx²与(Σx)²/n各自带着几十的绝对误差。 - 真方差是 0.25。两个带几十误差的数相减，得到的 −32.0 从头到尾都是噪声，符号也是任意的。
三个整数的例子更直白：Σx² 与 (Σx)²/n 都约 3×10¹⁸，这个量级的 ulp 是 512；真 S = 2，两边转成 f64 再相除，算出来是 0.0。
负方差在数学上不可能出现，所以它其实是个很好的报警信号。Chan 等人的原话是 "a blessing in disguise since this at least alerts the programmer that disastrous cancellation has occurred"——可惜 math stddev 把它变成了一个安静的 NaN。
二、把原点挪到数据旁边
合入的代码是这样的（variance.rs）：
let denom = variance_denominator(nums.len(), sample, head, span)? as f64;
// Shift the origin to the smallest value, then average and accumulate the
// deviations of the residuals.
let origin = nums.iter().copied().fold(f64::INFINITY, f64::min);
let mean = nums.iter().map(|x| x - origin).sum::<f64>() / nums.len() as f64;
let ss = nums.iter().map(|x| {
let d = (x - origin) - mean;
d * d
}).sum::<f64>();
Ok(ss / denom)
也就是 r = x − min、m = Σr/n、var = Σ(r − m)² / denom。实现上是三遍扫描：取 min、求残差均值、累加平方和。
这个平移在数学上是恒等的：Σ(x−d − (x̄−d))² = Σ(x−x̄)²，减去同一个常数不改变方差。它在浮点里换来四件事：
- 条件数与数据的绝对大小脱钩。 平移后
κ_d = ‖x−d‖₂/√S。Chan 等人给了个很干净的界：只要d落在[min, max]内，就有κ_d² < 1 + N。取d = min是其中最省事的一种：它是数据点本身，残差非负，且以range为上界。数据的量级a从公式里消失了。 - 常数序列精确为 0。 每个残差都是
x − x = 0，均值 0，平方和 0；与a是 0.3 还是 1e150 无关。 - 两遍求和都没有相消。 残差非负，求均值那一遍是正数相加；第二遍加的是平方，更是如此。误差按项数线性累积，不再被放大
κ²倍。 - 溢出的门槛从
a降到range。 残差以range为上界，Σr ≤ n·range，而不是n·a。但这只是把门槛降低，不是免死金牌：平方和要装进 f64，而ss ≈ n·var，即使var本身可表示，某个d²或ss仍可能先溢出。range超过约1.3×10¹⁵⁴、或者n·range²逼近 f64 上限，就要当心（第七节第 2 条）。
为什么是 min，不是 mean
"先减去均值再做两遍"看起来更自然，代价是均值本身要先求和：
[1e308 1e308]的和是inf，均值inf，偏差−inf，结果是NaN；- 常数序列的和带着
~n·u·a的绝对误差，均值偏离真值~u·a，残差不再是 0，平方和后放大成~(n·u·a)²。
Chan 等人也提醒过这一点：拿算出来的均值当平移量，条件数的界会退化成 1 + N²κ²u²，而且"算出的均值甚至可能落在 min 和 max 之外"（"One can easily construct examples where the computed mean does not even lie between min xᵢ and max xᵢ"）。min 不需要先算任何东西，它是输入里现成的、精确的值。
选 min 而不是首元素，是为了残差非负、上界最小；原点的选取也不依赖元素顺序。注意这只说原点：求和仍是顺序累加，换个顺序结果会在末位差几个 ulp（第三节 Welford 那条有例子）。取任意一个数据点作平移，κ² < 1+N 同样成立。
三、为什么不是"随便哪个两遍法"
这个 PR 的第一版，只是把 Number 路径接上同文件里 duration/filesize 已在用的普通两遍法。报告里的三条命令确实修好了，但 review 中发现它弄坏了一个原本正确的场景：常数序列。
对完全相同的值，单遍恒等式常常精确抵消成 0（两项由同一批数算出，容易落到同一个浮点数上）；普通两遍法的均值被求和舍入污染，反而给出非零。1..64 | each { 1e150 } | math variance 在 main 上是 0，在第一版里是 2.97×10²⁶⁹；0.3 重复 10 次：main 上是 0，第一版是 3.08×10⁻³³。修复不能用一个失败模式换另一个。
方案最后由维护者 fdncred 在 review 里拍板：r = x − min、m = Σr/n、Σ(r − m)²/denom，并把常数序列的用例补成测试。
于是把候选方案都实测了一遍。下表是 variance_demo 的输出（population variance，f64；粗体为与真值不符；×k 表示重复 k 次）：
| 用例 | 真值 | 单遍 | 两遍 | Σ(x/n) | corrected | Welford | 平移 |
|---|---|---|---|---|---|---|---|
[506250000.0 506250001.0] | 0.25 | −32 | 0.25 | 0.25 | 0.25 | 0.25 | 0.25 |
[1000000000 1000000001 1000000002] | 2/3 | 0 | 2/3 | 2/3 | 2/3 | 2/3 | 2/3 |
1..64 | each { 1e150 } | 0 | 0 | 2.97e269 | 2.97e269 | 0 | 0 | 0 |
[0.3]×10 | 0 | 0 | 3.08e−33 | 3.08e−33 | 0 | 0 | 0 |
[2⁶²]×7 | 0 | 0 | 0 | 1048576 | 0 | 0 | 0 |
[1e308 1e308] | 0 | NaN | inf | 0 | NaN | 0 | 0 |
[2⁵² 2⁵²+1] | 0.25 | 4.5e15 | 0.5 | 0.5 | 0.25 | 0.5 | 0.25 |
[2⁵³+1 2⁵³+2] | 1（f64） | 1.8e16 | 2 | 2 | 1 | 2 | 1 |
[4000000000 4000000001] | 0.25 | 0 | 0.25 | 0.25 | 0.25 | 0.25 | 0.25 |
[1e300 1e300+1e285] | 2.5e569（超出 f64） | NaN | inf | inf | NaN | inf | inf |
[−1e308 1e308 1e308 1e308] | 7.5e615（超出 f64） | NaN | inf | inf | NaN | NaN | NaN |
表里
[4000000000 4000000001]一行是 f64 算术的复刻：nushell 的旧整数路径在Value::mul做平方时就会OperatorOverflow，根本走不到最后的减法。
逐条看，每个变体都有自己的坑：
- 普通两遍：消去问题解决了，但常数序列上出现了回归（
1e150重复 64 次、0.3重复 10 次），均值在[1e308 1e308]上溢出；[2⁵², 2⁵²+1]也算成 0.5——两数之和2⁵³+1舍入到2⁵³，均值成了2⁵²，第二个偏差被算成 1。 Σ(x/n)：把均值改成逐项除以 n 再求和，[1e308 1e308]从inf变成 0。但x/n不再精确，常数序列的偏差带上了同号误差，1e150重复 64 次仍然错，[2⁶²]×7给出 1048576。- corrected two-pass（
Σd² − (Σd)²/n，d = x − x̄）：很接近正确答案。它修好了报告里的用例、常数序列和[2⁵², 2⁵²+1]；但修正项用到的还是普通均值，[1e308 1e308]仍是NaN。它靠"第二个误差抵消第一个"来处理常数序列，平移则是从结构上让残差为 0。 - Welford：单遍、稳定、常见，但依赖输入顺序。
[2⁵², 2⁵²+1]给出 0.5，换成[2⁵²+1, 2⁵²]给出 0：后者的均值从2⁵²+1减 0.5，恰好落在两个可表示数中间，按偶舍入回到2⁵²，第二个偏差成了 0。同一个多重集，换个顺序结果不同——对 shell 命令来说，这是最难解释的一类行为。平移法并没有消除顺序依赖，只是把差别压到末位：[0, 1e16, 1, 1]与[0, 1, 1, 1e16]相差 1 个 ulp（1.875×10³¹对1.8749999999999998×10³¹），不会把 0.5 变成 0。 - 平移原点：以上全部命中。
还有一个容易忽略的对比：单遍公式在常数序列上反而常常是精确的。这是它唯一"赢"的场景，也是为什么不能简单换成普通两遍法——修复必须同时满足"大偏移不消去"和"常数序列仍为 0"，平移 min 是同时做到这两点的最短路径。
四、实测：误差怎么随量级长起来
构造一组数据：a, a+1, …, a+63，真方差固定为 (64²−1)/12 = 341.25。让 a 从 10⁰ 走到 10¹⁵，每个数量级取 16 个尾数（a = 10^t·(1 + j/16)），一共 256 组样本，比较单遍公式与平移两遍法的相对误差。
| 量级 a | 单遍相对误差（16 个样本的中位数） | 恰好为 0 的样本 | 为负的样本 |
|---|---|---|---|
| ≤ 10⁶ | 0（精确） | 0 | 0 |
| 10⁷ | 4.6×10⁻⁴ | 0 | 0 |
| 10⁸ | 8.1×10⁻³ | 0 | 0 |
| 10⁹ | 0.56 | 6/16 | 0 |
| 10¹⁰ | 24 | 8/16 | 3/16 |
| 10¹¹ | 1.2×10⁴ | 4/16 | 6/16 |
| 10¹² | 7.9×10⁵ | 6/16 | 7/16 |
| 10¹³ | 1.0×10⁸ | 2/16 | 7/16 |
| 10¹⁴ | 1.6×10¹⁰ | 0 | 5/16 |
| 10¹⁵ | 1.6×10¹² | 5/16 | 4/16 |
256 组样本里，单遍公式精确的只有 116 组（全部落在 1.2×10⁷ 以下），31 组恰好为 0，32 组为负，最大相对误差 4.9×10¹²；平移两遍法的最大相对误差是 0。
误差不是"偶尔差一点"。从 10⁹ 起，每升一个数量级，中位误差大致放大两个数量级，和 κ² 的二次增长对得上：到 10⁹ 已经是 56%，到 10¹⁰ 大多数结果要么是 0、要么是负的。
拐点为什么在 10⁷ 附近
这一小节解释的是上面这组整数样本，不是一般数据的精确性定理。
Σx² 是 64 个平方的累加，每个约 a²。在这张网格上（t ≥ 4 时 a 是整数），只要累加过程不越过 2⁵³，每一步都是精确整数：64a² < 2⁵³ ≈ 9.0×10¹⁵ 给出 a 约小于 1.186×10⁷。实测拐点在 1.19×10⁷ 附近，比这条线略高：2⁵³ 到 2⁵⁴ 之间 ulp 是 2，整数仍落在格点上，个别样本碰巧没有舍入。116 组精确样本里，112 组在 10⁷ 以下，另外 4 组落在 10⁷ 到 1.19×10⁷ 之间。
减数 (Σx)²/64 在这一段也精确，但原因不是大小：a 为整数时 Σx = 64a + 2016 是 32 的倍数，平方是 1024 的倍数，而该量级的 ulp 最多 128，所以平方本身没有舍入。两个中间项都精确，差值才是精确的 21840。
换一组一般浮点数据，这两条都不自动成立：平方和会舍入，(Σx)² 也会舍入，误差界仍是 N·κ²·u，随均值/标准差的平方增长。这条线解释的是实验现象，不是普遍保证。
三个边界
[2⁵², 2⁵²+1]→ 0.25。 两个数都能被 f64 精确表示，平移后残差是 0 和 1，均值 0.5，方差 0.25，全程精确。测试里专门钉了这一条——它是"平移让减法变精确"的最小证明。[2⁵³+1, 2⁵³+2]→ 1.0。 作为i64它们是精确的，转成f64后落到2⁵³和2⁵³+2，真方差 1.0；原始整数的真方差是 0.25。这是转f64的取舍，不是回归，release note 和代码注释里都写明了（math avg早就在做同样的取舍）。[4000000000, 4000000001]→ 0.25。 旧路径用Value::mul（checked_mul）算平方，1.6×10¹⁹超过i64，直接报OperatorOverflow；现在按 f64 算。这是一处用户可见的行为变化：从"报错"变成"给出近似值"。
救不回来的是真值本身就超出 f64 的输入：[1e300, 1e300+1e285] 的真方差约 2.5×10⁵⁶⁹，结果只能是 inf；[−1e308, 1e308, 1e308, 1e308] 的真值约 7.5×10⁶¹⁵，NaN 或 inf 都可以接受。这一点在 PR 讨论里明确划成了"不追"。
还有一类更容易被忽略：真值可表示，中间平方先溢出。99 个 0 加一个 1e155，真方差 9.9×10³⁰⁷ 在 f64 内，但最大偏差 9.9×10¹⁵⁴ 的平方是 9.8×10³⁰⁹，平移后仍返回 inf；详见第七节第 2 条。
五、十个测试钉住了什么
修复带了 9 个 math variance 集成测试和 1 个 math stddev 测试：
| 测试 | 钉住的行为 |
|---|---|
variance_large_close_floats | 开篇的 −32.0 → 0.25 |
variance_large_close_ints | 曾返回 0.0 的三个大整数 → 2/3 |
variance_of_equal_huge_floats_is_exact_zero | [1e308 1e308] → 0（平移前是 NaN） |
variance_constant_huge_floats_is_exact_zero | 1..64 | each { 1e150 } → 0（普通两遍法的回归点） |
variance_large_close_ints_just_above_2_pow_52 | [2⁵², 2⁵²+1] → 0.25（减法精确） |
variance_ints_whose_square_overflows_i64 | [4000000000, 4000000001] → 0.25（旧路径报溢出） |
variance_ints_above_2_pow_53_are_approximate | [2⁵³+1, 2⁵³+2] → 1.0，并注明原整数真值 0.25 |
sample_variance_large_close_floats | --sample 的 n−1 分母：0.5 |
variance_of_ints_is_still_float | 结果类型仍是 float |
stddev_large_close_values | 同一输入的 math stddev → 0.5（曾是 NaN） |
浮点测试的写法有两点值得说：
- 期望值全部选可精确表示的数：0.25、0.5、0、1.0。这些用例的正确结果本身就是精确的，所以用精确相等来断言，不用 epsilon。用
1e-9之类的容差，会把真正的回归一起放过去。 2/3这个用例写成2.0 / 3.0，和代码里的ss / denom走同一条舍入路径——两边都是把 2 除以 3 一次，结果必然相同。这不是"近似相等"，是同一个运算。
另外，边界的两侧都钉了：[4000000000, 4000000001] 是旧路径报错、新路径精确的一侧；[2⁵³+1, 2⁵³+2] 是转 f64 后近似的一侧。只钉一侧的话，边界移动了也不容易被发现；两侧都钉，边界才算被围住。
六、别人的方差是怎么算的
这不是 nushell 独有的难题，也不是我们发明的方法。几个主流实现的取舍：
| 实现 | 公式 | 数值手段 |
|---|---|---|
CPython statistics | 同一个单遍恒等式 Σx² − (Σx)²/n | 用 Fraction 做精确有理数运算，最后才转 float；源码注释直言 "This formula has poor numeric properties for floats, but with fractions it is exact." |
NumPy _var | 两遍：先均值，再 (x − mean)² 求和 | 归约用 pairwise summation |
R var/cov | 两遍：先均值，再中心化叉积 | 均值再做一遍修正（tmp + Σ(x−tmp)/n），累加器用 long double |
Julia var | 数组：两遍（mean + centralize_sumabs2）；迭代器：Welford | 数组走 pairwise 归约 |
nushell math variance | 平移原点 + 两遍 | 全程 f64，无额外依赖 |
共同点是：没有一个是"单遍恒等式 + 顺序求和 + f64"。CPython 保留了那个公式，但换成了精确算术；R 用 long double 和修正均值；NumPy、Julia 用两遍加 pairwise。nushell 的约束不同：数据已经是内存里的 &[Value]，两遍不额外占空间；Rust 没有 long double；为一个 math 命令引入有理数或补偿求和，不划算。在这个约束下，平移 min 加 f64 两遍，是最小改动。
七、边界、措辞与存疑
- "平移对每个输入都精确"这句话要打个折。 合并后的注释写的是 "The shift is exact for every input"。严格说，
x − min只在min > 0且x ≤ 2·min时由 Sterbenz 引理保证无舍入；一般情况下残差会有一次舍入，误差以range为尺度。真正被用到的两条性质不受影响：常数序列的残差恰为 0；残差以range为上界，所以κ不随a增长。这句注释值得在后续 PR 里改准确——它不影响结果，但会影响下一个读代码的人。 - 平方和仍会溢出，即使真方差可表示。 99 个
0加一个1e155：真方差9.9×10³⁰⁷在 f64 内，平移后最大偏差9.9×10¹⁵⁴，平方溢出，结果是inf。ss ≈ n·var，所以var得小到约MAX/n才稳。这是所有用 f64 先累加平方和、再除以 n 的实现共同的上限；#19108 只是把触发条件从"数据离零远"改成"range大"。 - 真值超出 f64 的输入不做承诺。
1e300与1e300+1e285的方差是2.5×10⁵⁶⁹，f64 装不下；inf或NaN都是合理结果。 2⁵³以上的整数是近似的。 与math avg一致，已写进 release note。- N 很大时误差界不是零。 按 Chan 等人 Table 1，两遍法的界是
N·u + N²·κ²·u²；平移到min之后代入κ_d² < 1 + N，第二项变成N³·u²。N = 10⁷ 时约1.2×10⁻¹¹，还小于第一项N·u ≈ 1.1×10⁻⁹；元素数到一亿上下，第二项就追平第一项了。本文的用例都在 64 个元素以内。 - 单遍公式在常数序列上常常是精确的，这是它最后一块阵地；任何替换方案，都得先证明自己在这块阵地上不倒退。
复现
验证工程在 variance_demo/，只依赖标准库：
cd variance_demo
cargo test --release # 7 个测试：合并后的行为、各候选方案的反例、两个边界
cargo run --release # 打印文中全部表格与扫描数据，并生成 out/*.csv
naive / two_pass_plain / mean_div_n / corrected_two_pass / welford / shifted 六个函数分别对应文中的六个变体，naive 的运算顺序与旧的 sum_of_squares 一致。文中 256 组扫描数据由 cargo run 生成，图由 plot_fig1.py 从 CSV 画出。
nushell 侧可以用这三条命令验证（需要包含 251b3927 的构建）：
[506250000.0 506250001.0] | math variance # 0.25
[506250000.0 506250001.0] | math stddev # 0.5
[1000000000 1000000001 1000000002] | math variance # 0.6666666666666666
参考
- Tony F. Chan, Gene H. Golub, Randall J. LeVeque, Algorithms for Computing the Sample Variance: Analysis and Recommendations, The American Statistician 37(3), 1983. doi:10.1080/00031305.1983.10483115【正文读过；条件数、Table 1 的误差界、§3 的平移分析均出自此文】
- B. P. Welford, Note on a Method for Calculating Corrected Sums of Squares and Products, Technometrics 4(3), 1962. doi:10.1080/00401706.1962.10490022【仅经由 Chan 等人的综述引用；本文没有单独复现其误差证明】
- Sterbenz 引理（1974）：若
a ≤ b ≤ 2a，则b − a在同精度浮点中精确。Wikipedia【教科书结论】 - nushell PR #19108；合并提交
251b3927；variance.rs；variance 测试 - CPython
Lib/statistics.py（_ss）；NumPy_methods.py（_var）；Rcov.c；JuliaStatistics.jl。以上源码均为 2026-10 快照。
脚注
① 方法说明：文中的 nushell 输出（−32.0、NaN、0.0、0.25、0.5、0.666…）取自 PR 里对本地构建的验证记录；数值实验由 variance_demo 在 Rust 1.98 / f64 上重跑，运算顺序与旧实现一致。所有表格数字来自 cargo run --release 的实际输出，没有手算。
② 成本：10M 个 f64 的纯循环，单遍 5.2–5.9 ms，平移两遍 16–18 ms（20 次运行取均值；i9-13900HX，rustc 1.98.1，--release，数据为 5.0625×10⁸ + i mod 1000）。平移在实现上是三遍扫描（取 min、残差均值、平方和），多出来的是常数遍数的内存访问，复杂度仍是 O(n)。在 nushell 里每个元素还要先过一次 to_unit_f64，这笔固定开销两种算法都要付。
③ 修订中修正过四处说法：其一，最初把"普通两遍法会弄坏常数序列"归因于"均值不准"，更准确的说法是——均值的求和舍入误差以数据量级 a 为尺度，残差因此不再为 0；平移把这一步消掉之后，常数序列的 0 才是结构性的。其二，§7 的误差界早先写作 κ²u + N²κ²u²，对照 Table 1 第 3 行，正确的是 N·u + N²·κ²·u²。其三，§2 早先说"只有 range 接近 f64 上限才会溢出"，实际 99 个 0 加一个 1e155 这类输入真方差可表示、平方项先溢出。其四，§4 早先把拐点归因于"两边中间结果都精确"，实际被减数先越过 2⁵³，减数是靠整除结构才保持精确。
