fix(math variance): use a numerically stable two-pass formula for numbers (#19108)

## Description

`math variance` (and `math stddev`, which reuses it) computed the sum of
squared deviations with the single-pass identity `SS = sum(x²) −
(sum(x))²/n`. When the values are large but close together, the two
terms nearly cancel and `f64` loses the low-order digits, so the result
is wrong (0.0) or even **negative**, and `math stddev` returns **NaN**:

```nushell
[506250000.0 506250001.0] | math variance   # -32.0   (should be 0.25)
[506250000.0 506250001.0] | math stddev     # NaN
[1000000000 1000000001 1000000002] | math variance   # 0.0  (should be 0.667)
```

The duration/filesize path in the **same file** already used a
numerically stable two-pass form in `f64`. This change routes every
numeric input through `variance_unit_f64`, which now measures deviations
from a **shifted origin**: residuals `r = x - min`, their mean, then
`Sum(r - mean_r)^2 / denom`. Subtracting the smallest value is exact for
every input, so a constant sequence still comes out exactly 0, and
summing residuals cannot overflow where summing x does. The result type
is unchanged (a plain float, as before - integer inputs already produced
floats via the division).

## User-facing changes (Release notes)

- Fixed `math variance` and `math stddev` returning incorrect (zero or
negative) values, or `NaN`, for large values that are close together, by
using a numerically stable two-pass formula.
- `math variance` and `math stddev` no longer fail with
`OperatorOverflow` on integers whose square exceeds `i64` (roughly above
3×10⁹, depending on how many values are summed); they are computed as
`f64`. Integers above 2⁵³ are exact in `i64` but not in `f64`, so such
inputs now return an approximate result instead of an error — the same
precision trade-off `math avg` already makes.

## Additional notes

- The stable two-pass form was already present in the file for the
duration/filesize units; this routes the `Number` path through it too,
shifts the origin to the smallest value before averaging so the mean
cannot overflow, and removes the now-unused single-pass `sum_of_squares`
helper.
- Added a regression example for large, close values, plus nine
integration tests in `crates/nu-command/tests/commands/math/variance.rs`
and one in `stddev.rs`: the float case (which previously returned
`-32.0`, a negative variance), the integer case (previously `0.0`), the
`--sample` denominator, that integer input still yields `float`, a
`describe` check on `math stddev`, both sides of the integer boundary -
`[4000000000 4000000001]`, which used to fail with `OperatorOverflow`,
and `[9007199254740993 9007199254740994]`, which as `f64` lands on 2^53
and 2^53 + 2 and therefore has an exact variance of `1.0` - and three
that pin the shifted origin itself: `1..64 | each { 1e150 }` -> 0,
`[1e308 1e308]` -> 0, `[4503599627370496 4503599627370497]` -> 0.25.
- `Number` results now carry the input span instead of the call head,
matching the duration/filesize path already in this file.
- Checked locally against the pinned toolchain: `cargo fmt --all
--check`, `cargo clippy --workspace --all-targets --exclude nu-explore
--exclude nu_plugin_* -- -D warnings`, `cargo test -p nu-command --lib`
(822 cases), `cargo test -p nu-command --test tests --
commands::math::variance commands::math::stddev` (21 cases) and `typos`
all pass.
- This change was drafted with AI assistance (Claude); I verified the
reproduction and the fix against a local build.

crates/nu-command/src/math/variance.rs

@@ -1,6 +1,6 @@
 use crate::math::utils::{
-    NUMBER_INPUT_TYPES, NumericUnit, expand_range_input, run_with_function,
-    run_with_function_and_cell_paths, to_unit_f64, variance_denominator,
+    NumericUnit, expand_range_input, run_with_function, run_with_function_and_cell_paths,
+    to_unit_f64, variance_denominator,
 };
 use nu_engine::command_prelude::*;
 
@@ -141,42 +141,29 @@ impl Command for MathVariance {
                 example: "[1KB 3KB] | math variance",
                 result: Some(Value::test_float(1_000_000.0)),
             },
+            Example {
+                description: "Variance of large, close values is computed accurately.",
+                example: "[506250000.0 506250001.0] | math variance",
+                result: Some(Value::test_float(0.25)),
+            },
         ]
     }
 }
 
