From 890870cd57167fda190716689b64897fcb16596e Mon Sep 17 00:00:00 2001 From: vangel76 Date: Mon, 26 Oct 2026 22:11:08 +0100 Subject: [PATCH 1/3] perf: select instead of sort in MEDIAN, PERCENTILE, QUARTILE, LARGE, SMALL Read the needed order statistics with select_nth_unstable_by(total_cmp) instead of a full sort. Results are bit-identical; a test compares them against the previous full-sort code on randomised data. Fixes #2679 Co-Authored-By: Claude Sonnet 5.5 --- crates/functions/src/stat.rs | 149 +++++++++++++++++++++++++++++++---- 1 file changed, 132 insertions(+), 17 deletions(-) diff --git a/crates/functions/src/stat.rs b/crates/functions/src/stat.rs index 344ba9ff8..4ef2625cc 100644 --- a/crates/functions/src/stat.rs +++ b/crates/functions/src/stat.rs @@ -51,10 +51,11 @@ fn median(v: Vec) -> R { 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, k: f64, largest: bool) -> R { 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 { +/// Value at fractional rank `pos` of the sorted data, linearly interpolated, found by selection. +fn interp(mut v: Vec, pos: f64) -> R { 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, k: f64) -> R { 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, k: f64) -> R { @@ -159,8 +169,7 @@ fn percentile_exc(v: Vec, k: f64) -> R { 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, q: f64) -> R { @@ -2274,7 +2283,9 @@ pub(crate) fn specs() -> Vec { #[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 { + 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) -> R { + 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, k: f64, largest: bool) -> R { + 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, k: f64) -> R { + 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, k: f64) -> R { + 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, b: R) -> 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 = (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); From 4523b817ce53be54a9ce9d5b5b49761c3907a75d Mon Sep 17 00:00:00 2001 From: vangel76 Date: Thu, 29 Oct 2026 00:51:12 +0100 Subject: [PATCH 2/3] Read already-sorted input directly instead of selecting A full sort recognises ascending and descending runs in O(n), which a selection cannot beat, so sorted, reverse-sorted and tie-block input was up to 1.8x slower than before. Check for a monotone run first (one pass that stops at the first out-of-order pair) and index into it. Extend the bit-for-bit oracle test with ascending, descending and nearly ascending data. Co-Authored-By: Claude Sonnet 5.5 --- crates/functions/src/stat/descriptive.rs | 32 +++++++++++++++++++++++- crates/functions/src/stat/tests.rs | 22 +++++++++++++--- 2 files changed, 50 insertions(+), 4 deletions(-) diff --git a/crates/functions/src/stat/descriptive.rs b/crates/functions/src/stat/descriptive.rs index 49d30d42c..08938d039 100644 --- a/crates/functions/src/stat/descriptive.rs +++ b/crates/functions/src/stat/descriptive.rs @@ -38,12 +38,31 @@ pub(super) fn sorted(mut v: Vec) -> Vec { v } +/// Sorts `v` if it is already monotone in `total_cmp` order (a descending slice is reversed) and +/// says whether it did. One cheap pass that stops at the first out-of-order pair, so unordered +/// data costs next to nothing. Needed because the std sort recognises ascending and descending +/// runs in O(n), which a selection cannot beat. +fn sorted_if_monotone(v: &mut [f64]) -> bool { + if v.is_sorted_by(|a, b| a.total_cmp(b).is_le()) { + return true; + } + if v.is_sorted_by(|a, b| a.total_cmp(b).is_ge()) { + v.reverse(); + return true; + } + false +} + pub(super) fn median(v: Vec) -> R { if v.is_empty() { return Err(CellError::Num); } let mut v = v; let n = v.len(); + if sorted_if_monotone(&mut v) { + let hi = v.get(n / 2).copied().unwrap_or(0.0); + return if n % 2 == 1 { Ok(hi) } else { Ok(v.get(n / 2 - 1).copied().unwrap_or(0.0).midpoint(hi)) }; + } 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)) } @@ -120,7 +139,11 @@ pub(super) fn kth(v: Vec, k: f64, largest: bool) -> R { let mut v = v; let k = k as usize; 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]`. + // Selection instead of a full sort (monotone input is read directly); same total order, so the + // same bit pattern as `sorted(v)[idx]`. + if sorted_if_monotone(&mut v) { + return v.get(idx).copied().ok_or(CellError::Num); + } Ok(*v.select_nth_unstable_by(idx, f64::total_cmp).1) } @@ -132,6 +155,13 @@ fn interp(mut v: Vec, pos: f64) -> R { if i >= v.len() { return Err(CellError::Num); } + if sorted_if_monotone(&mut v) { + let lo = v.get(i).copied().ok_or(CellError::Num)?; + return match v.get(i + 1) { + Some(hi) if frac > 0.0 => Ok(lo + frac * (hi - lo)), + _ => Ok(lo), + }; + } 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. diff --git a/crates/functions/src/stat/tests.rs b/crates/functions/src/stat/tests.rs index 60da9f59e..f11629369 100644 --- a/crates/functions/src/stat/tests.rs +++ b/crates/functions/src/stat/tests.rs @@ -263,8 +263,8 @@ fn selection_matches_full_sort_bit_for_bit() { }; 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 = (0..len) + for round in 0..24 { + let mut v: Vec = (0..len) .map(|_| match (round % 4, next() % 8) { // Few distinct values: many ties. (0, r) => (r % 3) as f64, @@ -277,7 +277,23 @@ fn selection_matches_full_sort_bit_for_bit() { _ => (next() % 50) as f64, }) .collect(); - assert!(same(median(v.clone()), old_median(v.clone())), "median len={len}"); + // Rounds 12..24 repeat the data as already ascending, already descending and nearly + // ascending (ends swapped): the monotone shortcut and its fallback to selection. + match round / 4 { + 3 => v.sort_by(f64::total_cmp), + 4 => { + v.sort_by(f64::total_cmp); + v.reverse(); + } + 5 => { + v.sort_by(f64::total_cmp); + if let Some(last) = v.len().checked_sub(1) { + v.swap(0, last); + } + } + _ => {} + } + assert!(same(median(v.clone()), old_median(v.clone())), "median len={len} round={round}"); 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] { From 8b330ac9e213c05906b4dee0468444e558f6e4a0 Mon Sep 17 00:00:00 2001 From: vangel76 Date: Thu, 29 Oct 2026 00:52:59 +0100 Subject: [PATCH 3/3] Tighten selection comment and test index arithmetic Co-Authored-By: Claude Sonnet 5.5 --- crates/functions/src/stat/descriptive.rs | 7 ++++--- crates/functions/src/stat/tests.rs | 2 +- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/crates/functions/src/stat/descriptive.rs b/crates/functions/src/stat/descriptive.rs index 08938d039..3c8f78f69 100644 --- a/crates/functions/src/stat/descriptive.rs +++ b/crates/functions/src/stat/descriptive.rs @@ -39,9 +39,10 @@ pub(super) fn sorted(mut v: Vec) -> Vec { } /// Sorts `v` if it is already monotone in `total_cmp` order (a descending slice is reversed) and -/// says whether it did. One cheap pass that stops at the first out-of-order pair, so unordered -/// data costs next to nothing. Needed because the std sort recognises ascending and descending -/// runs in O(n), which a selection cannot beat. +/// says whether it did. At most two passes, each stopping at the first out-of-order pair, so +/// unordered data costs next to nothing. Needed because the std sort recognises ascending and +/// descending runs in O(n) with a small constant; selection is O(n) too but measured 1.2-1.8x +/// slower on such input. fn sorted_if_monotone(v: &mut [f64]) -> bool { if v.is_sorted_by(|a, b| a.total_cmp(b).is_le()) { return true; diff --git a/crates/functions/src/stat/tests.rs b/crates/functions/src/stat/tests.rs index f11629369..023b36815 100644 --- a/crates/functions/src/stat/tests.rs +++ b/crates/functions/src/stat/tests.rs @@ -269,7 +269,7 @@ fn selection_matches_full_sort_bit_for_bit() { // 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()], + (1, r) => specials[(next() as usize).wrapping_add(r as usize) % specials.len()], // Wide random range. (2, _) => (next() as i64 as f64) / 1e3, // Integers with some specials mixed in.