From 5c383ae1854198c06db2ca8b0f43215b5b8e20fa Mon Sep 17 00:00:00 2001 From: Felix Agene Date: Sun, 26 Jul 2026 09:49:34 -0500 Subject: [PATCH 1/3] fix: replace hand-rolled quickselect with select_nth_unstable_by `select_inplace` implemented the Numerical Recipes quickselect by hand. With NaN in the slice its partition loop walks past the end of the array, which is the crash reported in #163 (a segfault under the older unsafe version, an index-out-of-bounds panic today): let mut v: Vec = (0..1000).map(|i| (i as f64 * 0.7).sin()).collect(); for i in (0..1000).step_by(7) { v[i] = f64::NAN; } Data::new(v).quantile(0.5); // panics on main `<[T]>::select_nth_unstable_by` with `total_cmp` fixes this three ways: NaN gets a defined position in a total order, the worst case is O(n) rather than quadratic, and it is ~4.4x faster at 1e6 elements. Closes #163 --- src/statistics/slice_statistics.rs | 98 +++++++++++++----------------- 1 file changed, 41 insertions(+), 57 deletions(-) diff --git a/src/statistics/slice_statistics.rs b/src/statistics/slice_statistics.rs index 96b83e7b..28d00641 100644 --- a/src/statistics/slice_statistics.rs +++ b/src/statistics/slice_statistics.rs @@ -64,8 +64,12 @@ impl + AsRef<[f64]>> Data { self.0.as_ref().iter() } - // Selection algorithm from Numerical Recipes - // See: https://en.wikipedia.org/wiki/Selection_algorithm + // Selection via `<[T]>::select_nth_unstable_by`, which is an introselect + // with an O(n) worst case. The previous hand-rolled Numerical Recipes + // quickselect was ~4x slower on large uniform data and quadratic on + // adversarial inputs; `total_cmp` also gives NaN a defined ordering + // (positive NaNs sort past +inf) instead of the arbitrary result the NR + // partition produced. fn select_inplace(&mut self, rank: usize) -> f64 { if rank == 0 { return self.min(); @@ -74,61 +78,11 @@ impl + AsRef<[f64]>> Data { return self.max(); } - let mut low = 0; - let mut high = self.len() - 1; - loop { - if high <= low + 1 { - if high == low + 1 && self[high] < self[low] { - self.swap(low, high) - } - return self[rank]; - } - - let middle = (low + high) / 2; - self.swap(middle, low + 1); - - if self[low] > self[high] { - self.swap(low, high); - } - if self[low + 1] > self[high] { - self.swap(low + 1, high); - } - if self[low] > self[low + 1] { - self.swap(low, low + 1); - } - - let mut begin = low + 1; - let mut end = high; - let pivot = self[begin]; - loop { - loop { - begin += 1; - if self[begin] >= pivot { - break; - } - } - loop { - end -= 1; - if self[end] <= pivot { - break; - } - } - if end < begin { - break; - } - self.swap(begin, end); - } - - self[low + 1] = self[end]; - self[end] = pivot; - - if end >= rank { - high = end - 1; - } - if end <= rank { - low = begin; - } - } + *self + .0 + .as_mut() + .select_nth_unstable_by(rank, |a, b| a.total_cmp(b)) + .1 } } @@ -550,6 +504,36 @@ mod tests { // TODO: test codeplex issue 5667 (Math.NET) + /// `select_inplace` used to implement the Numerical Recipes quickselect by + /// hand, whose partition loop walks off the end of the array when the slice + /// contains NaN (statrs-dev/statrs#163). `total_cmp` gives NaN a defined + /// place in a total order, so the search stays in bounds. + #[test] + fn test_order_statistics_with_nan_do_not_panic() { + let mut v: Vec = (0..1000).map(|i| (i as f64 * 0.7).sin()).collect(); + for i in (0..1000).step_by(7) { + v[i] = f64::NAN; + } + let mut d = Data::new(v); + // the values themselves are unspecified once NaN is present; what is + // being pinned is that these terminate in bounds + let _ = d.quantile(0.5); + let _ = d.quantile(0.0); + let _ = d.quantile(1.0); + let _ = d.order_statistic(500); + let _ = OrderStatistics::median(&mut d); + let _ = d.interquartile_range(); + + // all-NaN, and NaN at the exact positions the old partition tripped on + let mut all_nan = Data::new(vec![f64::NAN; 64]); + let _ = all_nan.quantile(0.5); + for pos in [0usize, 1, 31, 62, 63] { + let mut v = vec![0.0f64; 64]; + v[pos] = f64::NAN; + let _ = Data::new(v).quantile(0.5); + } + } + #[test] fn test_median_robust_on_infinities() { let data3 = [2.0, f64::NEG_INFINITY, f64::INFINITY]; From 79943e2f18744780179ce46573cec04682411914 Mon Sep 17 00:00:00 2001 From: Felix Agene Date: Sun, 26 Jul 2026 18:43:16 -0500 Subject: [PATCH 2/3] test: build the NaN order-statistic test under no_std The test allocated its input with `Vec`, which does not exist under `--no-default-features`: the crate is `no_std` without `alloc`, and there is no `extern crate alloc`. `core::array::from_fn` and array literals give the same data with no allocation. CI did not catch this because the test job runs only default features, and the feature-powerset job passes `--no-dev-deps`, so test code is never compiled without std. Co-Authored-By: Claude Opus 5 (1M context) --- src/statistics/slice_statistics.rs | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/statistics/slice_statistics.rs b/src/statistics/slice_statistics.rs index 28d00641..41e49ace 100644 --- a/src/statistics/slice_statistics.rs +++ b/src/statistics/slice_statistics.rs @@ -510,7 +510,8 @@ mod tests { /// place in a total order, so the search stays in bounds. #[test] fn test_order_statistics_with_nan_do_not_panic() { - let mut v: Vec = (0..1000).map(|i| (i as f64 * 0.7).sin()).collect(); + // Arrays rather than `Vec`, so the test also builds under `no_std`. + let mut v: [f64; 1000] = core::array::from_fn(|i| (i as f64 * 0.7).sin()); for i in (0..1000).step_by(7) { v[i] = f64::NAN; } @@ -525,10 +526,10 @@ mod tests { let _ = d.interquartile_range(); // all-NaN, and NaN at the exact positions the old partition tripped on - let mut all_nan = Data::new(vec![f64::NAN; 64]); + let mut all_nan = Data::new([f64::NAN; 64]); let _ = all_nan.quantile(0.5); for pos in [0usize, 1, 31, 62, 63] { - let mut v = vec![0.0f64; 64]; + let mut v = [0.0f64; 64]; v[pos] = f64::NAN; let _ = Data::new(v).quantile(0.5); } From 3ff51d7439d5c932e4d0fe5e9db9d69d842b1f0e Mon Sep 17 00:00:00 2001 From: Felix Agene Date: Mon, 27 Jul 2026 12:07:43 -0500 Subject: [PATCH 3/3] perf: add a selection benchmark, and ordered fast paths it exposed Follow-up to @day01's review on #407. Adds a `selection scaling` group to benches/order_statistics.rs sweeping size against input shape, written against only the public `Data` API so the same file compiles on any branch and `--save-baseline` works across a checkout. The existing group measures one fixed size on random data, which is not enough to reproduce a performance claim. Running it against main showed the claim in this PR was wrong. The shapes are why: the Numerical Recipes quickselect took its pivot from a median of three, which on sorted or reverse-sorted input lands on the true median and finishes in one partition. It beat plain introselect on exactly those shapes -- by 75 to 89% on sorted input at every size measured -- while losing badly on large random, organ-pipe and duplicate-heavy input. The honest summary of the original change was "much faster on some shapes, nearly twice as slow on others", not "~4x faster". So test for order first, which recovers what median-of-three was doing and improves on it: the answer becomes an index rather than a partition. Both checks stop at the first out-of-order pair, so unordered input pays a couple of comparisons, and `median` and `quantile` benefit twice since they select twice. `<=` is used rather than `total_cmp`, which is cheaper per element and still correct -- any NaN makes some comparison false, so a slice containing one never takes a fast path. With that, 14 of 15 cases beat main (medians of the Criterion intervals): shape 1024 65536 1048576 random +19% -28% -76% sorted -34% -30% -24% reversed -56% -51% -52% organ_pipe -46% -65% -71% duplicates -15% -63% -78% The exception is small random input, where introselect's setup costs more than the old quickselect's; 2.0us against 2.4us at n = 1024. Left alone, since the reason for this PR is that the NR partition walks off the end of the array on NaN, not throughput. Also adds the positive-value tests asked for. The existing ones cover NaN, infinities and constant sequences -- they pin that nothing panics but say nothing about the values. The new ones check every order statistic, the median, and the quantile endpoints against a sorted reference over seven shapes and nine sizes, plus permutation invariance. These cover the new fast paths, which return an index rather than partitioning: mutating the descending path to index from the wrong end fails both. Co-Authored-By: Claude Opus 5 (1M context) --- benches/order_statistics.rs | 63 ++++++++++++- src/statistics/slice_statistics.rs | 141 ++++++++++++++++++++++++++--- 2 files changed, 192 insertions(+), 12 deletions(-) diff --git a/benches/order_statistics.rs b/benches/order_statistics.rs index d1771ff4..c103143a 100644 --- a/benches/order_statistics.rs +++ b/benches/order_statistics.rs @@ -94,5 +94,66 @@ fn bench_order_statistic(c: &mut Criterion) { group.finish(); } -criterion_group!(benches, bench_order_statistic); +/// Selection cost across input sizes and shapes. +/// +/// Kept separate from the fixed-size group above so the two can evolve +/// independently, and written against only the public `Data` API so the same +/// file compiles on any branch. To compare two implementations: +/// +/// ```text +/// git checkout main && cargo bench --bench order_statistics -- --save-baseline main +/// git checkout && cargo bench --bench order_statistics -- --baseline main +/// ``` +/// +/// The shapes matter as much as the sizes. A median-of-three pivot is fine on +/// random data and degrades on input built to defeat it, so a random-only +/// benchmark would understate the difference between a quickselect and an +/// introselect with an O(n) fallback. +fn bench_selection_scaling(c: &mut Criterion) { + fn shaped(shape: &str, n: usize) -> Vec { + let mut rng = StdRng::seed_from_u64(0xB0A7); + match shape { + // shuffled, the ordinary case + "random" => { + let mut v: Vec = (0..n).map(|i| i as f64).collect(); + v.shuffle(&mut rng); + v + } + // already ordered, both directions + "sorted" => (0..n).map(|i| i as f64).collect(), + "reversed" => (0..n).rev().map(|i| i as f64).collect(), + // ascends then descends: the classic median-of-three adversary, + // since the first, middle and last elements are unrepresentative + "organ_pipe" => (0..n) + .map(|i| if i < n / 2 { i } else { n - i } as f64) + .collect(), + // few distinct values, so partitions are heavily unbalanced + "duplicates" => { + let distinct = ((n as f64).sqrt() as usize).max(1); + let mut v: Vec = (0..n).map(|i| (i % distinct) as f64).collect(); + v.shuffle(&mut rng); + v + } + other => panic!("unknown shape {other}"), + } + } + + let mut group = c.benchmark_group("selection scaling"); + for &n in &[1_024usize, 65_536, 1_048_576] { + for shape in ["random", "sorted", "reversed", "organ_pipe", "duplicates"] { + let data = shaped(shape, n); + group.throughput(Throughput::Elements(n as u64)); + group.bench_function(format!("median/{shape}/{n}"), |b| { + b.iter_batched( + || Data::new(data.clone()), + |data| black_box(data.median()), + BatchSize::LargeInput, + ) + }); + } + } + group.finish(); +} + +criterion_group!(benches, bench_order_statistic, bench_selection_scaling); criterion_main!(benches); diff --git a/src/statistics/slice_statistics.rs b/src/statistics/slice_statistics.rs index 41e49ace..62840b25 100644 --- a/src/statistics/slice_statistics.rs +++ b/src/statistics/slice_statistics.rs @@ -64,12 +64,19 @@ impl + AsRef<[f64]>> Data { self.0.as_ref().iter() } - // Selection via `<[T]>::select_nth_unstable_by`, which is an introselect - // with an O(n) worst case. The previous hand-rolled Numerical Recipes - // quickselect was ~4x slower on large uniform data and quadratic on - // adversarial inputs; `total_cmp` also gives NaN a defined ordering - // (positive NaNs sort past +inf) instead of the arbitrary result the NR - // partition produced. + // Selection via `<[T]>::select_nth_unstable_by`, an introselect with an + // O(n) worst case. `total_cmp` also gives NaN a defined position (positive + // NaNs sort past +inf) rather than the arbitrary result the previous + // Numerical Recipes partition produced, which is the point of the change. + // + // The ordered fast paths are not incidental. The NR quickselect took its + // pivot from a median of three, which on sorted or reverse-sorted input + // lands on the true median immediately and finishes in a single partition + // -- so it beat plain introselect on exactly those shapes. Testing for + // order first recovers that, and improves on it: the answer becomes an + // index rather than a partition. Both checks stop at the first out-of-order + // pair, costing a couple of comparisons on unordered input, and they also + // pay off in `median` and `quantile`, which call this twice. fn select_inplace(&mut self, rank: usize) -> f64 { if rank == 0 { return self.min(); @@ -78,11 +85,19 @@ impl + AsRef<[f64]>> Data { return self.max(); } - *self - .0 - .as_mut() - .select_nth_unstable_by(rank, |a, b| a.total_cmp(b)) - .1 + // `<=` rather than `total_cmp` here, which is both cheaper per element + // and still correct: any NaN makes some comparison false, so a slice + // containing one never takes a fast path and falls through to the + // general case, and without NaN the two orderings agree. + let slice = self.0.as_mut(); + if slice.is_sorted_by(|a, b| a <= b) { + return slice[rank]; + } + if slice.is_sorted_by(|a, b| a >= b) { + return slice[slice.len() - 1 - rank]; + } + + *slice.select_nth_unstable_by(rank, |a, b| a.total_cmp(b)).1 } } @@ -504,6 +519,110 @@ mod tests { // TODO: test codeplex issue 5667 (Math.NET) + /// Selection agrees with sorting, on ordinary values across a range of + /// shapes and sizes. + /// + /// The counterpart to the NaN and infinity tests below, which pin that + /// nothing panics but say nothing about the values returned. This also + /// covers the ordered fast paths in `select_inplace`: those return an + /// index instead of partitioning, so they have to produce exactly what the + /// general path would, including for the descending case where the index + /// is counted from the far end. + #[test] + #[cfg(feature = "std")] + fn test_order_statistics_match_sorted_reference() { + fn shaped(shape: &str, n: usize) -> Vec { + match shape { + // deterministic pseudo-random, no rand dependency needed + "scattered" => (0..n).map(|i| ((i * 7919 + 13) % 1000) as f64 * 0.5).collect(), + "sorted" => (0..n).map(|i| i as f64 * 1.5).collect(), + "reversed" => (0..n).rev().map(|i| i as f64 * 1.5).collect(), + "organ_pipe" => (0..n).map(|i| if i < n / 2 { i } else { n - i } as f64).collect(), + "duplicates" => (0..n).map(|i| (i % 3) as f64).collect(), + "constant" => vec![4.25; n], + "negatives" => (0..n).map(|i| -((i * 37 % 100) as f64)).collect(), + other => panic!("unknown shape {other}"), + } + } + + for n in [1usize, 2, 3, 4, 5, 8, 17, 101, 1000] { + for shape in [ + "scattered", "sorted", "reversed", "organ_pipe", "duplicates", "constant", + "negatives", + ] { + let data = shaped(shape, n); + let mut sorted = data.clone(); + sorted.sort_by(f64::total_cmp); + + // every order statistic, 1-based + for order in 1..=n { + let mut d = Data::new(data.clone()); + assert_eq!( + d.order_statistic(order), + sorted[order - 1], + "{shape}/{n}: order_statistic({order})" + ); + } + + // out of range on either side + let mut d = Data::new(data.clone()); + assert!(d.order_statistic(0).is_nan(), "{shape}/{n}: order 0"); + assert!(d.order_statistic(n + 1).is_nan(), "{shape}/{n}: order n+1"); + + // median, including the even case that averages two selections + let expected_median = if n % 2 == 1 { + sorted[n / 2] + } else { + (sorted[n / 2 - 1] + sorted[n / 2]) / 2.0 + }; + let mut d = Data::new(data.clone()); + assert_eq!( + OrderStatistics::median(&mut d), + expected_median, + "{shape}/{n}: median" + ); + + // the quantile endpoints, and that interior values stay in range + let mut d = Data::new(data.clone()); + assert_eq!(d.quantile(0.0), sorted[0], "{shape}/{n}: quantile(0)"); + assert_eq!(d.quantile(1.0), sorted[n - 1], "{shape}/{n}: quantile(1)"); + for &tau in &[0.1, 0.25, 0.5, 0.75, 0.9] { + let q = d.quantile(tau); + assert!( + q >= sorted[0] && q <= sorted[n - 1], + "{shape}/{n}: quantile({tau}) = {q} outside [{}, {}]", + sorted[0], + sorted[n - 1] + ); + } + } + } + } + + /// Selection must not depend on the order the input happens to arrive in. + #[test] + #[cfg(feature = "std")] + fn test_order_statistic_is_permutation_invariant() { + let base: Vec = (0..64).map(|i| ((i * 31 + 7) % 41) as f64).collect(); + let mut sorted = base.clone(); + sorted.sort_by(f64::total_cmp); + + // a few fixed permutations, including the ones that hit the fast paths + let mut rotated = base.clone(); + rotated.rotate_left(23); + let mut ascending = base.clone(); + ascending.sort_by(f64::total_cmp); + let mut descending = ascending.clone(); + descending.reverse(); + + for variant in [base.clone(), rotated, ascending, descending] { + for order in 1..=variant.len() { + let mut d = Data::new(variant.clone()); + assert_eq!(d.order_statistic(order), sorted[order - 1]); + } + } + } + /// `select_inplace` used to implement the Numerical Recipes quickselect by /// hand, whose partition loop walks off the end of the array when the slice /// contains NaN (statrs-dev/statrs#163). `total_cmp` gives NaN a defined