diff --git a/crates/functions/src/stat.rs b/crates/functions/src/stat.rs
--- a/crates/functions/src/stat.rs
+++ b/crates/functions/src/stat.rs
@@ -51,10 +51,11 @@ 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)) }
+    let mut v = v;
+    let n = v.len();
+    let (lower, hi, _) = v.select_nth_unstable_by(n / 2, f64::total_cmp);
+    let hi = *hi;
+    if n % 2 == 1 { Ok(hi) } else { Ok(lower.iter().copied().max_by(f64::total_cmp).unwrap_or(0.0).midpoint(hi)) }
 }
 
 /// Values occurring at least twice, ordered by first occurrence, with their (equal, maximal) count.
@@ -125,29 +126,38 @@ fn kth(v: Vec<f64>, k: f64, largest: bool) -> R<f64> {
     if k.is_nan() || k < 1.0 || k > v.len() as f64 {
         return Err(CellError::Num);
     }
-    let s = sorted(v);
+    let mut v = v;
     let k = k as usize;
-    let idx = if largest { s.len() - k } else { k - 1 };
-    s.get(idx).copied().ok_or(CellError::Num)
+    let idx = if largest { v.len() - k } else { k - 1 };
+    // Selection instead of a full sort; same total order, so the same bit pattern as `sorted(v)[idx]`.
+    Ok(*v.select_nth_unstable_by(idx, f64::total_cmp).1)
 }
 
-fn interp(s: &[f64], pos: f64) -> R<f64> {
+/// Value at fractional rank `pos` of the sorted data, linearly interpolated, found by selection.
+fn interp(mut v: Vec<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),
+    if i >= v.len() {
+        return Err(CellError::Num);
+    }
+    let (_, lo, upper) = v.select_nth_unstable_by(i, f64::total_cmp);
+    let lo = *lo;
+    // The next order statistic is the minimum of the upper partition.
+    if frac > 0.0
+        && let Some(hi) = upper.iter().copied().min_by(f64::total_cmp)
+    {
+        return 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)
+    let pos = k * (v.len() - 1) as f64;
+    interp(v, pos)
 }
 
 fn percentile_exc(v: Vec<f64>, k: f64) -> R<f64> {
@@ -159,8 +169,7 @@ fn percentile_exc(v: Vec<f64>, k: f64) -> R<f64> {
     if pos < 1.0 || pos > n {
         return Err(CellError::Num);
     }
-    let s = sorted(v);
-    interp(&s, pos - 1.0)
+    interp(v, pos - 1.0)
 }
 
 fn quartile_inc(v: Vec<f64>, q: f64) -> R<f64> {
@@ -2274,7 +2283,9 @@ pub(crate) fn specs() -> Vec<FnSpec> {
 
 #[cfg(test)]
 mod tests {
-    use super::{MAX_REGRESSORS, aggregate_values, aggregate_values_k, check_regression_size};
+    use super::{
+        MAX_REGRESSORS, R, aggregate_values, aggregate_values_k, check_regression_size, kth, median, percentile_exc, percentile_inc, sorted,
+    };
     use crate::util::testutil::*;
     use gridcraft_core::{CellError, Value};
 
@@ -2463,6 +2474,110 @@ mod tests {
         close(ev("RANK.AVG", vec![n(3.5), r2, n(1.0)]), 3.5);
     }
 
+    // Reference implementations: the previous full-sort versions, as the oracle for the selection ones.
+    fn old_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 old_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)) }
+    }
+
+    fn old_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 old_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);
+        old_interp(&s, k * (s.len() - 1) as f64)
+    }
+
+    fn old_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);
+        old_interp(&s, pos - 1.0)
+    }
+
+    /// Same Ok bit pattern or same error.
+    fn same(a: R<f64>, b: R<f64>) -> bool {
+        match (a, b) {
+            (Ok(x), Ok(y)) => x.to_bits() == y.to_bits(),
+            (Err(x), Err(y)) => x == y,
+            _ => false,
+        }
+    }
+
+    #[test]
+    fn selection_matches_full_sort_bit_for_bit() {
+        let mut state = 0x9E37_79B9_7F4A_7C15_u64;
+        let mut next = move || {
+            state ^= state << 13;
+            state ^= state >> 7;
+            state ^= state << 17;
+            state
+        };
+        let specials = [0.0, -0.0, f64::INFINITY, f64::NEG_INFINITY, f64::NAN, -f64::NAN, f64::MAX, f64::MIN_POSITIVE, 1.0, -1.0];
+        for len in (0..=40).chain([63, 64, 65, 100, 257, 1000, 4096]) {
+            for round in 0..12 {
+                let v: Vec<f64> = (0..len)
+                    .map(|_| match (round % 4, next() % 8) {
+                        // Few distinct values: many ties.
+                        (0, r) => (r % 3) as f64,
+                        // Specials, including both zeros and NaNs.
+                        (1, r) => specials[(next() as usize + r as usize) % specials.len()],
+                        // Wide random range.
+                        (2, _) => (next() as i64 as f64) / 1e3,
+                        // Integers with some specials mixed in.
+                        (_, 0) => specials[next() as usize % specials.len()],
+                        _ => (next() % 50) as f64,
+                    })
+                    .collect();
+                assert!(same(median(v.clone()), old_median(v.clone())), "median len={len}");
+                let ks = [-1.0, 0.0, 0.2, 1.0, 1.5, 2.0, len as f64 / 2.0, len as f64 - 0.5, len as f64, len as f64 + 1.0, f64::NAN];
+                for &k in &ks {
+                    for largest in [true, false] {
+                        assert!(same(kth(v.clone(), k, largest), old_kth(v.clone(), k, largest)), "kth len={len} k={k}");
+                    }
+                }
+                let ps = [-0.1, 0.0, 0.01, 0.1, 0.25, 0.3, 0.5, 0.75, 0.9, 0.99, 1.0, 1.1, f64::NAN, (next() % 1001) as f64 / 1000.0];
+                for &p in &ps {
+                    assert!(same(percentile_inc(v.clone(), p), old_percentile_inc(v.clone(), p)), "pinc len={len} p={p}");
+                    assert!(same(percentile_exc(v.clone(), p), old_percentile_exc(v.clone(), p)), "pexc len={len} p={p}");
+                }
+            }
+        }
+    }
+
     #[test]
     fn percentiles() {
         close(ev("PERCENTILE.INC", vec![c(&[1.0, 3.0, 2.0, 4.0]), n(0.3)]), 1.9);