-fn sum_of_squares(values: &[Value], span: Span, head: Span) -> Result<Value, ShellError> {
-    let n = Value::int(values.len() as i64, span);
-    let mut sum_x = Value::int(0, span);
-    let mut sum_x2 = Value::int(0, span);
-    for value in values {
-        let v = match &value {
-            Value::Int { .. } | Value::Float { .. } => value.clone(),
-            Value::Error { error, .. } => return Err(*error.clone()),
-            other => {
-                return Err(ShellError::OnlySupportsThisInputType {
-                    exp_input_type: NUMBER_INPUT_TYPES.into(),
-                    wrong_type: other.get_type().to_string(),
-                    dst_span: head,
-                    src_span: other.span(),
-                });
-            }
-        };
-        let v_squared = &v.mul(span, &v, span)?;
-        sum_x2 = sum_x2.add(span, v_squared, span)?;
-        sum_x = sum_x.add(span, &v, span)?;
-    }
-
-    let sum_x_squared = sum_x.mul(span, &sum_x, span)?;
-    let sum_x_squared_div_n = sum_x_squared.div(span, &n, span)?;
-
-    let ss = sum_x2.sub(span, &sum_x_squared_div_n, span)?;
-
-    Ok(ss)
-}
-
-/// Variance for duration/filesize via `f64` units (ns / bytes) to avoid i64 overflow
-/// when squaring large values such as multi-second durations.
+/// Numerically-stable two-pass variance in `f64` base units.
+///
+/// All numeric types (int/float/duration/filesize) are converted to `f64` and
+/// the mean-deviation form `Σ(r − mean_r)² / denom` is used, where `denom` is
+/// `n` for the population variance and `n − 1` with `--sample`. The deviations
+/// are taken from residuals `r = x − min` rather than from `x` itself: shifting
+/// the origin to the smallest value keeps every residual of a constant sequence
+/// exactly `0`, and summing residuals cannot overflow the way summing `x` does
+/// for large-magnitude floats. This avoids the catastrophic cancellation of the
+/// single-pass `Σx² − (Σx)²/n` formula, which loses precision (and can even go
+/// negative) when values are large but close together. The result is a plain
+/// `f64`: squared base units (B² / ns²) for filesize/duration, or a plain number
+/// for int/float. Integer inputs above `2^53` are exact as `i64` but not as
+/// `f64`, so their variance is approximate rather than an error.
 fn variance_unit_f64(
     values: &[Value],
     sample: bool,
@@ -189,11 +176,15 @@ fn variance_unit_f64(
         nums.push(n);
     }
     let denom = variance_denominator(nums.len(), sample, head, span)? as f64;
-    let mean = nums.iter().sum::<f64>() / nums.len() as f64;
+    // Shift the origin to the smallest value, then average and accumulate the
+    // deviations of the residuals. The shift is exact for every input, so a
+    // constant sequence yields all-zero residuals and an exact variance of 0.
+    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 - mean;
+            let d = (x - origin) - mean;
             d * d
         })
         .sum::<f64>();
@@ -204,22 +195,12 @@ pub fn compute_variance(
     sample: bool,
 ) -> impl Fn(&[Value], Span, Span) -> Result<Value, ShellError> {
     move |values: &[Value], span: Span, head: Span| {
-        let unit = values_unit(values, head)?;
-        let denom = variance_denominator(values.len(), sample, head, span)?;
-        match unit {
-            // Duration/filesize: compute in f64 so large units (e.g. seconds as ns) don't overflow.
-            // Result is a plain number (squared units).
-            NumericUnit::Duration | NumericUnit::Filesize => {
-                let var = variance_unit_f64(values, sample, span, head)?;
-                Ok(Value::float(var, span))
-            }
-            NumericUnit::Number => {
-                // Arithmetic uses the original value span; errors point at the call head.
-                let ss = sum_of_squares(values, span, head)?;
-                let n = Value::int(denom as i64, head);
-                ss.div(head, &n, head)
-            }
-        }
+        // Validate that all values share one numeric unit (int/float may mix;
+        // duration/filesize must be homogeneous).
+        let _ = values_unit(values, head)?;
+        // Two-pass f64 variance for all units (numerically stable).
+        let var = variance_unit_f64(values, sample, span, head)?;
+        Ok(Value::float(var, span))
     }
 }
 

