diff --git a/crates/functions/src/stat.rs b/crates/functions/src/stat.rs
index 554b75271..270bd2373 100644
--- a/crates/functions/src/stat.rs
+++ b/crates/functions/src/stat.rs
@@ -7,1593 +7,28 @@ use crate::special::{erf, gamma, invert, ln_beta, ln_gamma, norm_cdf, norm_inv,
 use crate::util::{A, R, S, arg, array_numbers, array_val, as_array, boolean, has, num, num_val, numbers, numbers_a, opt_num};
 use crate::{Arg, Ctx, FnSpec, VAR};
 
+mod descriptive;
+mod distributions;
+mod regression;
+
+pub use descriptive::{aggregate_values, aggregate_values_k};
+use descriptive::{
+    arg_nums, arg_pairs, avedev, correl, count_impl, counta_impl, countblank, devsq, frequency, geomean, harmean, intercept, kurt, large, max_of,
+    mean, median, min_of, mode_mult, mode_numbers, mode_sngl, percentile_exc, percentile_inc, percentrank_impl, prob, quartile_exc, quartile_inc,
+    rank_impl, skew, skew_p, slope, small, standardize, trimmean, var_p, var_s,
+};
+use distributions::{
+    beta_dist, beta_inv, betadist_legacy, binom_dist, binom_dist_range, binom_inv, chisq_dist, chisq_dist_rt, chisq_inv, chisq_inv_rt, chisq_test,
+    confidence_norm, confidence_t, expon_dist, f_dist, f_dist_rt, f_inv, f_inv_rt, f_test, fisher, gamma_dist, gamma_fn, gamma_inv, gammaln, gauss,
+    hypgeom_dist, lognorm_dist, lognorm_inv, negbinom_dist, norm_dist, norm_inv_fn, norm_s_dist, norm_s_inv, poisson_dist, t_dist, t_dist_2t,
+    t_dist_rt, t_inv, t_inv_2t, t_test, tdist_legacy, weibull_dist, z_test,
+};
+#[cfg(test)]
+use regression::{MAX_REGRESSORS, check_regression_size};
+use regression::{covariance, forecast, linest_impl, steyx, trend_impl};
+
 type C<'a> = &'a mut dyn Ctx;
 
-// ---------------------------------------------------------------------------------------------
-// Basic helpers over number lists
-// ---------------------------------------------------------------------------------------------
-
-fn mean(x: &[f64]) -> R<f64> {
-    if x.is_empty() {
-        return Err(CellError::Div0);
-    }
-    Ok(x.iter().sum::<f64>() / x.len() as f64)
-}
-
-fn devsq(x: &[f64]) -> f64 {
-    if x.is_empty() {
-        return 0.0;
-    }
-    let m = x.iter().sum::<f64>() / x.len() as f64;
-    x.iter().map(|v| (v - m) * (v - m)).sum()
-}
-
-fn var_s(x: &[f64]) -> R<f64> {
-    if x.len() < 2 {
-        return Err(CellError::Div0);
-    }
-    Ok(devsq(x) / (x.len() - 1) as f64)
-}
-
-fn var_p(x: &[f64]) -> R<f64> {
-    if x.is_empty() {
-        return Err(CellError::Div0);
-    }
-    Ok(devsq(x) / x.len() as f64)
-}
-
-fn sorted(mut v: Vec<f64>) -> Vec<f64> {
-    v.sort_by(f64::total_cmp);
-    v
-}
-
-fn median(v: Vec<f64>) -> R<f64> {
-    if v.is_empty() {
-        return Err(CellError::Num);
-    }
-    let s = sorted(v);
-    let n = s.len();
-    let hi = s.get(n / 2).copied().unwrap_or(0.0);
-    if n % 2 == 1 { Ok(hi) } else { Ok(s.get(n / 2 - 1).copied().unwrap_or(0.0).midpoint(hi)) }
-}
-
-/// Values occurring at least twice, ordered by first occurrence, with their (equal, maximal) count.
-fn modes(v: &[f64]) -> Vec<f64> {
-    let mut counts: Vec<(f64, usize)> = Vec::new();
-    for &x in v {
-        if let Some(e) = counts.iter_mut().find(|e| e.0 == x) {
-            e.1 += 1;
-        } else {
-            counts.push((x, 1));
-        }
-    }
-    let best = counts.iter().map(|e| e.1).max().unwrap_or(0);
-    if best < 2 {
-        return Vec::new();
-    }
-    counts.into_iter().filter(|e| e.1 == best).map(|e| e.0).collect()
-}
-
-fn mode_sngl(v: &[f64]) -> R<f64> {
-    modes(v).first().copied().ok_or(CellError::NA)
-}
-
-/// Numbers for MODE / MODE.SNGL / MODE.MULT. Inside multi-cell ranges and arrays only numbers
-/// count, like [`numbers`]. Scalar arguments differ from SUM: in Excel 365, text and logical
-/// values are skipped whether typed or referenced (even numeric text such as `"1"`, so
-/// `MODE(1,"1",2,3)` is `#N/A`, and `MODE(1,TRUE,TRUE,2)` is `#N/A`), but a single-cell reference
-/// to a blank cell is `#VALUE!` (`MODE(B9,2,D7)` with D7 blank). Arguments are processed left to
-/// right, so the first error wins. The results cached in POI's FormulaEvalTestData come from an
-/// older Excel that answered `#VALUE!` for typed text; they are not the reference here.
-fn mode_numbers(args: &[Arg]) -> R<Vec<f64>> {
-    let mut out = Vec::new();
-    for a in args {
-        match &a.value {
-            Value::Array(arr) => {
-                for v in arr.iter() {
-                    match v {
-                        Value::Number(n) => out.push(*n),
-                        Value::Error(e) => return Err(*e),
-                        _ => {}
-                    }
-                }
-            }
-            Value::Number(n) => out.push(*n),
-            Value::Error(e) => return Err(*e),
-            Value::Text(_) | Value::Bool(_) => {}
-            Value::Empty => return Err(CellError::Value),
-        }
-    }
-    Ok(out)
-}
-
-fn product(v: &[f64]) -> f64 {
-    if v.is_empty() { 0.0 } else { v.iter().product() }
-}
-
-fn min_of(v: &[f64]) -> f64 {
-    if v.is_empty() { 0.0 } else { v.iter().copied().fold(f64::INFINITY, f64::min) }
-}
-
-fn max_of(v: &[f64]) -> f64 {
-    if v.is_empty() { 0.0 } else { v.iter().copied().fold(f64::NEG_INFINITY, f64::max) }
-}
-
-/// k-th largest (k is rounded up like Excel).
-fn kth(v: Vec<f64>, k: f64, largest: bool) -> R<f64> {
-    let k = k.ceil();
-    if k.is_nan() || k < 1.0 || k > v.len() as f64 {
-        return Err(CellError::Num);
-    }
-    let s = sorted(v);
-    let k = k as usize;
-    let idx = if largest { s.len() - k } else { k - 1 };
-    s.get(idx).copied().ok_or(CellError::Num)
-}
-
-fn interp(s: &[f64], pos: f64) -> R<f64> {
-    let i = pos.floor();
-    let frac = pos - i;
-    let i = i as usize;
-    let lo = *s.get(i).ok_or(CellError::Num)?;
-    match s.get(i + 1) {
-        Some(hi) if frac > 0.0 => Ok(lo + frac * (hi - lo)),
-        _ => Ok(lo),
-    }
-}
-
-fn percentile_inc(v: Vec<f64>, k: f64) -> R<f64> {
-    if v.is_empty() || !(0.0..=1.0).contains(&k) {
-        return Err(CellError::Num);
-    }
-    let s = sorted(v);
-    interp(&s, k * (s.len() - 1) as f64)
-}
-
-fn percentile_exc(v: Vec<f64>, k: f64) -> R<f64> {
-    let n = v.len() as f64;
-    if v.is_empty() || !(k > 0.0 && k < 1.0) {
-        return Err(CellError::Num);
-    }
-    let pos = k * (n + 1.0);
-    if pos < 1.0 || pos > n {
-        return Err(CellError::Num);
-    }
-    let s = sorted(v);
-    interp(&s, pos - 1.0)
-}
-
-fn quartile_inc(v: Vec<f64>, q: f64) -> R<f64> {
-    let q = q.trunc();
-    if !(0.0..=4.0).contains(&q) {
-        return Err(CellError::Num);
-    }
-    percentile_inc(v, q / 4.0)
-}
-
-fn quartile_exc(v: Vec<f64>, q: f64) -> R<f64> {
-    let q = q.trunc();
-    if !(1.0..=3.0).contains(&q) {
-        return Err(CellError::Num);
-    }
-    percentile_exc(v, q / 4.0)
-}
-
-fn arg_nums(a: &[Arg], i: usize) -> R<Vec<f64>> {
-    array_numbers(&arg(a, i)?.value)
-}
-
-/// Pairs of numbers at matching positions of two arrays (pairs with a non-number skipped).
-/// Different sizes give `#N/A`; errors propagate.
-fn pairs(x: &Value, y: &Value) -> R<(Vec<f64>, Vec<f64>)> {
-    let xa = as_array(x);
-    let ya = as_array(y);
-    if xa.data.len() != ya.data.len() {
-        return Err(CellError::NA);
-    }
-    let mut xs = Vec::new();
-    let mut ys = Vec::new();
-    for (a, b) in xa.data.iter().zip(ya.data.iter()) {
-        if let Value::Error(e) = a {
-            return Err(*e);
-        }
-        if let Value::Error(e) = b {
-            return Err(*e);
-        }
-        if let (Value::Number(p), Value::Number(q)) = (a, b) {
-            xs.push(*p);
-            ys.push(*q);
-        }
-    }
-    Ok((xs, ys))
-}
-
-fn arg_pairs(a: &[Arg], i: usize, j: usize) -> R<(Vec<f64>, Vec<f64>)> {
-    pairs(&arg(a, i)?.value, &arg(a, j)?.value)
-}
-
-/// (sum of cross deviations, devsq x, devsq y, mean x, mean y).
-fn moments(x: &[f64], y: &[f64]) -> R<(f64, f64, f64, f64, f64)> {
-    let mx = mean(x)?;
-    let my = mean(y)?;
-    let mut sxy = 0.0;
-    let mut sxx = 0.0;
-    let mut syy = 0.0;
-    for (a, b) in x.iter().zip(y) {
-        sxy += (a - mx) * (b - my);
-        sxx += (a - mx) * (a - mx);
-        syy += (b - my) * (b - my);
-    }
-    Ok((sxy, sxx, syy, mx, my))
-}
-
-fn correl(x: &[f64], y: &[f64]) -> R<f64> {
-    if x.len() < 2 {
-        return Err(CellError::Div0);
-    }
-    let (sxy, sxx, syy, _, _) = moments(x, y)?;
-    if sxx == 0.0 || syy == 0.0 {
-        return Err(CellError::Div0);
-    }
-    Ok(sxy / (sxx * syy).sqrt())
-}
-
-fn slope(y: &[f64], x: &[f64]) -> R<f64> {
-    let (sxy, sxx, _, _, _) = moments(x, y)?;
-    if sxx == 0.0 {
-        return Err(CellError::Div0);
-    }
-    Ok(sxy / sxx)
-}
-
-fn intercept(y: &[f64], x: &[f64]) -> R<f64> {
-    let (sxy, sxx, _, mx, my) = moments(x, y)?;
-    if sxx == 0.0 {
-        return Err(CellError::Div0);
-    }
-    Ok(my - sxy / sxx * mx)
-}
-
-fn trunc_digits(r: f64, sig: f64) -> f64 {
-    let f = 10f64.powi(sig as i32);
-    if !f.is_finite() {
-        return r;
-    }
-    ((r * f) + 1e-9).floor() / f
-}
-
-fn col_value(v: Vec<f64>) -> R<Value> {
-    let n = v.len();
-    array_val(n, 1, v.into_iter().map(Value::Number).collect())
-}
-
-// ---------------------------------------------------------------------------------------------
-// Aggregates
-// ---------------------------------------------------------------------------------------------
-
-fn count_impl(a: &[Arg]) -> f64 {
-    let mut c = 0usize;
-    for x in a {
-        match &x.value {
-            Value::Array(arr) => c += arr.iter().filter(|v| v.is_number()).count(),
-            Value::Number(_) => c += 1,
-            _ if x.from_ref => {}
-            Value::Bool(_) => c += 1,
-            Value::Text(t) if gridcraft_core::parse::parse_number_text(t).is_some() => c += 1,
-            _ => {}
-        }
-    }
-    c as f64
-}
-
-fn counta_impl(a: &[Arg]) -> f64 {
-    let mut c = 0usize;
-    for x in a {
-        match &x.value {
-            Value::Array(arr) => c += arr.iter().filter(|v| !v.is_empty()).count(),
-            Value::Empty if x.from_ref => {}
-            _ => c += 1,
-        }
-    }
-    c as f64
-}
-
-fn countblank(a: &[Arg], _: C) -> R<Value> {
-    let v = &arg(a, 0)?.value;
-    let blank = |x: &Value| matches!(x, Value::Empty) || matches!(x, Value::Text(t) if t.is_empty());
-    let c = match v {
-        Value::Array(arr) => arr.iter().filter(|x| blank(x)).count(),
-        x => usize::from(blank(x)),
-    };
-    num_val(c as f64)
-}
-
-fn mode_mult(a: &[Arg], _: C) -> R<Value> {
-    let v = mode_numbers(a)?;
-    let m = modes(&v);
-    if m.is_empty() {
-        return Err(CellError::NA);
-    }
-    col_value(m)
-}
-
-fn large(a: &[Arg], _: C) -> R<Value> {
-    num_val(kth(arg_nums(a, 0)?, num(a, 1)?, true)?)
-}
-
-fn small(a: &[Arg], _: C) -> R<Value> {
-    num_val(kth(arg_nums(a, 0)?, num(a, 1)?, false)?)
-}
-
-fn rank_impl(a: &[Arg], avg: bool) -> R<Value> {
-    let x = num(a, 0)?;
-    let v = arg_nums(a, 1)?;
-    let asc = opt_num(a, 2, 0.0)? != 0.0;
-    let eq = v.iter().filter(|&&y| y == x).count();
-    if eq == 0 {
-        return Err(CellError::NA);
-    }
-    let before = v.iter().filter(|&&y| if asc { y < x } else { y > x }).count();
-    let r = before as f64 + 1.0;
-    num_val(if avg { r + (eq as f64 - 1.0) / 2.0 } else { r })
-}
-
-fn percentrank_impl(a: &[Arg], exc: bool) -> R<Value> {
-    let s = sorted(arg_nums(a, 0)?);
-    let x = num(a, 1)?;
-    let sig = opt_num(a, 2, 3.0)?.trunc();
-    if sig < 1.0 || s.is_empty() {
-        return Err(CellError::Num);
-    }
-    let (first, last) = (s.first().copied().unwrap_or(0.0), s.last().copied().unwrap_or(0.0));
-    if x < first || x > last {
-        return Err(CellError::NA);
-    }
-    let n = s.len() as f64;
-    let rank_of = |v: f64| -> f64 {
-        let less = s.iter().filter(|&&y| y < v).count() as f64;
-        if exc {
-            (less + 1.0) / (n + 1.0)
-        } else if n <= 1.0 {
-            1.0
-        } else {
-            less / (n - 1.0)
-        }
-    };
-    let r = if s.contains(&x) {
-        rank_of(x)
-    } else {
-        let lo = s.iter().copied().filter(|&y| y < x).fold(f64::NEG_INFINITY, f64::max);
-        let hi = s.iter().copied().filter(|&y| y > x).fold(f64::INFINITY, f64::min);
-        let (rl, rh) = (rank_of(lo), rank_of(hi));
-        rl + (x - lo) / (hi - lo) * (rh - rl)
-    };
-    num_val(trunc_digits(r, sig))
-}
-
-fn trimmean(a: &[Arg], _: C) -> R<Value> {
-    let v = arg_nums(a, 0)?;
-    let p = num(a, 1)?;
-    if !(0.0..1.0).contains(&p) || v.is_empty() {
-        return Err(CellError::Num);
-    }
-    let n = v.len();
-    let per_side = ((n as f64 * p) / 2.0).floor() as usize;
-    let s = sorted(v);
-    let kept = s.get(per_side..n - per_side).unwrap_or(&[]);
-    num_val(mean(kept)?)
-}
-
-fn kurt(v: &[f64]) -> R<f64> {
-    let n = v.len() as f64;
-    if v.len() < 4 {
-        return Err(CellError::Div0);
-    }
-    let m = mean(v)?;
-    let s = var_s(v)?.sqrt();
-    if s == 0.0 {
-        return Err(CellError::Div0);
-    }
-    let sum4: f64 = v.iter().map(|x| ((x - m) / s).powi(4)).sum();
-    Ok(n * (n + 1.0) / ((n - 1.0) * (n - 2.0) * (n - 3.0)) * sum4 - 3.0 * (n - 1.0).powi(2) / ((n - 2.0) * (n - 3.0)))
-}
-
-fn skew(v: &[f64]) -> R<f64> {
-    let n = v.len() as f64;
-    if v.len() < 3 {
-        return Err(CellError::Div0);
-    }
-    let m = mean(v)?;
-    let s = var_s(v)?.sqrt();
-    if s == 0.0 {
-        return Err(CellError::Div0);
-    }
-    let sum3: f64 = v.iter().map(|x| ((x - m) / s).powi(3)).sum();
-    Ok(n / ((n - 1.0) * (n - 2.0)) * sum3)
-}
-
-fn skew_p(v: &[f64]) -> R<f64> {
-    let m = mean(v)?;
-    let s = var_p(v)?.sqrt();
-    if s == 0.0 {
-        return Err(CellError::Div0);
-    }
-    Ok(v.iter().map(|x| ((x - m) / s).powi(3)).sum::<f64>() / v.len() as f64)
-}
-
-fn geomean(v: &[f64]) -> R<f64> {
-    if v.is_empty() || v.iter().any(|&x| x <= 0.0) {
-        return Err(CellError::Num);
-    }
-    Ok((v.iter().map(|x| x.ln()).sum::<f64>() / v.len() as f64).exp())
-}
-
-fn harmean(v: &[f64]) -> R<f64> {
-    if v.is_empty() || v.iter().any(|&x| x <= 0.0) {
-        return Err(CellError::Num);
-    }
-    Ok(v.len() as f64 / v.iter().map(|x| 1.0 / x).sum::<f64>())
-}
-
-fn avedev(v: &[f64]) -> R<f64> {
-    if v.is_empty() {
-        return Err(CellError::Num);
-    }
-    let m = mean(v)?;
-    Ok(v.iter().map(|x| (x - m).abs()).sum::<f64>() / v.len() as f64)
-}
-
-fn frequency(a: &[Arg], _: C) -> R<Value> {
-    let data = arg_nums(a, 0)?;
-    let bins = arg_nums(a, 1)?;
-    if bins.is_empty() {
-        return col_value(vec![data.len() as f64]);
-    }
-    let mut order: Vec<usize> = (0..bins.len()).collect();
-    order.sort_by(|&i, &j| bins.get(i).copied().unwrap_or(0.0).total_cmp(&bins.get(j).copied().unwrap_or(0.0)));
-    let mut counts = vec![0.0; bins.len() + 1];
-    for x in data {
-        let slot = order.iter().find(|&&i| x <= bins.get(i).copied().unwrap_or(f64::INFINITY)).copied().unwrap_or(bins.len());
-        if let Some(c) = counts.get_mut(slot) {
-            *c += 1.0;
-        }
-    }
-    col_value(counts)
-}
-
-fn prob(a: &[Arg], _: C) -> R<Value> {
-    let (xs, ps) = arg_pairs(a, 0, 1)?;
-    let lo = num(a, 2)?;
-    let hi = if has(a, 3) { num(a, 3)? } else { lo };
-    if ps.iter().any(|&p| p <= 0.0 || p > 1.0) {
-        return Err(CellError::Num);
-    }
-    let total: f64 = ps.iter().sum();
-    if (total - 1.0).abs() > 1e-7 {
-        return Err(CellError::Num);
-    }
-    num_val(xs.iter().zip(&ps).filter(|(x, _)| **x >= lo && **x <= hi).map(|(_, p)| p).sum())
-}
-
-fn standardize(a: &[Arg], _: C) -> R<Value> {
-    let (x, m, sd) = (num(a, 0)?, num(a, 1)?, num(a, 2)?);
-    if sd <= 0.0 {
-        return Err(CellError::Num);
-    }
-    num_val((x - m) / sd)
-}
-
-fn forecast(a: &[Arg], _: C) -> R<Value> {
-    let x = num(a, 0)?;
-    let (ys, xs) = arg_pairs(a, 1, 2)?;
-    if ys.is_empty() {
-        return Err(CellError::Div0);
-    }
-    num_val(intercept(&ys, &xs)? + slope(&ys, &xs)? * x)
-}
-
-fn steyx(a: &[Arg], _: C) -> R<Value> {
-    let (ys, xs) = arg_pairs(a, 0, 1)?;
-    let n = ys.len() as f64;
-    if ys.len() < 3 {
-        return Err(CellError::Div0);
-    }
-    let (sxy, sxx, syy, _, _) = moments(&xs, &ys)?;
-    if sxx == 0.0 {
-        return Err(CellError::Div0);
-    }
-    num_val(((syy - sxy * sxy / sxx) / (n - 2.0)).max(0.0).sqrt())
-}
-
-fn covariance(a: &[Arg], sample: bool) -> R<Value> {
-    let (xs, ys) = arg_pairs(a, 0, 1)?;
-    let n = xs.len();
-    if n == 0 || (sample && n < 2) {
-        return Err(CellError::Div0);
-    }
-    let (sxy, ..) = moments(&xs, &ys)?;
-    num_val(sxy / if sample { (n - 1) as f64 } else { n as f64 })
-}
-
-// ---------------------------------------------------------------------------------------------
-// Regression (LINEST, LOGEST, TREND, GROWTH)
-// ---------------------------------------------------------------------------------------------
-
-struct Xy {
-    y: Vec<f64>,
-    /// n rows of k predictor values.
-    x: Vec<Vec<f64>>,
-    k: usize,
-    y_is_row: bool,
-    y_rows: usize,
-    y_cols: usize,
-}
-
-fn strict_num(v: &Value) -> R<f64> {
-    match v {
-        Value::Number(n) => Ok(*n),
-        Value::Error(e) => Err(*e),
-        _ => Err(CellError::Value),
-    }
-}
-
-/// Most predictor columns LINEST, LOGEST, TREND and GROWTH accept. The normal equations are a
-/// k×k matrix inverted in O(k³), so an unbounded k (one `SEQUENCE(1,16384)` argument) allocated
-/// gigabytes and could abort the process.
-const MAX_REGRESSORS: usize = 255;
-/// Most multiply-adds (rows × k²) a regression may spend building its normal equations.
-const MAX_REGRESSION_WORK: usize = 1_000_000_000;
-
-/// Rejects regressions too large to fit: more than [`MAX_REGRESSORS`] predictors is `#VALUE!`,
-/// more than [`MAX_REGRESSION_WORK`] multiply-adds is `#NUM!`.
-fn check_regression_size(n: usize, k: usize) -> R<()> {
-    if k > MAX_REGRESSORS {
-        return Err(CellError::Value);
-    }
-    match n.checked_mul(k).and_then(|nk| nk.checked_mul(k)) {
-        Some(work) if work <= MAX_REGRESSION_WORK => Ok(()),
-        _ => Err(CellError::Num),
-    }
-}
-
-fn parse_xy(y: &Value, x: Option<&Value>) -> R<Xy> {
-    let ya = as_array(y);
-    let n = ya.data.len();
-    let y_is_row = ya.rows == 1 && ya.cols > 1;
-    let yv: Vec<f64> = ya.data.iter().map(strict_num).collect::<R<_>>()?;
-    let (xs, k) = match x {
-        None => ((1..=n).map(|i| vec![i as f64]).collect(), 1),
-        Some(xv) => {
-            let xa = as_array(xv);
-            if xa.rows == ya.rows && xa.cols == ya.cols {
-                (xa.data.iter().map(|v| strict_num(v).map(|f| vec![f])).collect::<R<Vec<_>>>()?, 1)
-            } else if !y_is_row && xa.rows == n {
-                check_regression_size(n, xa.cols)?;
-                let mut rows = Vec::with_capacity(n);
-                for r in 0..n {
-                    rows.push((0..xa.cols).map(|c| strict_num(xa.get(r, c).unwrap_or(&Value::Empty))).collect::<R<Vec<_>>>()?);
-                }
-                (rows, xa.cols)
-            } else if y_is_row && xa.cols == n {
-                check_regression_size(n, xa.rows)?;
-                let mut rows = Vec::with_capacity(n);
-                for c in 0..n {
-                    rows.push((0..xa.rows).map(|r| strict_num(xa.get(r, c).unwrap_or(&Value::Empty))).collect::<R<Vec<_>>>()?);
-                }
-                (rows, xa.rows)
-            } else {
-                return Err(CellError::Ref);
-            }
-        }
-    };
-    Ok(Xy { y: yv, x: xs, k, y_is_row, y_rows: ya.rows, y_cols: ya.cols })
-}
-
-struct Fit {
-    coef: Vec<f64>,
-    b: f64,
-    se: Vec<f64>,
-    se_b: Option<f64>,
-    r2: R<f64>,
-    sey: R<f64>,
-    f: R<f64>,
-    df: f64,
-    ssreg: f64,
-    ssresid: f64,
-}
-
-/// Inverts a symmetric positive semi-definite matrix by Gauss–Jordan on the diagonal, dropping
-/// (zeroing) collinear variables. Returns the inverse and the dropped flags.
-fn sym_inverse(m: &[Vec<f64>]) -> (Vec<Vec<f64>>, Vec<bool>) {
-    let p = m.len();
-    let mut a: Vec<Vec<f64>> = m.to_vec();
-    let mut inv: Vec<Vec<f64>> = (0..p).map(|i| (0..p).map(|j| if i == j { 1.0 } else { 0.0 }).collect()).collect();
-    let mut dropped = vec![false; p];
-    let diag: Vec<f64> = (0..p).map(|i| m.get(i).and_then(|r| r.get(i)).copied().unwrap_or(0.0).abs()).collect();
-    for j in 0..p {
-        let piv = a.get(j).and_then(|r| r.get(j)).copied().unwrap_or(0.0);
-        let scale = diag.get(j).copied().unwrap_or(0.0);
-        if piv.abs() <= 1e-12 * scale || piv.abs() < 1e-300 {
-            if let Some(d) = dropped.get_mut(j) {
-                *d = true;
-            }
-            for row in a.iter_mut() {
-                if let Some(v) = row.get_mut(j) {
-                    *v = 0.0;
-                }
-            }
-            if let Some(r) = a.get_mut(j) {
-                r.iter_mut().for_each(|v| *v = 0.0);
-            }
-            if let Some(r) = inv.get_mut(j) {
-                r.iter_mut().for_each(|v| *v = 0.0);
-            }
-            continue;
-        }
-        let arow: Vec<f64> = a.get(j).map(|r| r.iter().map(|v| v / piv).collect()).unwrap_or_default();
-        let irow: Vec<f64> = inv.get(j).map(|r| r.iter().map(|v| v / piv).collect()).unwrap_or_default();
-        for i in 0..p {
-            if i == j {
-                continue;
-            }
-            let f = a.get(i).and_then(|r| r.get(j)).copied().unwrap_or(0.0);
-            if f == 0.0 {
-                continue;
-            }
-            if let Some(r) = a.get_mut(i) {
-                r.iter_mut().zip(&arow).for_each(|(v, w)| *v -= f * w);
-            }
-            if let Some(r) = inv.get_mut(i) {
-                r.iter_mut().zip(&irow).for_each(|(v, w)| *v -= f * w);
-            }
-        }
-        if let Some(r) = a.get_mut(j) {
-            *r = arow;
-        }
-        if let Some(r) = inv.get_mut(j) {
-            *r = irow;
-        }
-    }
-    (inv, dropped)
-}
-
-fn regress(y: &[f64], x: &[Vec<f64>], k: usize, constant: bool) -> R<Fit> {
-    let n = y.len();
-    if n == 0 || x.len() != n {
-        return Err(CellError::Value);
-    }
-    check_regression_size(n, k)?;
-    let nf = n as f64;
-    let (xm, ym) = if constant {
-        let xm: Vec<f64> = (0..k).map(|j| x.iter().map(|r| r.get(j).copied().unwrap_or(0.0)).sum::<f64>() / nf).collect();
-        (xm, y.iter().sum::<f64>() / nf)
-    } else {
-        (vec![0.0; k], 0.0)
-    };
-    // Normal equations on (centred) data.
-    let mut xtx = vec![vec![0.0; k]; k];
-    let mut xty = vec![0.0; k];
-    for (row, yv) in x.iter().zip(y) {
-        let d: Vec<f64> = (0..k).map(|j| row.get(j).copied().unwrap_or(0.0) - xm.get(j).copied().unwrap_or(0.0)).collect();
-        let yd = yv - ym;
-        for (i, di) in d.iter().enumerate() {
-            if let Some(t) = xty.get_mut(i) {
-                *t += di * yd;
-            }
-            if let Some(r) = xtx.get_mut(i) {
-                for (j, dj) in d.iter().enumerate() {
-                    if let Some(v) = r.get_mut(j) {
-                        *v += di * dj;
-                    }
-                }
-            }
-        }
-    }
-    let (inv, dropped) = sym_inverse(&xtx);
-    let coef: Vec<f64> = (0..k).map(|i| inv.get(i).map(|r| r.iter().zip(&xty).map(|(a, b)| a * b).sum()).unwrap_or(0.0)).collect();
-    let b = if constant { ym - coef.iter().zip(&xm).map(|(c, m)| c * m).sum::<f64>() } else { 0.0 };
-    let mut ssresid = 0.0;
-    let mut sstot = 0.0;
-    for (row, yv) in x.iter().zip(y) {
-        let pred = b + coef.iter().zip(row).map(|(c, v)| c * v).sum::<f64>();
-        ssresid += (yv - pred).powi(2);
-        sstot += if constant { (yv - ym).powi(2) } else { yv * yv };
-    }
-    let used = dropped.iter().filter(|d| !**d).count() + usize::from(constant);
-    let df = nf - used as f64;
-    let ssreg = sstot - ssresid;
-    let r2 = if sstot == 0.0 { Err(CellError::Num) } else { Ok(ssreg / sstot) };
-    let sey2 = if df > 0.0 { Ok(ssresid / df) } else { Err(CellError::Num) };
-    let sey = sey2.map(f64::sqrt);
-    let v1 = if constant { nf - df - 1.0 } else { nf - df };
-    let f = match sey2 {
-        Ok(s2) if s2 > 0.0 && v1 > 0.0 => Ok((ssreg / v1) / s2),
-        _ => Err(CellError::Num),
-    };
-    let s2 = sey2.unwrap_or(f64::NAN);
-    let se: Vec<f64> = (0..k)
-        .map(|i| if dropped.get(i).copied().unwrap_or(true) { 0.0 } else { (inv.get(i).and_then(|r| r.get(i)).copied().unwrap_or(0.0) * s2).sqrt() })
-        .collect();
-    let se_b = if constant {
-        let mut q = 0.0;
-        for (i, mi) in xm.iter().enumerate() {
-            for (j, mj) in xm.iter().enumerate() {
-                q += mi * mj * inv.get(i).and_then(|r| r.get(j)).copied().unwrap_or(0.0);
-            }
-        }
-        Some((s2 * (1.0 / nf + q)).sqrt())
-    } else {
-        None
-    };
-    Ok(Fit { coef, b, se, se_b, r2, sey, f, df, ssreg, ssresid })
-}
-
-fn const_arg(a: &[Arg], i: usize) -> R<bool> {
-    if has(a, i) { boolean(a, i) } else { Ok(true) }
-}
-
-fn opt_value(a: &[Arg], i: usize) -> Option<&Value> {
-    if has(a, i) { a.get(i).map(|x| &x.value) } else { None }
-}
-
-fn linest_impl(a: &[Arg], log: bool) -> R<Value> {
-    let mut xy = parse_xy(&arg(a, 0)?.value, opt_value(a, 1))?;
-    if log {
-        if xy.y.iter().any(|&v| v <= 0.0) {
-            return Err(CellError::Num);
-        }
-        xy.y.iter_mut().for_each(|v| *v = v.ln());
-    }
-    let constant = const_arg(a, 2)?;
-    let stats = if has(a, 3) { boolean(a, 3)? } else { false };
-    let fit = regress(&xy.y, &xy.x, xy.k, constant)?;
-    let tr = |v: f64| if log { v.exp() } else { v };
-    let cols = xy.k + 1;
-    let mut data: Vec<Value> = Vec::new();
-    for c in fit.coef.iter().rev() {
-        data.push(Value::number(tr(*c)));
-    }
-    data.push(Value::number(tr(fit.b)));
-    if stats {
-        let rv = |r: R<f64>| r.map(Value::number).unwrap_or_else(Value::Error);
-        let na = Value::Error(CellError::NA);
-        for s in fit.se.iter().rev() {
-            data.push(Value::number(*s));
-        }
-        data.push(fit.se_b.map(Value::number).unwrap_or(na.clone()));
-        let rows: [[Value; 2]; 3] =
-            [[rv(fit.r2), rv(fit.sey)], [rv(fit.f), Value::number(fit.df)], [Value::number(fit.ssreg), Value::number(fit.ssresid)]];
-        for [p, q] in rows {
-            data.push(p);
-            if cols > 1 {
-                data.push(q);
-            }
-            for _ in 2..cols {
-                data.push(na.clone());
-            }
-        }
-        return array_val(5, cols, data);
-    }
-    array_val(1, cols, data)
-}
-
-fn trend_impl(a: &[Arg], log: bool) -> R<Value> {
-    let mut xy = parse_xy(&arg(a, 0)?.value, opt_value(a, 1))?;
-    if log {
-        if xy.y.iter().any(|&v| v <= 0.0) {
-            return Err(CellError::Num);
-        }
-        xy.y.iter_mut().for_each(|v| *v = v.ln());
-    }
-    let constant = const_arg(a, 3)?;
-    let fit = regress(&xy.y, &xy.x, xy.k, constant)?;
-    let predict = |row: &[f64]| -> Value {
-        let v = fit.b + fit.coef.iter().zip(row).map(|(c, x)| c * x).sum::<f64>();
-        Value::number(if log { v.exp() } else { v })
-    };
-    match opt_value(a, 2) {
-        None => {
-            let data: Vec<Value> = xy.x.iter().map(|r| predict(r)).collect();
-            array_val(xy.y_rows, xy.y_cols, data)
-        }
-        Some(nv) => {
-            let na = as_array(nv);
-            if xy.k == 1 {
-                let data = na.data.iter().map(|v| strict_num(v).map(|x| predict(&[x]))).collect::<R<Vec<_>>>()?;
-                array_val(na.rows, na.cols, data)
-            } else if !xy.y_is_row && na.cols == xy.k {
-                let mut data = Vec::with_capacity(na.rows);
-                for r in 0..na.rows {
-                    let row = (0..na.cols).map(|c| strict_num(na.get(r, c).unwrap_or(&Value::Empty))).collect::<R<Vec<_>>>()?;
-                    data.push(predict(&row));
-                }
-                array_val(na.rows, 1, data)
-            } else if xy.y_is_row && na.rows == xy.k {
-                let mut data = Vec::with_capacity(na.cols);
-                for c in 0..na.cols {
-                    let row = (0..na.rows).map(|r| strict_num(na.get(r, c).unwrap_or(&Value::Empty))).collect::<R<Vec<_>>>()?;
-                    data.push(predict(&row));
-                }
-                array_val(1, na.cols, data)
-            } else {
-                Err(CellError::Ref)
-            }
-        }
-    }
-}
-
-// ---------------------------------------------------------------------------------------------
-// Distributions
-// ---------------------------------------------------------------------------------------------
-
-fn ln_comb(n: f64, k: f64) -> f64 {
-    ln_gamma(n + 1.0) - ln_gamma(k + 1.0) - ln_gamma(n - k + 1.0)
-}
-
-fn t_upper(t: f64, v: f64) -> f64 {
-    // P(T > t) for t ≥ 0.
-    0.5 * reg_inc_beta(v / (v + t * t), v / 2.0, 0.5)
-}
-
-fn t_cdf(t: f64, v: f64) -> f64 {
-    if t >= 0.0 { 1.0 - t_upper(t, v) } else { t_upper(-t, v) }
-}
-
-fn t_pdf(t: f64, v: f64) -> f64 {
-    (ln_gamma((v + 1.0) / 2.0) - ln_gamma(v / 2.0) - 0.5 * (v * std::f64::consts::PI).ln() - (v + 1.0) / 2.0 * (1.0 + t * t / v).ln()).exp()
-}
-
-/// t ≥ 0 with P(T > t) = q, q in (0, 0.5].
-fn t_inv_upper(q: f64, v: f64) -> f64 {
-    if q >= 0.5 {
-        return 0.0;
-    }
-    invert(|t| t_upper(t, v), q, 0.0, 1.0, 1e60)
-}
-
-fn chi_inv(p: f64, v: f64) -> f64 {
-    if p <= 0.0 {
-        return 0.0;
-    }
-    invert(|x| reg_gamma_p(v / 2.0, x / 2.0), p, 0.0, v.max(1.0), 1e60)
-}
-
-fn chi_inv_rt(q: f64, v: f64) -> f64 {
-    if q >= 1.0 {
-        return 0.0;
-    }
-    invert(|x| reg_gamma_q(v / 2.0, x / 2.0), q, 0.0, v.max(1.0), 1e60)
-}
-
-fn f_cdf(x: f64, d1: f64, d2: f64) -> f64 {
-    if x <= 0.0 {
-        return 0.0;
-    }
-    reg_inc_beta(d1 * x / (d1 * x + d2), d1 / 2.0, d2 / 2.0)
-}
-
-fn f_upper(x: f64, d1: f64, d2: f64) -> f64 {
-    if x <= 0.0 {
-        return 1.0;
-    }
-    reg_inc_beta(d2 / (d2 + d1 * x), d2 / 2.0, d1 / 2.0)
-}
-
-fn f_pdf(x: f64, d1: f64, d2: f64) -> f64 {
-    if x == 0.0 {
-        return if d1 < 2.0 {
-            f64::INFINITY
-        } else if d1 == 2.0 {
-            1.0
-        } else {
-            0.0
-        };
-    }
-    ((d1 / 2.0) * (d1 / d2).ln() + (d1 / 2.0 - 1.0) * x.ln() - ((d1 + d2) / 2.0) * (1.0 + d1 * x / d2).ln() - ln_beta(d1 / 2.0, d2 / 2.0)).exp()
-}
-
-fn binom_pmf(k: f64, n: f64, p: f64) -> f64 {
-    if k < 0.0 || k > n {
-        return 0.0;
-    }
-    if p == 0.0 {
-        return if k == 0.0 { 1.0 } else { 0.0 };
-    }
-    if p == 1.0 {
-        return if k == n { 1.0 } else { 0.0 };
-    }
-    (ln_comb(n, k) + k * p.ln() + (n - k) * (1.0 - p).ln()).exp()
-}
-
-fn binom_cdf(k: f64, n: f64, p: f64) -> f64 {
-    if k < 0.0 {
-        return 0.0;
-    }
-    if k >= n {
-        return 1.0;
-    }
-    if n <= 10_000.0 {
-        let mut s = 0.0;
-        let mut i = 0.0;
-        while i <= k {
-            s += binom_pmf(i, n, p);
-            i += 1.0;
-        }
-        return s.min(1.0);
-    }
-    if p == 0.0 {
-        return 1.0;
-    }
-    if p == 1.0 {
-        return 0.0;
-    }
-    reg_inc_beta(1.0 - p, n - k, k + 1.0)
-}
-
-fn poisson_pmf(x: f64, m: f64) -> f64 {
-    if m == 0.0 {
-        return if x == 0.0 { 1.0 } else { 0.0 };
-    }
-    (-m + x * m.ln() - ln_gamma(x + 1.0)).exp()
-}
-
-fn poisson_cdf(x: f64, m: f64) -> f64 {
-    if m == 0.0 {
-        return 1.0;
-    }
-    if x <= 1000.0 {
-        let mut s = 0.0;
-        let mut i = 0.0;
-        while i <= x {
-            s += poisson_pmf(i, m);
-            i += 1.0;
-        }
-        return s.min(1.0);
-    }
-    reg_gamma_q(x + 1.0, m)
-}
-
-fn norm_dist(a: &[Arg], _: C) -> R<Value> {
-    let (x, m, sd, cum) = (num(a, 0)?, num(a, 1)?, num(a, 2)?, boolean(a, 3)?);
-    if sd <= 0.0 {
-        return Err(CellError::Num);
-    }
-    let z = (x - m) / sd;
-    num_val(if cum { norm_cdf(z) } else { norm_pdf(z) / sd })
-}
-
-fn norm_inv_fn(a: &[Arg], _: C) -> R<Value> {
-    let (p, m, sd) = (num(a, 0)?, num(a, 1)?, num(a, 2)?);
-    if !(p > 0.0 && p < 1.0) || sd <= 0.0 {
-        return Err(CellError::Num);
-    }
-    num_val(m + sd * norm_inv(p))
-}
-
-fn norm_s_dist(a: &[Arg], _: C) -> R<Value> {
-    let z = num(a, 0)?;
-    let cum = if a.len() > 1 { boolean(a, 1)? } else { true };
-    num_val(if cum { norm_cdf(z) } else { norm_pdf(z) })
-}
-
-fn norm_s_inv(a: &[Arg], _: C) -> R<Value> {
-    let p = num(a, 0)?;
-    if !(p > 0.0 && p < 1.0) {
-        return Err(CellError::Num);
-    }
-    num_val(norm_inv(p))
-}
-
-fn df_arg(a: &[Arg], i: usize) -> R<f64> {
-    let v = num(a, i)?.trunc();
-    if v < 1.0 {
-        return Err(CellError::Num);
-    }
-    Ok(v)
-}
-
-fn t_dist(a: &[Arg], _: C) -> R<Value> {
-    let x = num(a, 0)?;
-    let v = df_arg(a, 1)?;
-    let cum = boolean(a, 2)?;
-    num_val(if cum { t_cdf(x, v) } else { t_pdf(x, v) })
-}
-
-fn t_dist_2t(a: &[Arg], _: C) -> R<Value> {
-    let x = num(a, 0)?;
-    let v = df_arg(a, 1)?;
-    if x < 0.0 {
-        return Err(CellError::Num);
-    }
-    num_val(2.0 * t_upper(x, v))
-}
-
-fn t_dist_rt(a: &[Arg], _: C) -> R<Value> {
-    let x = num(a, 0)?;
-    let v = df_arg(a, 1)?;
-    num_val(1.0 - t_cdf(x, v))
-}
-
-fn tdist_legacy(a: &[Arg], _: C) -> R<Value> {
-    let x = num(a, 0)?;
-    let v = df_arg(a, 1)?;
-    let tails = num(a, 2)?.trunc();
-    if x < 0.0 || !(tails == 1.0 || tails == 2.0) {
-        return Err(CellError::Num);
-    }
-    num_val(tails * t_upper(x, v))
-}
-
-fn t_inv(a: &[Arg], _: C) -> R<Value> {
-    let p = num(a, 0)?;
-    let v = df_arg(a, 1)?;
-    if !(p > 0.0 && p < 1.0) {
-        return Err(CellError::Num);
-    }
-    num_val(if p < 0.5 { -t_inv_upper(p, v) } else { t_inv_upper(1.0 - p, v) })
-}
-
-fn t_inv_2t(a: &[Arg], _: C) -> R<Value> {
-    let p = num(a, 0)?;
-    let v = df_arg(a, 1)?;
-    if !(p > 0.0 && p <= 1.0) {
-        return Err(CellError::Num);
-    }
-    num_val(t_inv_upper(p / 2.0, v))
-}
-
-fn t_test(a: &[Arg], _: C) -> R<Value> {
-    let tails = num(a, 2)?.trunc();
-    let ty = num(a, 3)?.trunc();
-    if !(tails == 1.0 || tails == 2.0) || !(1.0..=3.0).contains(&ty) {
-        return Err(CellError::Num);
-    }
-    let (t, v) = if ty == 1.0 {
-        let (x, y) = arg_pairs(a, 0, 1)?;
-        let d: Vec<f64> = x.iter().zip(&y).map(|(p, q)| p - q).collect();
-        let n = d.len() as f64;
-        let sd = var_s(&d)?.sqrt();
-        if sd == 0.0 {
-            return Err(CellError::Div0);
-        }
-        (mean(&d)? / (sd / n.sqrt()), n - 1.0)
-    } else {
-        let x = arg_nums(a, 0)?;
-        let y = arg_nums(a, 1)?;
-        let (n1, n2) = (x.len() as f64, y.len() as f64);
-        let (v1, v2) = (var_s(&x)?, var_s(&y)?);
-        let diff = mean(&x)? - mean(&y)?;
-        if ty == 2.0 {
-            let df = n1 + n2 - 2.0;
-            let sp = ((n1 - 1.0) * v1 + (n2 - 1.0) * v2) / df;
-            let se = (sp * (1.0 / n1 + 1.0 / n2)).sqrt();
-            if se == 0.0 {
-                return Err(CellError::Div0);
-            }
-            (diff / se, df)
-        } else {
-            let (a1, a2) = (v1 / n1, v2 / n2);
-            let se = (a1 + a2).sqrt();
-            if se == 0.0 {
-                return Err(CellError::Div0);
-            }
-            let df = (a1 + a2).powi(2) / (a1 * a1 / (n1 - 1.0) + a2 * a2 / (n2 - 1.0));
-            (diff / se, df)
-        }
-    };
-    num_val(tails * t_upper(t.abs(), v))
-}
-
-fn chisq_dist(a: &[Arg], _: C) -> R<Value> {
-    let x = num(a, 0)?;
-    let v = df_arg(a, 1)?;
-    let cum = boolean(a, 2)?;
-    if x < 0.0 || v > 1e10 {
-        return Err(CellError::Num);
-    }
-    if cum {
-        return num_val(reg_gamma_p(v / 2.0, x / 2.0));
-    }
-    if x == 0.0 {
-        return if v == 2.0 {
-            num_val(0.5)
-        } else if v < 2.0 {
-            Err(CellError::Num)
-        } else {
-            num_val(0.0)
-        };
-    }
-    num_val(((v / 2.0 - 1.0) * x.ln() - x / 2.0 - (v / 2.0) * 2f64.ln() - ln_gamma(v / 2.0)).exp())
-}
-
-fn chisq_dist_rt(a: &[Arg], _: C) -> R<Value> {
-    let x = num(a, 0)?;
-    let v = df_arg(a, 1)?;
-    if x < 0.0 || v > 1e10 {
-        return Err(CellError::Num);
-    }
-    num_val(reg_gamma_q(v / 2.0, x / 2.0))
-}
-
-fn chisq_inv(a: &[Arg], _: C) -> R<Value> {
-    let p = num(a, 0)?;
-    let v = df_arg(a, 1)?;
-    if !(0.0..1.0).contains(&p) || v > 1e10 {
-        return Err(CellError::Num);
-    }
-    num_val(chi_inv(p, v))
-}
-
-fn chisq_inv_rt(a: &[Arg], _: C) -> R<Value> {
-    let p = num(a, 0)?;
-    let v = df_arg(a, 1)?;
-    if !(p > 0.0 && p <= 1.0) || v > 1e10 {
-        return Err(CellError::Num);
-    }
-    num_val(chi_inv_rt(p, v))
-}
-
-fn chisq_test(a: &[Arg], _: C) -> R<Value> {
-    let act = as_array(&arg(a, 0)?.value).into_owned();
-    let exp = as_array(&arg(a, 1)?.value).into_owned();
-    if act.rows != exp.rows || act.cols != exp.cols {
-        return Err(CellError::NA);
-    }
-    let mut x2 = 0.0;
-    for (o, e) in act.data.iter().zip(&exp.data) {
-        if let Value::Error(err) = o {
-            return Err(*err);
-        }
-        if let Value::Error(err) = e {
-            return Err(*err);
-        }
-        if let (Value::Number(o), Value::Number(e)) = (o, e) {
-            if *e == 0.0 {
-                return Err(CellError::Div0);
-            }
-            x2 += (o - e).powi(2) / e;
-        }
-    }
-    let df = if act.rows > 1 && act.cols > 1 { ((act.rows - 1) * (act.cols - 1)) as f64 } else { (act.data.len().saturating_sub(1)) as f64 };
-    if df < 1.0 {
-        return Err(CellError::NA);
-    }
-    num_val(reg_gamma_q(df / 2.0, x2 / 2.0))
-}
-
-fn f_args(a: &[Arg]) -> R<(f64, f64, f64)> {
-    let x = num(a, 0)?;
-    let d1 = num(a, 1)?.trunc();
-    let d2 = num(a, 2)?.trunc();
-    if d1 < 1.0 || d2 < 1.0 || d1 >= 1e10 || d2 >= 1e10 {
-        return Err(CellError::Num);
-    }
-    Ok((x, d1, d2))
-}
-
-fn f_dist(a: &[Arg], _: C) -> R<Value> {
-    let (x, d1, d2) = f_args(a)?;
-    let cum = boolean(a, 3)?;
-    if x < 0.0 {
-        return Err(CellError::Num);
-    }
-    num_val(if cum { f_cdf(x, d1, d2) } else { f_pdf(x, d1, d2) })
-}
-
-fn f_dist_rt(a: &[Arg], _: C) -> R<Value> {
-    let (x, d1, d2) = f_args(a)?;
-    if x < 0.0 {
-        return Err(CellError::Num);
-    }
-    num_val(f_upper(x, d1, d2))
-}
-
-fn f_inv(a: &[Arg], _: C) -> R<Value> {
-    let (p, d1, d2) = f_args(a)?;
-    if !(0.0..1.0).contains(&p) {
-        return Err(CellError::Num);
-    }
-    if p == 0.0 {
-        return num_val(0.0);
-    }
-    num_val(invert(|x| f_cdf(x, d1, d2), p, 0.0, 1.0, 1e60))
-}
-
-fn f_inv_rt(a: &[Arg], _: C) -> R<Value> {
-    let (p, d1, d2) = f_args(a)?;
-    if !(p > 0.0 && p <= 1.0) {
-        return Err(CellError::Num);
-    }
-    if p == 1.0 {
-        return num_val(0.0);
-    }
-    num_val(invert(|x| f_upper(x, d1, d2), p, 0.0, 1.0, 1e60))
-}
-
-fn f_test(a: &[Arg], _: C) -> R<Value> {
-    let x = arg_nums(a, 0)?;
-    let y = arg_nums(a, 1)?;
-    let (v1, v2) = (var_s(&x)?, var_s(&y)?);
-    if v1 == 0.0 || v2 == 0.0 {
-        return Err(CellError::Div0);
-    }
-    let f = v1 / v2;
-    let (d1, d2) = ((x.len() - 1) as f64, (y.len() - 1) as f64);
-    let c = f_cdf(f, d1, d2);
-    num_val((2.0 * c.min(1.0 - c)).min(1.0))
-}
-
-fn binom_dist(a: &[Arg], _: C) -> R<Value> {
-    let k = num(a, 0)?.trunc();
-    let n = num(a, 1)?.trunc();
-    let p = num(a, 2)?;
-    let cum = boolean(a, 3)?;
-    if k < 0.0 || k > n || !(0.0..=1.0).contains(&p) {
-        return Err(CellError::Num);
-    }
-    num_val(if cum { binom_cdf(k, n, p) } else { binom_pmf(k, n, p) })
-}
-
-fn binom_dist_range(a: &[Arg], _: C) -> R<Value> {
-    let n = num(a, 0)?.trunc();
-    let p = num(a, 1)?;
-    let s1 = num(a, 2)?.trunc();
-    let s2 = if has(a, 3) { num(a, 3)?.trunc() } else { s1 };
-    if !(0.0..=1e15).contains(&n) || !(0.0..=1.0).contains(&p) || s1 < 0.0 || s1 > n || s2 < s1 || s2 > n {
-        return Err(CellError::Num);
-    }
-    if s2 - s1 > 10_000.0 {
-        return num_val(binom_cdf(s2, n, p) - binom_cdf(s1 - 1.0, n, p));
-    }
-    let mut s = 0.0;
-    let mut k = s1;
-    while k <= s2 {
-        s += binom_pmf(k, n, p);
-        k += 1.0;
-    }
-    num_val(s)
-}
-
-fn binom_inv(a: &[Arg], _: C) -> R<Value> {
-    let n = num(a, 0)?.trunc();
-    let p = num(a, 1)?;
-    let alpha = num(a, 2)?;
-    if n < 0.0 || !(0.0..=1.0).contains(&p) || !(0.0..=1.0).contains(&alpha) {
-        return Err(CellError::Num);
-    }
-    // Smallest k with cdf(k) ≥ alpha, by bisection on the integer k.
-    let (mut lo, mut hi) = (0.0f64, n);
-    if binom_cdf(0.0, n, p) >= alpha {
-        return num_val(0.0);
-    }
-    while hi - lo > 1.0 {
-        let mid = ((lo + hi) / 2.0).floor();
-        // Above 2^53 adjacent f64 values are more than 1 apart, so the gap can stop shrinking.
-        if mid <= lo || mid >= hi {
-            break;
-        }
-        if binom_cdf(mid, n, p) >= alpha {
-            hi = mid;
-        } else {
-            lo = mid;
-        }
-    }
-    num_val(hi)
-}
-
-fn poisson_dist(a: &[Arg], _: C) -> R<Value> {
-    let x = num(a, 0)?.trunc();
-    let m = num(a, 1)?;
-    let cum = boolean(a, 2)?;
-    if x < 0.0 || m < 0.0 {
-        return Err(CellError::Num);
-    }
-    num_val(if cum { poisson_cdf(x, m) } else { poisson_pmf(x, m) })
-}
-
-fn expon_dist(a: &[Arg], _: C) -> R<Value> {
-    let (x, l, cum) = (num(a, 0)?, num(a, 1)?, boolean(a, 2)?);
-    if x < 0.0 || l <= 0.0 {
-        return Err(CellError::Num);
-    }
-    num_val(if cum { -(-l * x).exp_m1() } else { l * (-l * x).exp() })
-}
-
-fn gamma_fn(a: &[Arg], _: C) -> R<Value> {
-    let x = num(a, 0)?;
-    let g = gamma(x);
-    if !g.is_finite() {
-        return Err(CellError::Num);
-    }
-    num_val(g)
-}
-
-fn gamma_dist(a: &[Arg], _: C) -> R<Value> {
-    let (x, al, be, cum) = (num(a, 0)?, num(a, 1)?, num(a, 2)?, boolean(a, 3)?);
-    if x < 0.0 || al <= 0.0 || be <= 0.0 {
-        return Err(CellError::Num);
-    }
-    if cum {
-        return num_val(reg_gamma_p(al, x / be));
-    }
-    if x == 0.0 {
-        return if al < 1.0 {
-            Err(CellError::Num)
-        } else if al == 1.0 {
-            num_val(1.0 / be)
-        } else {
-            num_val(0.0)
-        };
-    }
-    num_val(((al - 1.0) * x.ln() - x / be - al * be.ln() - ln_gamma(al)).exp())
-}
-
-fn gamma_inv(a: &[Arg], _: C) -> R<Value> {
-    let (p, al, be) = (num(a, 0)?, num(a, 1)?, num(a, 2)?);
-    if !(0.0..1.0).contains(&p) || al <= 0.0 || be <= 0.0 {
-        return Err(CellError::Num);
-    }
-    if p == 0.0 {
-        return num_val(0.0);
-    }
-    num_val(invert(|x| reg_gamma_p(al, x / be), p, 0.0, (al * be).max(1e-3), 1e60))
-}
-
-fn gammaln(a: &[Arg], _: C) -> R<Value> {
-    let x = num(a, 0)?;
-    if x <= 0.0 {
-        return Err(CellError::Num);
-    }
-    num_val(ln_gamma(x))
-}
-
-fn beta_bounds(a: &[Arg], i: usize) -> R<(f64, f64)> {
-    let lo = if has(a, i) { num(a, i)? } else { 0.0 };
-    let hi = if has(a, i + 1) { num(a, i + 1)? } else { 1.0 };
-    if lo >= hi {
-        return Err(CellError::Num);
-    }
-    Ok((lo, hi))
-}
-
-fn beta_dist(a: &[Arg], _: C) -> R<Value> {
-    let (x, al, be, cum) = (num(a, 0)?, num(a, 1)?, num(a, 2)?, boolean(a, 3)?);
-    let (lo, hi) = beta_bounds(a, 4)?;
-    if al <= 0.0 || be <= 0.0 || x < lo || x > hi {
-        return Err(CellError::Num);
-    }
-    let z = (x - lo) / (hi - lo);
-    if cum {
-        return num_val(reg_inc_beta(z, al, be));
-    }
-    if (z == 0.0 && al < 1.0) || (z == 1.0 && be < 1.0) {
-        return Err(CellError::Num);
-    }
-    let lp = if z == 0.0 || z == 1.0 {
-        let base = if z == 0.0 { al } else { be };
-        if base == 1.0 { -ln_beta(al, be) } else { f64::NEG_INFINITY }
-    } else {
-        (al - 1.0) * z.ln() + (be - 1.0) * (1.0 - z).ln() - ln_beta(al, be)
-    };
-    num_val(lp.exp() / (hi - lo))
-}
-
-fn betadist_legacy(a: &[Arg], _: C) -> R<Value> {
-    let (x, al, be) = (num(a, 0)?, num(a, 1)?, num(a, 2)?);
-    let (lo, hi) = beta_bounds(a, 3)?;
-    if al <= 0.0 || be <= 0.0 || x < lo || x > hi {
-        return Err(CellError::Num);
-    }
-    num_val(reg_inc_beta((x - lo) / (hi - lo), al, be))
-}
-
-fn beta_inv(a: &[Arg], _: C) -> R<Value> {
-    let (p, al, be) = (num(a, 0)?, num(a, 1)?, num(a, 2)?);
-    let (lo, hi) = beta_bounds(a, 3)?;
-    if al <= 0.0 || be <= 0.0 || !(0.0..=1.0).contains(&p) {
-        return Err(CellError::Num);
-    }
-    let z = if p == 0.0 {
-        0.0
-    } else if p == 1.0 {
-        1.0
-    } else {
-        invert(|x| reg_inc_beta(x, al, be), p, 0.0, 1.0, 1.0)
-    };
-    num_val(lo + z * (hi - lo))
-}
-
-fn lognorm_dist(a: &[Arg], _: C) -> R<Value> {
-    let (x, m, sd) = (num(a, 0)?, num(a, 1)?, num(a, 2)?);
-    let cum = if a.len() > 3 { boolean(a, 3)? } else { true };
-    if x <= 0.0 || sd <= 0.0 {
-        return Err(CellError::Num);
-    }
-    let z = (x.ln() - m) / sd;
-    num_val(if cum { norm_cdf(z) } else { norm_pdf(z) / (x * sd) })
-}
-
-fn lognorm_inv(a: &[Arg], _: C) -> R<Value> {
-    let (p, m, sd) = (num(a, 0)?, num(a, 1)?, num(a, 2)?);
-    if !(p > 0.0 && p < 1.0) || sd <= 0.0 {
-        return Err(CellError::Num);
-    }
-    num_val((m + sd * norm_inv(p)).exp())
-}
-
-fn weibull_dist(a: &[Arg], _: C) -> R<Value> {
-    let (x, al, be, cum) = (num(a, 0)?, num(a, 1)?, num(a, 2)?, boolean(a, 3)?);
-    if x < 0.0 || al <= 0.0 || be <= 0.0 {
-        return Err(CellError::Num);
-    }
-    let t = (x / be).powf(al);
-    num_val(if cum { -(-t).exp_m1() } else { al / be.powf(al) * x.powf(al - 1.0) * (-t).exp() })
-}
-
-fn hypgeom_dist(a: &[Arg], _: C) -> R<Value> {
-    let k = num(a, 0)?.trunc();
-    let n = num(a, 1)?.trunc();
-    let kk = num(a, 2)?.trunc();
-    let nn = num(a, 3)?.trunc();
-    let cum = if a.len() > 4 { boolean(a, 4)? } else { false };
-    if k < 0.0 || n <= 0.0 || kk <= 0.0 || nn <= 0.0 || n > nn || kk > nn || k > n.min(kk) || k < (n - nn + kk).max(0.0) {
-        return Err(CellError::Num);
-    }
-    let pmf = |i: f64| (ln_comb(kk, i) + ln_comb(nn - kk, n - i) - ln_comb(nn, n)).exp();
-    if !cum {
-        return num_val(pmf(k));
-    }
-    let mut s = 0.0;
-    let mut i = (n - nn + kk).max(0.0);
-    if k - i > 1e6 || i + 1.0 == i {
-        return Err(CellError::Num);
-    }
-    while i <= k {
-        s += pmf(i);
-        i += 1.0;
-    }
-    num_val(s.min(1.0))
-}
-
-fn negbinom_dist(a: &[Arg], _: C) -> R<Value> {
-    let f = num(a, 0)?.trunc();
-    let s = num(a, 1)?.trunc();
-    let p = num(a, 2)?;
-    let cum = if a.len() > 3 { boolean(a, 3)? } else { false };
-    if f < 0.0 || s < 1.0 || !(0.0..=1.0).contains(&p) {
-        return Err(CellError::Num);
-    }
-    if cum {
-        if p == 0.0 {
-            return num_val(0.0);
-        }
-        return num_val(reg_inc_beta(p, s, f + 1.0));
-    }
-    if p == 0.0 || p == 1.0 {
-        return num_val(if p == 1.0 && f == 0.0 { 1.0 } else { 0.0 });
-    }
-    num_val((ln_comb(f + s - 1.0, s - 1.0) + s * p.ln() + f * (1.0 - p).ln()).exp())
-}
-
-fn confidence_norm(a: &[Arg], _: C) -> R<Value> {
-    let (al, sd, n) = (num(a, 0)?, num(a, 1)?, num(a, 2)?.trunc());
-    if !(al > 0.0 && al < 1.0) || sd <= 0.0 || n < 1.0 {
-        return Err(CellError::Num);
-    }
-    num_val(norm_inv(1.0 - al / 2.0) * sd / n.sqrt())
-}
-
-fn confidence_t(a: &[Arg], _: C) -> R<Value> {
-    let (al, sd, n) = (num(a, 0)?, num(a, 1)?, num(a, 2)?.trunc());
-    if !(al > 0.0 && al < 1.0) || sd <= 0.0 || n < 1.0 {
-        return Err(CellError::Num);
-    }
-    if n == 1.0 {
-        return Err(CellError::Div0);
-    }
-    num_val(t_inv_upper(al / 2.0, n - 1.0) * sd / n.sqrt())
-}
-
-fn z_test(a: &[Arg], _: C) -> R<Value> {
-    let v = arg_nums(a, 0)?;
-    let x = num(a, 1)?;
-    if v.is_empty() {
-        return Err(CellError::NA);
-    }
-    let sd = if has(a, 2) { num(a, 2)? } else { var_s(&v)?.sqrt() };
-    if sd <= 0.0 {
-        return Err(if has(a, 2) { CellError::Num } else { CellError::Div0 });
-    }
-    let z = (mean(&v)? - x) / (sd / (v.len() as f64).sqrt());
-    num_val(1.0 - norm_cdf(z))
-}
-
-fn fisher(a: &[Arg], _: C) -> R<Value> {
-    let x = num(a, 0)?;
-    if x <= -1.0 || x >= 1.0 {
-        return Err(CellError::Num);
-    }
-    num_val(0.5 * ((1.0 + x) / (1.0 - x)).ln())
-}
-
-fn gauss(a: &[Arg], _: C) -> R<Value> {
-    let z = num(a, 0)?;
-    // Φ(z) − 0.5 = erf(z/√2)/2, accurate near zero.
-    num_val(0.5 * erf(z / std::f64::consts::SQRT_2))
-}
-
-// ---------------------------------------------------------------------------------------------
-// SUBTOTAL / AGGREGATE support
-// ---------------------------------------------------------------------------------------------
-
-/// Evaluates SUBTOTAL/AGGREGATE function `function_num` over already-filtered cell values.
-///
-/// Codes 1–11 (and 101–111) are AVERAGE, COUNT, COUNTA, MAX, MIN, PRODUCT, STDEV.S, STDEV.P,
-/// SUM, VAR.S, VAR.P; 12 = MEDIAN, 13 = MODE.SNGL. Codes 14–19 need `k`: use
-/// [`aggregate_values_k`] (they give `#VALUE!` here). Values that are not numbers are ignored
-/// (COUNTA counts non-empty values); the first error value present is returned.
-pub fn aggregate_values(function_num: u32, values: &[Value]) -> Value {
-    crate::util::finish(aggregate_impl(function_num, values, None))
-}
-
-/// AGGREGATE codes 14 = LARGE, 15 = SMALL, 16 = PERCENTILE.INC, 17 = QUARTILE.INC,
-/// 18 = PERCENTILE.EXC, 19 = QUARTILE.EXC with their `k` argument; codes 1–13 delegate to
-/// [`aggregate_values`].
-pub fn aggregate_values_k(function_num: u32, values: &[Value], k: &Value) -> Value {
-    crate::util::finish(aggregate_impl(function_num, values, Some(k)))
-}
-
-fn aggregate_impl(function_num: u32, values: &[Value], k: Option<&Value>) -> R<Value> {
-    let code = if (101..=111).contains(&function_num) { function_num - 100 } else { function_num };
-    let mut xs = Vec::new();
-    let mut nonempty = 0usize;
-    let mut push = |v: &Value| -> R<()> {
-        match v {
-            Value::Number(n) => xs.push(*n),
-            Value::Error(e) => return Err(*e),
-            _ => {}
-        }
-        if !v.is_empty() {
-            nonempty += 1;
-        }
-        Ok(())
-    };
-    for v in values {
-        match v {
-            Value::Array(arr) => arr.iter().try_for_each(&mut push)?,
-            other => push(other)?,
-        }
-    }
-    let k = match (code, k) {
-        (14..=19, Some(k)) => k.scalar().to_number()?,
-        (14..=19, None) => return Err(CellError::Value),
-        _ => 0.0,
-    };
-    let r = match code {
-        1 => mean(&xs)?,
-        2 => xs.len() as f64,
-        3 => nonempty as f64,
-        4 => max_of(&xs),
-        5 => min_of(&xs),
-        6 => product(&xs),
-        7 => var_s(&xs)?.sqrt(),
-        8 => var_p(&xs)?.sqrt(),
-        9 => xs.iter().sum(),
-        10 => var_s(&xs)?,
-        11 => var_p(&xs)?,
-        12 => median(xs)?,
-        13 => mode_sngl(&xs)?,
-        14 => kth(xs, k, true)?,
-        15 => kth(xs, k, false)?,
-        16 => percentile_inc(xs, k)?,
-        17 => quartile_inc(xs, k)?,
-        18 => percentile_exc(xs, k)?,
-        19 => quartile_exc(xs, k)?,
-        _ => return Err(CellError::Value),
-    };
-    num_val(r)
-}
-
 // ---------------------------------------------------------------------------------------------
 // Registration
 // ---------------------------------------------------------------------------------------------
diff --git a/crates/functions/src/stat/descriptive.rs b/crates/functions/src/stat/descriptive.rs
new file mode 100644
index 000000000..af4c53607
--- /dev/null
+++ b/crates/functions/src/stat/descriptive.rs
@@ -0,0 +1,547 @@
+use super::*;
+
+// ---------------------------------------------------------------------------------------------
+// Basic helpers over number lists
+// ---------------------------------------------------------------------------------------------
+
+pub(super) fn mean(x: &[f64]) -> R<f64> {
+    if x.is_empty() {
+        return Err(CellError::Div0);
+    }
+    Ok(x.iter().sum::<f64>() / x.len() as f64)
+}
+
+pub(super) fn devsq(x: &[f64]) -> f64 {
+    if x.is_empty() {
+        return 0.0;
+    }
+    let m = x.iter().sum::<f64>() / x.len() as f64;
+    x.iter().map(|v| (v - m) * (v - m)).sum()
+}
+
+pub(super) fn var_s(x: &[f64]) -> R<f64> {
+    if x.len() < 2 {
+        return Err(CellError::Div0);
+    }
+    Ok(devsq(x) / (x.len() - 1) as f64)
+}
+
+pub(super) fn var_p(x: &[f64]) -> R<f64> {
+    if x.is_empty() {
+        return Err(CellError::Div0);
+    }
+    Ok(devsq(x) / x.len() as f64)
+}
+
+fn sorted(mut v: Vec<f64>) -> Vec<f64> {
+    v.sort_by(f64::total_cmp);
+    v
+}
+
+pub(super) fn median(v: Vec<f64>) -> R<f64> {
+    if v.is_empty() {
+        return Err(CellError::Num);
+    }
+    let s = sorted(v);
+    let n = s.len();
+    let hi = s.get(n / 2).copied().unwrap_or(0.0);
+    if n % 2 == 1 { Ok(hi) } else { Ok(s.get(n / 2 - 1).copied().unwrap_or(0.0).midpoint(hi)) }
+}
+
+/// Values occurring at least twice, ordered by first occurrence, with their (equal, maximal) count.
+fn modes(v: &[f64]) -> Vec<f64> {
+    let mut counts: Vec<(f64, usize)> = Vec::new();
+    for &x in v {
+        if let Some(e) = counts.iter_mut().find(|e| e.0 == x) {
+            e.1 += 1;
+        } else {
+            counts.push((x, 1));
+        }
+    }
+    let best = counts.iter().map(|e| e.1).max().unwrap_or(0);
+    if best < 2 {
+        return Vec::new();
+    }
+    counts.into_iter().filter(|e| e.1 == best).map(|e| e.0).collect()
+}
+
+pub(super) fn mode_sngl(v: &[f64]) -> R<f64> {
+    modes(v).first().copied().ok_or(CellError::NA)
+}
+
+/// Numbers for MODE / MODE.SNGL / MODE.MULT. Inside multi-cell ranges and arrays only numbers
+/// count, like [`numbers`]. Scalar arguments differ from SUM: in Excel 365, text and logical
+/// values are skipped whether typed or referenced (even numeric text such as `"1"`, so
+/// `MODE(1,"1",2,3)` is `#N/A`, and `MODE(1,TRUE,TRUE,2)` is `#N/A`), but a single-cell reference
+/// to a blank cell is `#VALUE!` (`MODE(B9,2,D7)` with D7 blank). Arguments are processed left to
+/// right, so the first error wins. The results cached in POI's FormulaEvalTestData come from an
+/// older Excel that answered `#VALUE!` for typed text; they are not the reference here.
+pub(super) fn mode_numbers(args: &[Arg]) -> R<Vec<f64>> {
+    let mut out = Vec::new();
+    for a in args {
+        match &a.value {
+            Value::Array(arr) => {
+                for v in arr.iter() {
+                    match v {
+                        Value::Number(n) => out.push(*n),
+                        Value::Error(e) => return Err(*e),
+                        _ => {}
+                    }
+                }
+            }
+            Value::Number(n) => out.push(*n),
+            Value::Error(e) => return Err(*e),
+            Value::Text(_) | Value::Bool(_) => {}
+            Value::Empty => return Err(CellError::Value),
+        }
+    }
+    Ok(out)
+}
+
+fn product(v: &[f64]) -> f64 {
+    if v.is_empty() { 0.0 } else { v.iter().product() }
+}
+
+pub(super) fn min_of(v: &[f64]) -> f64 {
+    if v.is_empty() { 0.0 } else { v.iter().copied().fold(f64::INFINITY, f64::min) }
+}
+
+pub(super) fn max_of(v: &[f64]) -> f64 {
+    if v.is_empty() { 0.0 } else { v.iter().copied().fold(f64::NEG_INFINITY, f64::max) }
+}
+
+/// k-th largest (k is rounded up like Excel).
+fn kth(v: Vec<f64>, k: f64, largest: bool) -> R<f64> {
+    let k = k.ceil();
+    if k.is_nan() || k < 1.0 || k > v.len() as f64 {
+        return Err(CellError::Num);
+    }
+    let s = sorted(v);
+    let k = k as usize;
+    let idx = if largest { s.len() - k } else { k - 1 };
+    s.get(idx).copied().ok_or(CellError::Num)
+}
+
+fn interp(s: &[f64], pos: f64) -> R<f64> {
+    let i = pos.floor();
+    let frac = pos - i;
+    let i = i as usize;
+    let lo = *s.get(i).ok_or(CellError::Num)?;
+    match s.get(i + 1) {
+        Some(hi) if frac > 0.0 => Ok(lo + frac * (hi - lo)),
+        _ => Ok(lo),
+    }
+}
+
+pub(super) fn percentile_inc(v: Vec<f64>, k: f64) -> R<f64> {
+    if v.is_empty() || !(0.0..=1.0).contains(&k) {
+        return Err(CellError::Num);
+    }
+    let s = sorted(v);
+    interp(&s, k * (s.len() - 1) as f64)
+}
+
+pub(super) fn percentile_exc(v: Vec<f64>, k: f64) -> R<f64> {
+    let n = v.len() as f64;
+    if v.is_empty() || !(k > 0.0 && k < 1.0) {
+        return Err(CellError::Num);
+    }
+    let pos = k * (n + 1.0);
+    if pos < 1.0 || pos > n {
+        return Err(CellError::Num);
+    }
+    let s = sorted(v);
+    interp(&s, pos - 1.0)
+}
+
+pub(super) fn quartile_inc(v: Vec<f64>, q: f64) -> R<f64> {
+    let q = q.trunc();
+    if !(0.0..=4.0).contains(&q) {
+        return Err(CellError::Num);
+    }
+    percentile_inc(v, q / 4.0)
+}
+
+pub(super) fn quartile_exc(v: Vec<f64>, q: f64) -> R<f64> {
+    let q = q.trunc();
+    if !(1.0..=3.0).contains(&q) {
+        return Err(CellError::Num);
+    }
+    percentile_exc(v, q / 4.0)
+}
+
+pub(super) fn arg_nums(a: &[Arg], i: usize) -> R<Vec<f64>> {
+    array_numbers(&arg(a, i)?.value)
+}
+
+/// Pairs of numbers at matching positions of two arrays (pairs with a non-number skipped).
+/// Different sizes give `#N/A`; errors propagate.
+fn pairs(x: &Value, y: &Value) -> R<(Vec<f64>, Vec<f64>)> {
+    let xa = as_array(x);
+    let ya = as_array(y);
+    if xa.data.len() != ya.data.len() {
+        return Err(CellError::NA);
+    }
+    let mut xs = Vec::new();
+    let mut ys = Vec::new();
+    for (a, b) in xa.data.iter().zip(ya.data.iter()) {
+        if let Value::Error(e) = a {
+            return Err(*e);
+        }
+        if let Value::Error(e) = b {
+            return Err(*e);
+        }
+        if let (Value::Number(p), Value::Number(q)) = (a, b) {
+            xs.push(*p);
+            ys.push(*q);
+        }
+    }
+    Ok((xs, ys))
+}
+
+pub(super) fn arg_pairs(a: &[Arg], i: usize, j: usize) -> R<(Vec<f64>, Vec<f64>)> {
+    pairs(&arg(a, i)?.value, &arg(a, j)?.value)
+}
+
+/// (sum of cross deviations, devsq x, devsq y, mean x, mean y).
+pub(super) fn moments(x: &[f64], y: &[f64]) -> R<(f64, f64, f64, f64, f64)> {
+    let mx = mean(x)?;
+    let my = mean(y)?;
+    let mut sxy = 0.0;
+    let mut sxx = 0.0;
+    let mut syy = 0.0;
+    for (a, b) in x.iter().zip(y) {
+        sxy += (a - mx) * (b - my);
+        sxx += (a - mx) * (a - mx);
+        syy += (b - my) * (b - my);
+    }
+    Ok((sxy, sxx, syy, mx, my))
+}
+
+pub(super) fn correl(x: &[f64], y: &[f64]) -> R<f64> {
+    if x.len() < 2 {
+        return Err(CellError::Div0);
+    }
+    let (sxy, sxx, syy, _, _) = moments(x, y)?;
+    if sxx == 0.0 || syy == 0.0 {
+        return Err(CellError::Div0);
+    }
+    Ok(sxy / (sxx * syy).sqrt())
+}
+
+pub(super) fn slope(y: &[f64], x: &[f64]) -> R<f64> {
+    let (sxy, sxx, _, _, _) = moments(x, y)?;
+    if sxx == 0.0 {
+        return Err(CellError::Div0);
+    }
+    Ok(sxy / sxx)
+}
+
+pub(super) fn intercept(y: &[f64], x: &[f64]) -> R<f64> {
+    let (sxy, sxx, _, mx, my) = moments(x, y)?;
+    if sxx == 0.0 {
+        return Err(CellError::Div0);
+    }
+    Ok(my - sxy / sxx * mx)
+}
+
+fn trunc_digits(r: f64, sig: f64) -> f64 {
+    let f = 10f64.powi(sig as i32);
+    if !f.is_finite() {
+        return r;
+    }
+    ((r * f) + 1e-9).floor() / f
+}
+
+fn col_value(v: Vec<f64>) -> R<Value> {
+    let n = v.len();
+    array_val(n, 1, v.into_iter().map(Value::Number).collect())
+}
+
+// ---------------------------------------------------------------------------------------------
+// Aggregates
+// ---------------------------------------------------------------------------------------------
+
+pub(super) fn count_impl(a: &[Arg]) -> f64 {
+    let mut c = 0usize;
+    for x in a {
+        match &x.value {
+            Value::Array(arr) => c += arr.iter().filter(|v| v.is_number()).count(),
+            Value::Number(_) => c += 1,
+            _ if x.from_ref => {}
+            Value::Bool(_) => c += 1,
+            Value::Text(t) if gridcraft_core::parse::parse_number_text(t).is_some() => c += 1,
+            _ => {}
+        }
+    }
+    c as f64
+}
+
+pub(super) fn counta_impl(a: &[Arg]) -> f64 {
+    let mut c = 0usize;
+    for x in a {
+        match &x.value {
+            Value::Array(arr) => c += arr.iter().filter(|v| !v.is_empty()).count(),
+            Value::Empty if x.from_ref => {}
+            _ => c += 1,
+        }
+    }
+    c as f64
+}
+
+pub(super) fn countblank(a: &[Arg], _: C) -> R<Value> {
+    let v = &arg(a, 0)?.value;
+    let blank = |x: &Value| matches!(x, Value::Empty) || matches!(x, Value::Text(t) if t.is_empty());
+    let c = match v {
+        Value::Array(arr) => arr.iter().filter(|x| blank(x)).count(),
+        x => usize::from(blank(x)),
+    };
+    num_val(c as f64)
+}
+
+pub(super) fn mode_mult(a: &[Arg], _: C) -> R<Value> {
+    let v = mode_numbers(a)?;
+    let m = modes(&v);
+    if m.is_empty() {
+        return Err(CellError::NA);
+    }
+    col_value(m)
+}
+
+pub(super) fn large(a: &[Arg], _: C) -> R<Value> {
+    num_val(kth(arg_nums(a, 0)?, num(a, 1)?, true)?)
+}
+
+pub(super) fn small(a: &[Arg], _: C) -> R<Value> {
+    num_val(kth(arg_nums(a, 0)?, num(a, 1)?, false)?)
+}
+
+pub(super) fn rank_impl(a: &[Arg], avg: bool) -> R<Value> {
+    let x = num(a, 0)?;
+    let v = arg_nums(a, 1)?;
+    let asc = opt_num(a, 2, 0.0)? != 0.0;
+    let eq = v.iter().filter(|&&y| y == x).count();
+    if eq == 0 {
+        return Err(CellError::NA);
+    }
+    let before = v.iter().filter(|&&y| if asc { y < x } else { y > x }).count();
+    let r = before as f64 + 1.0;
+    num_val(if avg { r + (eq as f64 - 1.0) / 2.0 } else { r })
+}
+
+pub(super) fn percentrank_impl(a: &[Arg], exc: bool) -> R<Value> {
+    let s = sorted(arg_nums(a, 0)?);
+    let x = num(a, 1)?;
+    let sig = opt_num(a, 2, 3.0)?.trunc();
+    if sig < 1.0 || s.is_empty() {
+        return Err(CellError::Num);
+    }
+    let (first, last) = (s.first().copied().unwrap_or(0.0), s.last().copied().unwrap_or(0.0));
+    if x < first || x > last {
+        return Err(CellError::NA);
+    }
+    let n = s.len() as f64;
+    let rank_of = |v: f64| -> f64 {
+        let less = s.iter().filter(|&&y| y < v).count() as f64;
+        if exc {
+            (less + 1.0) / (n + 1.0)
+        } else if n <= 1.0 {
+            1.0
+        } else {
+            less / (n - 1.0)
+        }
+    };
+    let r = if s.contains(&x) {
+        rank_of(x)
+    } else {
+        let lo = s.iter().copied().filter(|&y| y < x).fold(f64::NEG_INFINITY, f64::max);
+        let hi = s.iter().copied().filter(|&y| y > x).fold(f64::INFINITY, f64::min);
+        let (rl, rh) = (rank_of(lo), rank_of(hi));
+        rl + (x - lo) / (hi - lo) * (rh - rl)
+    };
+    num_val(trunc_digits(r, sig))
+}
+
+pub(super) fn trimmean(a: &[Arg], _: C) -> R<Value> {
+    let v = arg_nums(a, 0)?;
+    let p = num(a, 1)?;
+    if !(0.0..1.0).contains(&p) || v.is_empty() {
+        return Err(CellError::Num);
+    }
+    let n = v.len();
+    let per_side = ((n as f64 * p) / 2.0).floor() as usize;
+    let s = sorted(v);
+    let kept = s.get(per_side..n - per_side).unwrap_or(&[]);
+    num_val(mean(kept)?)
+}
+
+pub(super) fn kurt(v: &[f64]) -> R<f64> {
+    let n = v.len() as f64;
+    if v.len() < 4 {
+        return Err(CellError::Div0);
+    }
+    let m = mean(v)?;
+    let s = var_s(v)?.sqrt();
+    if s == 0.0 {
+        return Err(CellError::Div0);
+    }
+    let sum4: f64 = v.iter().map(|x| ((x - m) / s).powi(4)).sum();
+    Ok(n * (n + 1.0) / ((n - 1.0) * (n - 2.0) * (n - 3.0)) * sum4 - 3.0 * (n - 1.0).powi(2) / ((n - 2.0) * (n - 3.0)))
+}
+
+pub(super) fn skew(v: &[f64]) -> R<f64> {
+    let n = v.len() as f64;
+    if v.len() < 3 {
+        return Err(CellError::Div0);
+    }
+    let m = mean(v)?;
+    let s = var_s(v)?.sqrt();
+    if s == 0.0 {
+        return Err(CellError::Div0);
+    }
+    let sum3: f64 = v.iter().map(|x| ((x - m) / s).powi(3)).sum();
+    Ok(n / ((n - 1.0) * (n - 2.0)) * sum3)
+}
+
+pub(super) fn skew_p(v: &[f64]) -> R<f64> {
+    let m = mean(v)?;
+    let s = var_p(v)?.sqrt();
+    if s == 0.0 {
+        return Err(CellError::Div0);
+    }
+    Ok(v.iter().map(|x| ((x - m) / s).powi(3)).sum::<f64>() / v.len() as f64)
+}
+
+pub(super) fn geomean(v: &[f64]) -> R<f64> {
+    if v.is_empty() || v.iter().any(|&x| x <= 0.0) {
+        return Err(CellError::Num);
+    }
+    Ok((v.iter().map(|x| x.ln()).sum::<f64>() / v.len() as f64).exp())
+}
+
+pub(super) fn harmean(v: &[f64]) -> R<f64> {
+    if v.is_empty() || v.iter().any(|&x| x <= 0.0) {
+        return Err(CellError::Num);
+    }
+    Ok(v.len() as f64 / v.iter().map(|x| 1.0 / x).sum::<f64>())
+}
+
+pub(super) fn avedev(v: &[f64]) -> R<f64> {
+    if v.is_empty() {
+        return Err(CellError::Num);
+    }
+    let m = mean(v)?;
+    Ok(v.iter().map(|x| (x - m).abs()).sum::<f64>() / v.len() as f64)
+}
+
+pub(super) fn frequency(a: &[Arg], _: C) -> R<Value> {
+    let data = arg_nums(a, 0)?;
+    let bins = arg_nums(a, 1)?;
+    if bins.is_empty() {
+        return col_value(vec![data.len() as f64]);
+    }
+    let mut order: Vec<usize> = (0..bins.len()).collect();
+    order.sort_by(|&i, &j| bins.get(i).copied().unwrap_or(0.0).total_cmp(&bins.get(j).copied().unwrap_or(0.0)));
+    let mut counts = vec![0.0; bins.len() + 1];
+    for x in data {
+        let slot = order.iter().find(|&&i| x <= bins.get(i).copied().unwrap_or(f64::INFINITY)).copied().unwrap_or(bins.len());
+        if let Some(c) = counts.get_mut(slot) {
+            *c += 1.0;
+        }
+    }
+    col_value(counts)
+}
+
+pub(super) fn prob(a: &[Arg], _: C) -> R<Value> {
+    let (xs, ps) = arg_pairs(a, 0, 1)?;
+    let lo = num(a, 2)?;
+    let hi = if has(a, 3) { num(a, 3)? } else { lo };
+    if ps.iter().any(|&p| p <= 0.0 || p > 1.0) {
+        return Err(CellError::Num);
+    }
+    let total: f64 = ps.iter().sum();
+    if (total - 1.0).abs() > 1e-7 {
+        return Err(CellError::Num);
+    }
+    num_val(xs.iter().zip(&ps).filter(|(x, _)| **x >= lo && **x <= hi).map(|(_, p)| p).sum())
+}
+
+pub(super) fn standardize(a: &[Arg], _: C) -> R<Value> {
+    let (x, m, sd) = (num(a, 0)?, num(a, 1)?, num(a, 2)?);
+    if sd <= 0.0 {
+        return Err(CellError::Num);
+    }
+    num_val((x - m) / sd)
+}
+
+// ---------------------------------------------------------------------------------------------
+// SUBTOTAL / AGGREGATE support
+// ---------------------------------------------------------------------------------------------
+
+/// Evaluates SUBTOTAL/AGGREGATE function `function_num` over already-filtered cell values.
+///
+/// Codes 1–11 (and 101–111) are AVERAGE, COUNT, COUNTA, MAX, MIN, PRODUCT, STDEV.S, STDEV.P,
+/// SUM, VAR.S, VAR.P; 12 = MEDIAN, 13 = MODE.SNGL. Codes 14–19 need `k`: use
+/// [`aggregate_values_k`] (they give `#VALUE!` here). Values that are not numbers are ignored
+/// (COUNTA counts non-empty values); the first error value present is returned.
+pub fn aggregate_values(function_num: u32, values: &[Value]) -> Value {
+    crate::util::finish(aggregate_impl(function_num, values, None))
+}
+
+/// AGGREGATE codes 14 = LARGE, 15 = SMALL, 16 = PERCENTILE.INC, 17 = QUARTILE.INC,
+/// 18 = PERCENTILE.EXC, 19 = QUARTILE.EXC with their `k` argument; codes 1–13 delegate to
+/// [`aggregate_values`].
+pub fn aggregate_values_k(function_num: u32, values: &[Value], k: &Value) -> Value {
+    crate::util::finish(aggregate_impl(function_num, values, Some(k)))
+}
+
+fn aggregate_impl(function_num: u32, values: &[Value], k: Option<&Value>) -> R<Value> {
+    let code = if (101..=111).contains(&function_num) { function_num - 100 } else { function_num };
+    let mut xs = Vec::new();
+    let mut nonempty = 0usize;
+    let mut push = |v: &Value| -> R<()> {
+        match v {
+            Value::Number(n) => xs.push(*n),
+            Value::Error(e) => return Err(*e),
+            _ => {}
+        }
+        if !v.is_empty() {
+            nonempty += 1;
+        }
+        Ok(())
+    };
+    for v in values {
+        match v {
+            Value::Array(arr) => arr.iter().try_for_each(&mut push)?,
+            other => push(other)?,
+        }
+    }
+    let k = match (code, k) {
+        (14..=19, Some(k)) => k.scalar().to_number()?,
+        (14..=19, None) => return Err(CellError::Value),
+        _ => 0.0,
+    };
+    let r = match code {
+        1 => mean(&xs)?,
+        2 => xs.len() as f64,
+        3 => nonempty as f64,
+        4 => max_of(&xs),
+        5 => min_of(&xs),
+        6 => product(&xs),
+        7 => var_s(&xs)?.sqrt(),
+        8 => var_p(&xs)?.sqrt(),
+        9 => xs.iter().sum(),
+        10 => var_s(&xs)?,
+        11 => var_p(&xs)?,
+        12 => median(xs)?,
+        13 => mode_sngl(&xs)?,
+        14 => kth(xs, k, true)?,
+        15 => kth(xs, k, false)?,
+        16 => percentile_inc(xs, k)?,
+        17 => quartile_inc(xs, k)?,
+        18 => percentile_exc(xs, k)?,
+        19 => quartile_exc(xs, k)?,
+        _ => return Err(CellError::Value),
+    };
+    num_val(r)
+}
diff --git a/crates/functions/src/stat/distributions.rs b/crates/functions/src/stat/distributions.rs
new file mode 100644
index 000000000..cb2f0398a
--- /dev/null
+++ b/crates/functions/src/stat/distributions.rs
@@ -0,0 +1,702 @@
+use super::descriptive::{arg_nums, arg_pairs, mean, var_s};
+use super::*;
+
+// ---------------------------------------------------------------------------------------------
+// Distributions
+// ---------------------------------------------------------------------------------------------
+
+fn ln_comb(n: f64, k: f64) -> f64 {
+    ln_gamma(n + 1.0) - ln_gamma(k + 1.0) - ln_gamma(n - k + 1.0)
+}
+
+fn t_upper(t: f64, v: f64) -> f64 {
+    // P(T > t) for t ≥ 0.
+    0.5 * reg_inc_beta(v / (v + t * t), v / 2.0, 0.5)
+}
+
+fn t_cdf(t: f64, v: f64) -> f64 {
+    if t >= 0.0 { 1.0 - t_upper(t, v) } else { t_upper(-t, v) }
+}
+
+fn t_pdf(t: f64, v: f64) -> f64 {
+    (ln_gamma((v + 1.0) / 2.0) - ln_gamma(v / 2.0) - 0.5 * (v * std::f64::consts::PI).ln() - (v + 1.0) / 2.0 * (1.0 + t * t / v).ln()).exp()
+}
+
+/// t ≥ 0 with P(T > t) = q, q in (0, 0.5].
+fn t_inv_upper(q: f64, v: f64) -> f64 {
+    if q >= 0.5 {
+        return 0.0;
+    }
+    invert(|t| t_upper(t, v), q, 0.0, 1.0, 1e60)
+}
+
+fn chi_inv(p: f64, v: f64) -> f64 {
+    if p <= 0.0 {
+        return 0.0;
+    }
+    invert(|x| reg_gamma_p(v / 2.0, x / 2.0), p, 0.0, v.max(1.0), 1e60)
+}
+
+fn chi_inv_rt(q: f64, v: f64) -> f64 {
+    if q >= 1.0 {
+        return 0.0;
+    }
+    invert(|x| reg_gamma_q(v / 2.0, x / 2.0), q, 0.0, v.max(1.0), 1e60)
+}
+
+fn f_cdf(x: f64, d1: f64, d2: f64) -> f64 {
+    if x <= 0.0 {
+        return 0.0;
+    }
+    reg_inc_beta(d1 * x / (d1 * x + d2), d1 / 2.0, d2 / 2.0)
+}
+
+fn f_upper(x: f64, d1: f64, d2: f64) -> f64 {
+    if x <= 0.0 {
+        return 1.0;
+    }
+    reg_inc_beta(d2 / (d2 + d1 * x), d2 / 2.0, d1 / 2.0)
+}
+
+fn f_pdf(x: f64, d1: f64, d2: f64) -> f64 {
+    if x == 0.0 {
+        return if d1 < 2.0 {
+            f64::INFINITY
+        } else if d1 == 2.0 {
+            1.0
+        } else {
+            0.0
+        };
+    }
+    ((d1 / 2.0) * (d1 / d2).ln() + (d1 / 2.0 - 1.0) * x.ln() - ((d1 + d2) / 2.0) * (1.0 + d1 * x / d2).ln() - ln_beta(d1 / 2.0, d2 / 2.0)).exp()
+}
+
+fn binom_pmf(k: f64, n: f64, p: f64) -> f64 {
+    if k < 0.0 || k > n {
+        return 0.0;
+    }
+    if p == 0.0 {
+        return if k == 0.0 { 1.0 } else { 0.0 };
+    }
+    if p == 1.0 {
+        return if k == n { 1.0 } else { 0.0 };
+    }
+    (ln_comb(n, k) + k * p.ln() + (n - k) * (1.0 - p).ln()).exp()
+}
+
+fn binom_cdf(k: f64, n: f64, p: f64) -> f64 {
+    if k < 0.0 {
+        return 0.0;
+    }
+    if k >= n {
+        return 1.0;
+    }
+    if n <= 10_000.0 {
+        let mut s = 0.0;
+        let mut i = 0.0;
+        while i <= k {
+            s += binom_pmf(i, n, p);
+            i += 1.0;
+        }
+        return s.min(1.0);
+    }
+    if p == 0.0 {
+        return 1.0;
+    }
+    if p == 1.0 {
+        return 0.0;
+    }
+    reg_inc_beta(1.0 - p, n - k, k + 1.0)
+}
+
+fn poisson_pmf(x: f64, m: f64) -> f64 {
+    if m == 0.0 {
+        return if x == 0.0 { 1.0 } else { 0.0 };
+    }
+    (-m + x * m.ln() - ln_gamma(x + 1.0)).exp()
+}
+
+fn poisson_cdf(x: f64, m: f64) -> f64 {
+    if m == 0.0 {
+        return 1.0;
+    }
+    if x <= 1000.0 {
+        let mut s = 0.0;
+        let mut i = 0.0;
+        while i <= x {
+            s += poisson_pmf(i, m);
+            i += 1.0;
+        }
+        return s.min(1.0);
+    }
+    reg_gamma_q(x + 1.0, m)
+}
+
+pub(super) fn norm_dist(a: &[Arg], _: C) -> R<Value> {
+    let (x, m, sd, cum) = (num(a, 0)?, num(a, 1)?, num(a, 2)?, boolean(a, 3)?);
+    if sd <= 0.0 {
+        return Err(CellError::Num);
+    }
+    let z = (x - m) / sd;
+    num_val(if cum { norm_cdf(z) } else { norm_pdf(z) / sd })
+}
+
+pub(super) fn norm_inv_fn(a: &[Arg], _: C) -> R<Value> {
+    let (p, m, sd) = (num(a, 0)?, num(a, 1)?, num(a, 2)?);
+    if !(p > 0.0 && p < 1.0) || sd <= 0.0 {
+        return Err(CellError::Num);
+    }
+    num_val(m + sd * norm_inv(p))
+}
+
+pub(super) fn norm_s_dist(a: &[Arg], _: C) -> R<Value> {
+    let z = num(a, 0)?;
+    let cum = if a.len() > 1 { boolean(a, 1)? } else { true };
+    num_val(if cum { norm_cdf(z) } else { norm_pdf(z) })
+}
+
+pub(super) fn norm_s_inv(a: &[Arg], _: C) -> R<Value> {
+    let p = num(a, 0)?;
+    if !(p > 0.0 && p < 1.0) {
+        return Err(CellError::Num);
+    }
+    num_val(norm_inv(p))
+}
+
+fn df_arg(a: &[Arg], i: usize) -> R<f64> {
+    let v = num(a, i)?.trunc();
+    if v < 1.0 {
+        return Err(CellError::Num);
+    }
+    Ok(v)
+}
+
+pub(super) fn t_dist(a: &[Arg], _: C) -> R<Value> {
+    let x = num(a, 0)?;
+    let v = df_arg(a, 1)?;
+    let cum = boolean(a, 2)?;
+    num_val(if cum { t_cdf(x, v) } else { t_pdf(x, v) })
+}
+
+pub(super) fn t_dist_2t(a: &[Arg], _: C) -> R<Value> {
+    let x = num(a, 0)?;
+    let v = df_arg(a, 1)?;
+    if x < 0.0 {
+        return Err(CellError::Num);
+    }
+    num_val(2.0 * t_upper(x, v))
+}
+
+pub(super) fn t_dist_rt(a: &[Arg], _: C) -> R<Value> {
+    let x = num(a, 0)?;
+    let v = df_arg(a, 1)?;
+    num_val(1.0 - t_cdf(x, v))
+}
+
+pub(super) fn tdist_legacy(a: &[Arg], _: C) -> R<Value> {
+    let x = num(a, 0)?;
+    let v = df_arg(a, 1)?;
+    let tails = num(a, 2)?.trunc();
+    if x < 0.0 || !(tails == 1.0 || tails == 2.0) {
+        return Err(CellError::Num);
+    }
+    num_val(tails * t_upper(x, v))
+}
+
+pub(super) fn t_inv(a: &[Arg], _: C) -> R<Value> {
+    let p = num(a, 0)?;
+    let v = df_arg(a, 1)?;
+    if !(p > 0.0 && p < 1.0) {
+        return Err(CellError::Num);
+    }
+    num_val(if p < 0.5 { -t_inv_upper(p, v) } else { t_inv_upper(1.0 - p, v) })
+}
+
+pub(super) fn t_inv_2t(a: &[Arg], _: C) -> R<Value> {
+    let p = num(a, 0)?;
+    let v = df_arg(a, 1)?;
+    if !(p > 0.0 && p <= 1.0) {
+        return Err(CellError::Num);
+    }
+    num_val(t_inv_upper(p / 2.0, v))
+}
+
+pub(super) fn t_test(a: &[Arg], _: C) -> R<Value> {
+    let tails = num(a, 2)?.trunc();
+    let ty = num(a, 3)?.trunc();
+    if !(tails == 1.0 || tails == 2.0) || !(1.0..=3.0).contains(&ty) {
+        return Err(CellError::Num);
+    }
+    let (t, v) = if ty == 1.0 {
+        let (x, y) = arg_pairs(a, 0, 1)?;
+        let d: Vec<f64> = x.iter().zip(&y).map(|(p, q)| p - q).collect();
+        let n = d.len() as f64;
+        let sd = var_s(&d)?.sqrt();
+        if sd == 0.0 {
+            return Err(CellError::Div0);
+        }
+        (mean(&d)? / (sd / n.sqrt()), n - 1.0)
+    } else {
+        let x = arg_nums(a, 0)?;
+        let y = arg_nums(a, 1)?;
+        let (n1, n2) = (x.len() as f64, y.len() as f64);
+        let (v1, v2) = (var_s(&x)?, var_s(&y)?);
+        let diff = mean(&x)? - mean(&y)?;
+        if ty == 2.0 {
+            let df = n1 + n2 - 2.0;
+            let sp = ((n1 - 1.0) * v1 + (n2 - 1.0) * v2) / df;
+            let se = (sp * (1.0 / n1 + 1.0 / n2)).sqrt();
+            if se == 0.0 {
+                return Err(CellError::Div0);
+            }
+            (diff / se, df)
+        } else {
+            let (a1, a2) = (v1 / n1, v2 / n2);
+            let se = (a1 + a2).sqrt();
+            if se == 0.0 {
+                return Err(CellError::Div0);
+            }
+            let df = (a1 + a2).powi(2) / (a1 * a1 / (n1 - 1.0) + a2 * a2 / (n2 - 1.0));
+            (diff / se, df)
+        }
+    };
+    num_val(tails * t_upper(t.abs(), v))
+}
+
+pub(super) fn chisq_dist(a: &[Arg], _: C) -> R<Value> {
+    let x = num(a, 0)?;
+    let v = df_arg(a, 1)?;
+    let cum = boolean(a, 2)?;
+    if x < 0.0 || v > 1e10 {
+        return Err(CellError::Num);
+    }
+    if cum {
+        return num_val(reg_gamma_p(v / 2.0, x / 2.0));
+    }
+    if x == 0.0 {
+        return if v == 2.0 {
+            num_val(0.5)
+        } else if v < 2.0 {
+            Err(CellError::Num)
+        } else {
+            num_val(0.0)
+        };
+    }
+    num_val(((v / 2.0 - 1.0) * x.ln() - x / 2.0 - (v / 2.0) * 2f64.ln() - ln_gamma(v / 2.0)).exp())
+}
+
+pub(super) fn chisq_dist_rt(a: &[Arg], _: C) -> R<Value> {
+    let x = num(a, 0)?;
+    let v = df_arg(a, 1)?;
+    if x < 0.0 || v > 1e10 {
+        return Err(CellError::Num);
+    }
+    num_val(reg_gamma_q(v / 2.0, x / 2.0))
+}
+
+pub(super) fn chisq_inv(a: &[Arg], _: C) -> R<Value> {
+    let p = num(a, 0)?;
+    let v = df_arg(a, 1)?;
+    if !(0.0..1.0).contains(&p) || v > 1e10 {
+        return Err(CellError::Num);
+    }
+    num_val(chi_inv(p, v))
+}
+
+pub(super) fn chisq_inv_rt(a: &[Arg], _: C) -> R<Value> {
+    let p = num(a, 0)?;
+    let v = df_arg(a, 1)?;
+    if !(p > 0.0 && p <= 1.0) || v > 1e10 {
+        return Err(CellError::Num);
+    }
+    num_val(chi_inv_rt(p, v))
+}
+
+pub(super) fn chisq_test(a: &[Arg], _: C) -> R<Value> {
+    let act = as_array(&arg(a, 0)?.value).into_owned();
+    let exp = as_array(&arg(a, 1)?.value).into_owned();
+    if act.rows != exp.rows || act.cols != exp.cols {
+        return Err(CellError::NA);
+    }
+    let mut x2 = 0.0;
+    for (o, e) in act.data.iter().zip(&exp.data) {
+        if let Value::Error(err) = o {
+            return Err(*err);
+        }
+        if let Value::Error(err) = e {
+            return Err(*err);
+        }
+        if let (Value::Number(o), Value::Number(e)) = (o, e) {
+            if *e == 0.0 {
+                return Err(CellError::Div0);
+            }
+            x2 += (o - e).powi(2) / e;
+        }
+    }
+    let df = if act.rows > 1 && act.cols > 1 { ((act.rows - 1) * (act.cols - 1)) as f64 } else { (act.data.len().saturating_sub(1)) as f64 };
+    if df < 1.0 {
+        return Err(CellError::NA);
+    }
+    num_val(reg_gamma_q(df / 2.0, x2 / 2.0))
+}
+
+fn f_args(a: &[Arg]) -> R<(f64, f64, f64)> {
+    let x = num(a, 0)?;
+    let d1 = num(a, 1)?.trunc();
+    let d2 = num(a, 2)?.trunc();
+    if d1 < 1.0 || d2 < 1.0 || d1 >= 1e10 || d2 >= 1e10 {
+        return Err(CellError::Num);
+    }
+    Ok((x, d1, d2))
+}
+
+pub(super) fn f_dist(a: &[Arg], _: C) -> R<Value> {
+    let (x, d1, d2) = f_args(a)?;
+    let cum = boolean(a, 3)?;
+    if x < 0.0 {
+        return Err(CellError::Num);
+    }
+    num_val(if cum { f_cdf(x, d1, d2) } else { f_pdf(x, d1, d2) })
+}
+
+pub(super) fn f_dist_rt(a: &[Arg], _: C) -> R<Value> {
+    let (x, d1, d2) = f_args(a)?;
+    if x < 0.0 {
+        return Err(CellError::Num);
+    }
+    num_val(f_upper(x, d1, d2))
+}
+
+pub(super) fn f_inv(a: &[Arg], _: C) -> R<Value> {
+    let (p, d1, d2) = f_args(a)?;
+    if !(0.0..1.0).contains(&p) {
+        return Err(CellError::Num);
+    }
+    if p == 0.0 {
+        return num_val(0.0);
+    }
+    num_val(invert(|x| f_cdf(x, d1, d2), p, 0.0, 1.0, 1e60))
+}
+
+pub(super) fn f_inv_rt(a: &[Arg], _: C) -> R<Value> {
+    let (p, d1, d2) = f_args(a)?;
+    if !(p > 0.0 && p <= 1.0) {
+        return Err(CellError::Num);
+    }
+    if p == 1.0 {
+        return num_val(0.0);
+    }
+    num_val(invert(|x| f_upper(x, d1, d2), p, 0.0, 1.0, 1e60))
+}
+
+pub(super) fn f_test(a: &[Arg], _: C) -> R<Value> {
+    let x = arg_nums(a, 0)?;
+    let y = arg_nums(a, 1)?;
+    let (v1, v2) = (var_s(&x)?, var_s(&y)?);
+    if v1 == 0.0 || v2 == 0.0 {
+        return Err(CellError::Div0);
+    }
+    let f = v1 / v2;
+    let (d1, d2) = ((x.len() - 1) as f64, (y.len() - 1) as f64);
+    let c = f_cdf(f, d1, d2);
+    num_val((2.0 * c.min(1.0 - c)).min(1.0))
+}
+
+pub(super) fn binom_dist(a: &[Arg], _: C) -> R<Value> {
+    let k = num(a, 0)?.trunc();
+    let n = num(a, 1)?.trunc();
+    let p = num(a, 2)?;
+    let cum = boolean(a, 3)?;
+    if k < 0.0 || k > n || !(0.0..=1.0).contains(&p) {
+        return Err(CellError::Num);
+    }
+    num_val(if cum { binom_cdf(k, n, p) } else { binom_pmf(k, n, p) })
+}
+
+pub(super) fn binom_dist_range(a: &[Arg], _: C) -> R<Value> {
+    let n = num(a, 0)?.trunc();
+    let p = num(a, 1)?;
+    let s1 = num(a, 2)?.trunc();
+    let s2 = if has(a, 3) { num(a, 3)?.trunc() } else { s1 };
+    if !(0.0..=1e15).contains(&n) || !(0.0..=1.0).contains(&p) || s1 < 0.0 || s1 > n || s2 < s1 || s2 > n {
+        return Err(CellError::Num);
+    }
+    if s2 - s1 > 10_000.0 {
+        return num_val(binom_cdf(s2, n, p) - binom_cdf(s1 - 1.0, n, p));
+    }
+    let mut s = 0.0;
+    let mut k = s1;
+    while k <= s2 {
+        s += binom_pmf(k, n, p);
+        k += 1.0;
+    }
+    num_val(s)
+}
+
+pub(super) fn binom_inv(a: &[Arg], _: C) -> R<Value> {
+    let n = num(a, 0)?.trunc();
+    let p = num(a, 1)?;
+    let alpha = num(a, 2)?;
+    if n < 0.0 || !(0.0..=1.0).contains(&p) || !(0.0..=1.0).contains(&alpha) {
+        return Err(CellError::Num);
+    }
+    // Smallest k with cdf(k) ≥ alpha, by bisection on the integer k.
+    let (mut lo, mut hi) = (0.0f64, n);
+    if binom_cdf(0.0, n, p) >= alpha {
+        return num_val(0.0);
+    }
+    while hi - lo > 1.0 {
+        let mid = ((lo + hi) / 2.0).floor();
+        // Above 2^53 adjacent f64 values are more than 1 apart, so the gap can stop shrinking.
+        if mid <= lo || mid >= hi {
+            break;
+        }
+        if binom_cdf(mid, n, p) >= alpha {
+            hi = mid;
+        } else {
+            lo = mid;
+        }
+    }
+    num_val(hi)
+}
+
+pub(super) fn poisson_dist(a: &[Arg], _: C) -> R<Value> {
+    let x = num(a, 0)?.trunc();
+    let m = num(a, 1)?;
+    let cum = boolean(a, 2)?;
+    if x < 0.0 || m < 0.0 {
+        return Err(CellError::Num);
+    }
+    num_val(if cum { poisson_cdf(x, m) } else { poisson_pmf(x, m) })
+}
+
+pub(super) fn expon_dist(a: &[Arg], _: C) -> R<Value> {
+    let (x, l, cum) = (num(a, 0)?, num(a, 1)?, boolean(a, 2)?);
+    if x < 0.0 || l <= 0.0 {
+        return Err(CellError::Num);
+    }
+    num_val(if cum { -(-l * x).exp_m1() } else { l * (-l * x).exp() })
+}
+
+pub(super) fn gamma_fn(a: &[Arg], _: C) -> R<Value> {
+    let x = num(a, 0)?;
+    let g = gamma(x);
+    if !g.is_finite() {
+        return Err(CellError::Num);
+    }
+    num_val(g)
+}
+
+pub(super) fn gamma_dist(a: &[Arg], _: C) -> R<Value> {
+    let (x, al, be, cum) = (num(a, 0)?, num(a, 1)?, num(a, 2)?, boolean(a, 3)?);
+    if x < 0.0 || al <= 0.0 || be <= 0.0 {
+        return Err(CellError::Num);
+    }
+    if cum {
+        return num_val(reg_gamma_p(al, x / be));
+    }
+    if x == 0.0 {
+        return if al < 1.0 {
+            Err(CellError::Num)
+        } else if al == 1.0 {
+            num_val(1.0 / be)
+        } else {
+            num_val(0.0)
+        };
+    }
+    num_val(((al - 1.0) * x.ln() - x / be - al * be.ln() - ln_gamma(al)).exp())
+}
+
+pub(super) fn gamma_inv(a: &[Arg], _: C) -> R<Value> {
+    let (p, al, be) = (num(a, 0)?, num(a, 1)?, num(a, 2)?);
+    if !(0.0..1.0).contains(&p) || al <= 0.0 || be <= 0.0 {
+        return Err(CellError::Num);
+    }
+    if p == 0.0 {
+        return num_val(0.0);
+    }
+    num_val(invert(|x| reg_gamma_p(al, x / be), p, 0.0, (al * be).max(1e-3), 1e60))
+}
+
+pub(super) fn gammaln(a: &[Arg], _: C) -> R<Value> {
+    let x = num(a, 0)?;
+    if x <= 0.0 {
+        return Err(CellError::Num);
+    }
+    num_val(ln_gamma(x))
+}
+
+fn beta_bounds(a: &[Arg], i: usize) -> R<(f64, f64)> {
+    let lo = if has(a, i) { num(a, i)? } else { 0.0 };
+    let hi = if has(a, i + 1) { num(a, i + 1)? } else { 1.0 };
+    if lo >= hi {
+        return Err(CellError::Num);
+    }
+    Ok((lo, hi))
+}
+
+pub(super) fn beta_dist(a: &[Arg], _: C) -> R<Value> {
+    let (x, al, be, cum) = (num(a, 0)?, num(a, 1)?, num(a, 2)?, boolean(a, 3)?);
+    let (lo, hi) = beta_bounds(a, 4)?;
+    if al <= 0.0 || be <= 0.0 || x < lo || x > hi {
+        return Err(CellError::Num);
+    }
+    let z = (x - lo) / (hi - lo);
+    if cum {
+        return num_val(reg_inc_beta(z, al, be));
+    }
+    if (z == 0.0 && al < 1.0) || (z == 1.0 && be < 1.0) {
+        return Err(CellError::Num);
+    }
+    let lp = if z == 0.0 || z == 1.0 {
+        let base = if z == 0.0 { al } else { be };
+        if base == 1.0 { -ln_beta(al, be) } else { f64::NEG_INFINITY }
+    } else {
+        (al - 1.0) * z.ln() + (be - 1.0) * (1.0 - z).ln() - ln_beta(al, be)
+    };
+    num_val(lp.exp() / (hi - lo))
+}
+
+pub(super) fn betadist_legacy(a: &[Arg], _: C) -> R<Value> {
+    let (x, al, be) = (num(a, 0)?, num(a, 1)?, num(a, 2)?);
+    let (lo, hi) = beta_bounds(a, 3)?;
+    if al <= 0.0 || be <= 0.0 || x < lo || x > hi {
+        return Err(CellError::Num);
+    }
+    num_val(reg_inc_beta((x - lo) / (hi - lo), al, be))
+}
+
+pub(super) fn beta_inv(a: &[Arg], _: C) -> R<Value> {
+    let (p, al, be) = (num(a, 0)?, num(a, 1)?, num(a, 2)?);
+    let (lo, hi) = beta_bounds(a, 3)?;
+    if al <= 0.0 || be <= 0.0 || !(0.0..=1.0).contains(&p) {
+        return Err(CellError::Num);
+    }
+    let z = if p == 0.0 {
+        0.0
+    } else if p == 1.0 {
+        1.0
+    } else {
+        invert(|x| reg_inc_beta(x, al, be), p, 0.0, 1.0, 1.0)
+    };
+    num_val(lo + z * (hi - lo))
+}
+
+pub(super) fn lognorm_dist(a: &[Arg], _: C) -> R<Value> {
+    let (x, m, sd) = (num(a, 0)?, num(a, 1)?, num(a, 2)?);
+    let cum = if a.len() > 3 { boolean(a, 3)? } else { true };
+    if x <= 0.0 || sd <= 0.0 {
+        return Err(CellError::Num);
+    }
+    let z = (x.ln() - m) / sd;
+    num_val(if cum { norm_cdf(z) } else { norm_pdf(z) / (x * sd) })
+}
+
+pub(super) fn lognorm_inv(a: &[Arg], _: C) -> R<Value> {
+    let (p, m, sd) = (num(a, 0)?, num(a, 1)?, num(a, 2)?);
+    if !(p > 0.0 && p < 1.0) || sd <= 0.0 {
+        return Err(CellError::Num);
+    }
+    num_val((m + sd * norm_inv(p)).exp())
+}
+
+pub(super) fn weibull_dist(a: &[Arg], _: C) -> R<Value> {
+    let (x, al, be, cum) = (num(a, 0)?, num(a, 1)?, num(a, 2)?, boolean(a, 3)?);
+    if x < 0.0 || al <= 0.0 || be <= 0.0 {
+        return Err(CellError::Num);
+    }
+    let t = (x / be).powf(al);
+    num_val(if cum { -(-t).exp_m1() } else { al / be.powf(al) * x.powf(al - 1.0) * (-t).exp() })
+}
+
+pub(super) fn hypgeom_dist(a: &[Arg], _: C) -> R<Value> {
+    let k = num(a, 0)?.trunc();
+    let n = num(a, 1)?.trunc();
+    let kk = num(a, 2)?.trunc();
+    let nn = num(a, 3)?.trunc();
+    let cum = if a.len() > 4 { boolean(a, 4)? } else { false };
+    if k < 0.0 || n <= 0.0 || kk <= 0.0 || nn <= 0.0 || n > nn || kk > nn || k > n.min(kk) || k < (n - nn + kk).max(0.0) {
+        return Err(CellError::Num);
+    }
+    let pmf = |i: f64| (ln_comb(kk, i) + ln_comb(nn - kk, n - i) - ln_comb(nn, n)).exp();
+    if !cum {
+        return num_val(pmf(k));
+    }
+    let mut s = 0.0;
+    let mut i = (n - nn + kk).max(0.0);
+    if k - i > 1e6 || i + 1.0 == i {
+        return Err(CellError::Num);
+    }
+    while i <= k {
+        s += pmf(i);
+        i += 1.0;
+    }
+    num_val(s.min(1.0))
+}
+
+pub(super) fn negbinom_dist(a: &[Arg], _: C) -> R<Value> {
+    let f = num(a, 0)?.trunc();
+    let s = num(a, 1)?.trunc();
+    let p = num(a, 2)?;
+    let cum = if a.len() > 3 { boolean(a, 3)? } else { false };
+    if f < 0.0 || s < 1.0 || !(0.0..=1.0).contains(&p) {
+        return Err(CellError::Num);
+    }
+    if cum {
+        if p == 0.0 {
+            return num_val(0.0);
+        }
+        return num_val(reg_inc_beta(p, s, f + 1.0));
+    }
+    if p == 0.0 || p == 1.0 {
+        return num_val(if p == 1.0 && f == 0.0 { 1.0 } else { 0.0 });
+    }
+    num_val((ln_comb(f + s - 1.0, s - 1.0) + s * p.ln() + f * (1.0 - p).ln()).exp())
+}
+
+pub(super) fn confidence_norm(a: &[Arg], _: C) -> R<Value> {
+    let (al, sd, n) = (num(a, 0)?, num(a, 1)?, num(a, 2)?.trunc());
+    if !(al > 0.0 && al < 1.0) || sd <= 0.0 || n < 1.0 {
+        return Err(CellError::Num);
+    }
+    num_val(norm_inv(1.0 - al / 2.0) * sd / n.sqrt())
+}
+
+pub(super) fn confidence_t(a: &[Arg], _: C) -> R<Value> {
+    let (al, sd, n) = (num(a, 0)?, num(a, 1)?, num(a, 2)?.trunc());
+    if !(al > 0.0 && al < 1.0) || sd <= 0.0 || n < 1.0 {
+        return Err(CellError::Num);
+    }
+    if n == 1.0 {
+        return Err(CellError::Div0);
+    }
+    num_val(t_inv_upper(al / 2.0, n - 1.0) * sd / n.sqrt())
+}
+
+pub(super) fn z_test(a: &[Arg], _: C) -> R<Value> {
+    let v = arg_nums(a, 0)?;
+    let x = num(a, 1)?;
+    if v.is_empty() {
+        return Err(CellError::NA);
+    }
+    let sd = if has(a, 2) { num(a, 2)? } else { var_s(&v)?.sqrt() };
+    if sd <= 0.0 {
+        return Err(if has(a, 2) { CellError::Num } else { CellError::Div0 });
+    }
+    let z = (mean(&v)? - x) / (sd / (v.len() as f64).sqrt());
+    num_val(1.0 - norm_cdf(z))
+}
+
+pub(super) fn fisher(a: &[Arg], _: C) -> R<Value> {
+    let x = num(a, 0)?;
+    if x <= -1.0 || x >= 1.0 {
+        return Err(CellError::Num);
+    }
+    num_val(0.5 * ((1.0 + x) / (1.0 - x)).ln())
+}
+
+pub(super) fn gauss(a: &[Arg], _: C) -> R<Value> {
+    let z = num(a, 0)?;
+    // Φ(z) − 0.5 = erf(z/√2)/2, accurate near zero.
+    num_val(0.5 * erf(z / std::f64::consts::SQRT_2))
+}
diff --git a/crates/functions/src/stat/regression.rs b/crates/functions/src/stat/regression.rs
new file mode 100644
index 000000000..410b086d4
--- /dev/null
+++ b/crates/functions/src/stat/regression.rs
@@ -0,0 +1,342 @@
+use super::descriptive::{arg_pairs, intercept, moments, slope};
+use super::*;
+
+pub(super) fn forecast(a: &[Arg], _: C) -> R<Value> {
+    let x = num(a, 0)?;
+    let (ys, xs) = arg_pairs(a, 1, 2)?;
+    if ys.is_empty() {
+        return Err(CellError::Div0);
+    }
+    num_val(intercept(&ys, &xs)? + slope(&ys, &xs)? * x)
+}
+
+pub(super) fn steyx(a: &[Arg], _: C) -> R<Value> {
+    let (ys, xs) = arg_pairs(a, 0, 1)?;
+    let n = ys.len() as f64;
+    if ys.len() < 3 {
+        return Err(CellError::Div0);
+    }
+    let (sxy, sxx, syy, _, _) = moments(&xs, &ys)?;
+    if sxx == 0.0 {
+        return Err(CellError::Div0);
+    }
+    num_val(((syy - sxy * sxy / sxx) / (n - 2.0)).max(0.0).sqrt())
+}
+
+pub(super) fn covariance(a: &[Arg], sample: bool) -> R<Value> {
+    let (xs, ys) = arg_pairs(a, 0, 1)?;
+    let n = xs.len();
+    if n == 0 || (sample && n < 2) {
+        return Err(CellError::Div0);
+    }
+    let (sxy, ..) = moments(&xs, &ys)?;
+    num_val(sxy / if sample { (n - 1) as f64 } else { n as f64 })
+}
+
+// ---------------------------------------------------------------------------------------------
+// Regression (LINEST, LOGEST, TREND, GROWTH)
+// ---------------------------------------------------------------------------------------------
+
+struct Xy {
+    y: Vec<f64>,
+    /// n rows of k predictor values.
+    x: Vec<Vec<f64>>,
+    k: usize,
+    y_is_row: bool,
+    y_rows: usize,
+    y_cols: usize,
+}
+
+fn strict_num(v: &Value) -> R<f64> {
+    match v {
+        Value::Number(n) => Ok(*n),
+        Value::Error(e) => Err(*e),
+        _ => Err(CellError::Value),
+    }
+}
+
+/// Most predictor columns LINEST, LOGEST, TREND and GROWTH accept. The normal equations are a
+/// k×k matrix inverted in O(k³), so an unbounded k (one `SEQUENCE(1,16384)` argument) allocated
+/// gigabytes and could abort the process.
+pub(super) const MAX_REGRESSORS: usize = 255;
+
+/// Most multiply-adds (rows × k²) a regression may spend building its normal equations.
+const MAX_REGRESSION_WORK: usize = 1_000_000_000;
+
+/// Rejects regressions too large to fit: more than [`MAX_REGRESSORS`] predictors is `#VALUE!`,
+/// more than [`MAX_REGRESSION_WORK`] multiply-adds is `#NUM!`.
+pub(super) fn check_regression_size(n: usize, k: usize) -> R<()> {
+    if k > MAX_REGRESSORS {
+        return Err(CellError::Value);
+    }
+    match n.checked_mul(k).and_then(|nk| nk.checked_mul(k)) {
+        Some(work) if work <= MAX_REGRESSION_WORK => Ok(()),
+        _ => Err(CellError::Num),
+    }
+}
+
+fn parse_xy(y: &Value, x: Option<&Value>) -> R<Xy> {
+    let ya = as_array(y);
+    let n = ya.data.len();
+    let y_is_row = ya.rows == 1 && ya.cols > 1;
+    let yv: Vec<f64> = ya.data.iter().map(strict_num).collect::<R<_>>()?;
+    let (xs, k) = match x {
+        None => ((1..=n).map(|i| vec![i as f64]).collect(), 1),
+        Some(xv) => {
+            let xa = as_array(xv);
+            if xa.rows == ya.rows && xa.cols == ya.cols {
+                (xa.data.iter().map(|v| strict_num(v).map(|f| vec![f])).collect::<R<Vec<_>>>()?, 1)
+            } else if !y_is_row && xa.rows == n {
+                check_regression_size(n, xa.cols)?;
+                let mut rows = Vec::with_capacity(n);
+                for r in 0..n {
+                    rows.push((0..xa.cols).map(|c| strict_num(xa.get(r, c).unwrap_or(&Value::Empty))).collect::<R<Vec<_>>>()?);
+                }
+                (rows, xa.cols)
+            } else if y_is_row && xa.cols == n {
+                check_regression_size(n, xa.rows)?;
+                let mut rows = Vec::with_capacity(n);
+                for c in 0..n {
+                    rows.push((0..xa.rows).map(|r| strict_num(xa.get(r, c).unwrap_or(&Value::Empty))).collect::<R<Vec<_>>>()?);
+                }
+                (rows, xa.rows)
+            } else {
+                return Err(CellError::Ref);
+            }
+        }
+    };
+    Ok(Xy { y: yv, x: xs, k, y_is_row, y_rows: ya.rows, y_cols: ya.cols })
+}
+
+struct Fit {
+    coef: Vec<f64>,
+    b: f64,
+    se: Vec<f64>,
+    se_b: Option<f64>,
+    r2: R<f64>,
+    sey: R<f64>,
+    f: R<f64>,
+    df: f64,
+    ssreg: f64,
+    ssresid: f64,
+}
+
+/// Inverts a symmetric positive semi-definite matrix by Gauss–Jordan on the diagonal, dropping
+/// (zeroing) collinear variables. Returns the inverse and the dropped flags.
+fn sym_inverse(m: &[Vec<f64>]) -> (Vec<Vec<f64>>, Vec<bool>) {
+    let p = m.len();
+    let mut a: Vec<Vec<f64>> = m.to_vec();
+    let mut inv: Vec<Vec<f64>> = (0..p).map(|i| (0..p).map(|j| if i == j { 1.0 } else { 0.0 }).collect()).collect();
+    let mut dropped = vec![false; p];
+    let diag: Vec<f64> = (0..p).map(|i| m.get(i).and_then(|r| r.get(i)).copied().unwrap_or(0.0).abs()).collect();
+    for j in 0..p {
+        let piv = a.get(j).and_then(|r| r.get(j)).copied().unwrap_or(0.0);
+        let scale = diag.get(j).copied().unwrap_or(0.0);
+        if piv.abs() <= 1e-12 * scale || piv.abs() < 1e-300 {
+            if let Some(d) = dropped.get_mut(j) {
+                *d = true;
+            }
+            for row in a.iter_mut() {
+                if let Some(v) = row.get_mut(j) {
+                    *v = 0.0;
+                }
+            }
+            if let Some(r) = a.get_mut(j) {
+                r.iter_mut().for_each(|v| *v = 0.0);
+            }
+            if let Some(r) = inv.get_mut(j) {
+                r.iter_mut().for_each(|v| *v = 0.0);
+            }
+            continue;
+        }
+        let arow: Vec<f64> = a.get(j).map(|r| r.iter().map(|v| v / piv).collect()).unwrap_or_default();
+        let irow: Vec<f64> = inv.get(j).map(|r| r.iter().map(|v| v / piv).collect()).unwrap_or_default();
+        for i in 0..p {
+            if i == j {
+                continue;
+            }
+            let f = a.get(i).and_then(|r| r.get(j)).copied().unwrap_or(0.0);
+            if f == 0.0 {
+                continue;
+            }
+            if let Some(r) = a.get_mut(i) {
+                r.iter_mut().zip(&arow).for_each(|(v, w)| *v -= f * w);
+            }
+            if let Some(r) = inv.get_mut(i) {
+                r.iter_mut().zip(&irow).for_each(|(v, w)| *v -= f * w);
+            }
+        }
+        if let Some(r) = a.get_mut(j) {
+            *r = arow;
+        }
+        if let Some(r) = inv.get_mut(j) {
+            *r = irow;
+        }
+    }
+    (inv, dropped)
+}
+
+fn regress(y: &[f64], x: &[Vec<f64>], k: usize, constant: bool) -> R<Fit> {
+    let n = y.len();
+    if n == 0 || x.len() != n {
+        return Err(CellError::Value);
+    }
+    check_regression_size(n, k)?;
+    let nf = n as f64;
+    let (xm, ym) = if constant {
+        let xm: Vec<f64> = (0..k).map(|j| x.iter().map(|r| r.get(j).copied().unwrap_or(0.0)).sum::<f64>() / nf).collect();
+        (xm, y.iter().sum::<f64>() / nf)
+    } else {
+        (vec![0.0; k], 0.0)
+    };
+    // Normal equations on (centred) data.
+    let mut xtx = vec![vec![0.0; k]; k];
+    let mut xty = vec![0.0; k];
+    for (row, yv) in x.iter().zip(y) {
+        let d: Vec<f64> = (0..k).map(|j| row.get(j).copied().unwrap_or(0.0) - xm.get(j).copied().unwrap_or(0.0)).collect();
+        let yd = yv - ym;
+        for (i, di) in d.iter().enumerate() {
+            if let Some(t) = xty.get_mut(i) {
+                *t += di * yd;
+            }
+            if let Some(r) = xtx.get_mut(i) {
+                for (j, dj) in d.iter().enumerate() {
+                    if let Some(v) = r.get_mut(j) {
+                        *v += di * dj;
+                    }
+                }
+            }
+        }
+    }
+    let (inv, dropped) = sym_inverse(&xtx);
+    let coef: Vec<f64> = (0..k).map(|i| inv.get(i).map(|r| r.iter().zip(&xty).map(|(a, b)| a * b).sum()).unwrap_or(0.0)).collect();
+    let b = if constant { ym - coef.iter().zip(&xm).map(|(c, m)| c * m).sum::<f64>() } else { 0.0 };
+    let mut ssresid = 0.0;
+    let mut sstot = 0.0;
+    for (row, yv) in x.iter().zip(y) {
+        let pred = b + coef.iter().zip(row).map(|(c, v)| c * v).sum::<f64>();
+        ssresid += (yv - pred).powi(2);
+        sstot += if constant { (yv - ym).powi(2) } else { yv * yv };
+    }
+    let used = dropped.iter().filter(|d| !**d).count() + usize::from(constant);
+    let df = nf - used as f64;
+    let ssreg = sstot - ssresid;
+    let r2 = if sstot == 0.0 { Err(CellError::Num) } else { Ok(ssreg / sstot) };
+    let sey2 = if df > 0.0 { Ok(ssresid / df) } else { Err(CellError::Num) };
+    let sey = sey2.map(f64::sqrt);
+    let v1 = if constant { nf - df - 1.0 } else { nf - df };
+    let f = match sey2 {
+        Ok(s2) if s2 > 0.0 && v1 > 0.0 => Ok((ssreg / v1) / s2),
+        _ => Err(CellError::Num),
+    };
+    let s2 = sey2.unwrap_or(f64::NAN);
+    let se: Vec<f64> = (0..k)
+        .map(|i| if dropped.get(i).copied().unwrap_or(true) { 0.0 } else { (inv.get(i).and_then(|r| r.get(i)).copied().unwrap_or(0.0) * s2).sqrt() })
+        .collect();
+    let se_b = if constant {
+        let mut q = 0.0;
+        for (i, mi) in xm.iter().enumerate() {
+            for (j, mj) in xm.iter().enumerate() {
+                q += mi * mj * inv.get(i).and_then(|r| r.get(j)).copied().unwrap_or(0.0);
+            }
+        }
+        Some((s2 * (1.0 / nf + q)).sqrt())
+    } else {
+        None
+    };
+    Ok(Fit { coef, b, se, se_b, r2, sey, f, df, ssreg, ssresid })
+}
+
+fn const_arg(a: &[Arg], i: usize) -> R<bool> {
+    if has(a, i) { boolean(a, i) } else { Ok(true) }
+}
+
+fn opt_value(a: &[Arg], i: usize) -> Option<&Value> {
+    if has(a, i) { a.get(i).map(|x| &x.value) } else { None }
+}
+
+pub(super) fn linest_impl(a: &[Arg], log: bool) -> R<Value> {
+    let mut xy = parse_xy(&arg(a, 0)?.value, opt_value(a, 1))?;
+    if log {
+        if xy.y.iter().any(|&v| v <= 0.0) {
+            return Err(CellError::Num);
+        }
+        xy.y.iter_mut().for_each(|v| *v = v.ln());
+    }
+    let constant = const_arg(a, 2)?;
+    let stats = if has(a, 3) { boolean(a, 3)? } else { false };
+    let fit = regress(&xy.y, &xy.x, xy.k, constant)?;
+    let tr = |v: f64| if log { v.exp() } else { v };
+    let cols = xy.k + 1;
+    let mut data: Vec<Value> = Vec::new();
+    for c in fit.coef.iter().rev() {
+        data.push(Value::number(tr(*c)));
+    }
+    data.push(Value::number(tr(fit.b)));
+    if stats {
+        let rv = |r: R<f64>| r.map(Value::number).unwrap_or_else(Value::Error);
+        let na = Value::Error(CellError::NA);
+        for s in fit.se.iter().rev() {
+            data.push(Value::number(*s));
+        }
+        data.push(fit.se_b.map(Value::number).unwrap_or(na.clone()));
+        let rows: [[Value; 2]; 3] =
+            [[rv(fit.r2), rv(fit.sey)], [rv(fit.f), Value::number(fit.df)], [Value::number(fit.ssreg), Value::number(fit.ssresid)]];
+        for [p, q] in rows {
+            data.push(p);
+            if cols > 1 {
+                data.push(q);
+            }
+            for _ in 2..cols {
+                data.push(na.clone());
+            }
+        }
+        return array_val(5, cols, data);
+    }
+    array_val(1, cols, data)
+}
+
+pub(super) fn trend_impl(a: &[Arg], log: bool) -> R<Value> {
+    let mut xy = parse_xy(&arg(a, 0)?.value, opt_value(a, 1))?;
+    if log {
+        if xy.y.iter().any(|&v| v <= 0.0) {
+            return Err(CellError::Num);
+        }
+        xy.y.iter_mut().for_each(|v| *v = v.ln());
+    }
+    let constant = const_arg(a, 3)?;
+    let fit = regress(&xy.y, &xy.x, xy.k, constant)?;
+    let predict = |row: &[f64]| -> Value {
+        let v = fit.b + fit.coef.iter().zip(row).map(|(c, x)| c * x).sum::<f64>();
+        Value::number(if log { v.exp() } else { v })
+    };
+    match opt_value(a, 2) {
+        None => {
+            let data: Vec<Value> = xy.x.iter().map(|r| predict(r)).collect();
+            array_val(xy.y_rows, xy.y_cols, data)
+        }
+        Some(nv) => {
+            let na = as_array(nv);
+            if xy.k == 1 {
+                let data = na.data.iter().map(|v| strict_num(v).map(|x| predict(&[x]))).collect::<R<Vec<_>>>()?;
+                array_val(na.rows, na.cols, data)
+            } else if !xy.y_is_row && na.cols == xy.k {
+                let mut data = Vec::with_capacity(na.rows);
+                for r in 0..na.rows {
+                    let row = (0..na.cols).map(|c| strict_num(na.get(r, c).unwrap_or(&Value::Empty))).collect::<R<Vec<_>>>()?;
+                    data.push(predict(&row));
+                }
+                array_val(na.rows, 1, data)
+            } else if xy.y_is_row && na.rows == xy.k {
+                let mut data = Vec::with_capacity(na.cols);
+                for c in 0..na.cols {
+                    let row = (0..na.rows).map(|r| strict_num(na.get(r, c).unwrap_or(&Value::Empty))).collect::<R<Vec<_>>>()?;
+                    data.push(predict(&row));
+                }
+                array_val(1, na.cols, data)
+            } else {
+                Err(CellError::Ref)
+            }
+        }
+    }
+}