crates/nu-command/tests/commands/math/stddev.rs

@@ -28,3 +28,12 @@ fn stddev_duration_keeps_unit() -> Result {
         .run("[1sec 3sec] | math stddev")
         .expect_value_eq(expected)
 }
+
+#[test]
+fn stddev_large_close_values() -> Result {
+    // Regression: the single-pass formula produced a negative variance here,
+    // so the square root was `NaN`.
+    test()
+        .run("[506250000.0 506250001.0] | math stddev")
+        .expect_value_eq(0.5)
+}

crates/nu-command/tests/commands/math/variance.rs

@@ -54,3 +54,80 @@ fn variance_filesize_returns_number_in_bytes_squared() -> Result {
         .run("[1KB 3KB] | math variance")
         .expect_value_eq(1_000_000.0)
 }
+
+#[test]
+fn variance_large_close_floats() -> Result {
+    // Regression: the single-pass formula returned -32.0 here (negative variance).
+    test()
+        .run("[506250000.0 506250001.0] | math variance")
+        .expect_value_eq(0.25)
+}
+
+#[test]
+fn variance_large_close_ints() -> Result {
+    // Regression: the single-pass formula returned 0.0 here.
+    test()
+        .run("[1000000000 1000000001 1000000002] | math variance")
+        .expect_value_eq(2.0 / 3.0)
+}
+
+#[test]
+fn variance_of_equal_huge_floats_is_exact_zero() -> Result {
+    test()
+        .run("[1e308 1e308] | math variance")
+        .expect_value_eq(0.0)
+}
+
+#[test]
+fn variance_constant_huge_floats_is_exact_zero() -> Result {
+    // Plain summation overflows here; summing shifted residuals cannot, and every
+    // residual of a constant sequence is exactly 0.
+    test()
+        .run("1..64 | each { 1e150 } | math variance")
+        .expect_value_eq(0.0)
+}
+
+#[test]
+fn variance_large_close_ints_just_above_2_pow_52() -> Result {
+    // Exact as f64 (2^52 and 2^52 + 1 are both representable), so the shifted
+    // two-pass form gives the exact variance rather than a rounded approximation.
+    test()
+        .run("[4503599627370496 4503599627370497] | math variance")
+        .expect_value_eq(0.25)
+}
+
+#[test]
+fn sample_variance_large_close_floats() -> Result {
+    // Same input, `--sample` divides by n - 1.
+    test()
+        .run("[506250000.0 506250001.0] | math variance --sample")
+        .expect_value_eq(0.5)
+}
+
+#[test]
+fn variance_of_ints_is_still_float() -> Result {
+    // The two-pass rewrite must not change the result type.
+    test()
+        .run("[1 2 3 4 5] | math variance | describe")
+        .expect_value_eq("float")
+}
+
+#[test]
+fn variance_ints_whose_square_overflows_i64() -> Result {
+    // The single-pass path used to fail with OperatorOverflow here because it
+    // multiplied the values; they are summed as f64 now.
+    test()
+        .run("[4000000000 4000000001] | math variance")
+        .expect_value_eq(0.25)
+}
+
+#[test]
+fn variance_ints_above_2_pow_53_are_approximate() -> Result {
+    // Documented trade-off, not a regression: these are exact as i64 but not as
+    // f64, where they land on 2^53 and 2^53 + 2 — whose variance is exactly 1.0.
+    // So the answer is approximate for the original integers (their true variance
+    // is 0.25) rather than an error. `math avg` already behaves the same way.
+    test()
+        .run("[9007199254740993 9007199254740994] | math variance")
+        .expect_value_eq(1.0)
+}