From 6be4cfd512ec609b57ac068dbd763f5d9c1248e1 Mon Sep 17 00:00:00 2001
From: Rouven Spreckels <rs@qu1x.dev>
Date: Mon, 21 Jun 2021 21:43:21 +0200
Subject: [PATCH 01/13] Implement Yaroslavskiy-Bentley-Bloch Quicksort.

---
 src/sort.rs   | 234 +++++++++++++++++++++++++++++---------------------
 tests/sort.rs |  23 +++--
 2 files changed, 149 insertions(+), 108 deletions(-)

diff --git a/src/sort.rs b/src/sort.rs
index f43a95b1..0e4631e4 100644
--- a/src/sort.rs
+++ b/src/sort.rs
@@ -1,8 +1,6 @@
 use indexmap::IndexMap;
 use ndarray::prelude::*;
 use ndarray::{Data, DataMut, Slice};
-use rand::prelude::*;
-use rand::thread_rng;
 
 /// Methods for sorting and partitioning 1-D arrays.
 pub trait Sort1dExt<A, S>
@@ -50,26 +48,21 @@ where
         S: DataMut,
         S2: Data<Elem = usize>;
 
-    /// Partitions the array in increasing order based on the value initially
-    /// located at `pivot_index` and returns the new index of the value.
+    /// Partitions the array in increasing order based on the values initially located at `0` and
+    /// `n - 1` where `n` is the number of elements in the array and returns the new indexes of the
+    /// values.
     ///
-    /// The elements are rearranged in such a way that the value initially
-    /// located at `pivot_index` is moved to the position it would be in an
-    /// array sorted in increasing order. The return value is the new index of
-    /// the value after rearrangement. All elements smaller than the value are
-    /// moved to its left and all elements equal or greater than the value are
-    /// moved to its right. The ordering of the elements in the two partitions
-    /// is undefined.
+    /// The elements are rearranged in such a way that the values initially located at `0` and
+    /// `n - 1` are moved to the position it would be in an array sorted in increasing order. The
+    /// return values are the new indexes of the values after rearrangement. All elements less than
+    /// the values are moved to their left and all elements equal or greater than the values are
+    /// moved to their right. The ordering of the elements in the three partitions is undefined.
     ///
-    /// `self` is shuffled **in place** to operate the desired partition:
-    /// no copy of the array is allocated.
+    /// The array is shuffled **in place**, no copy of the array is allocated.
     ///
-    /// The method uses Hoare's partition algorithm.
-    /// Complexity: O(`n`), where `n` is the number of elements in the array.
-    /// Average number of element swaps: n/6 - 1/3 (see
-    /// [link](https://cs.stackexchange.com/questions/11458/quicksort-partitioning-hoare-vs-lomuto/11550))
+    /// This method implements the partitioning scheme of [Yaroslavskiy-Bentley-Bloch Quicksort].
     ///
-    /// **Panics** if `pivot_index` is greater than or equal to `n`.
+    /// [Yaroslavskiy-Bentley-Bloch Quicksort]: https://api.semanticscholar.org/CorpusID:51871084
     ///
     /// # Example
     ///
@@ -78,23 +71,30 @@ where
     /// use ndarray_stats::Sort1dExt;
     ///
     /// let mut data = array![3, 1, 4, 5, 2];
-    /// let pivot_index = 2;
-    /// let pivot_value = data[pivot_index];
+    /// // Sorted pivot values.
+    /// let (lower_value, upper_value) = (data[data.len() - 1], data[0]);
     ///
-    /// // Partition by the value located at `pivot_index`.
-    /// let new_index = data.partition_mut(pivot_index);
-    /// // The pivot value is now located at `new_index`.
-    /// assert_eq!(data[new_index], pivot_value);
-    /// // Elements less than that value are moved to the left.
-    /// for i in 0..new_index {
-    ///     assert!(data[i] < pivot_value);
+    /// // Partitions by the values located at `0` and `data.len() - 1`.
+    /// let (lower_index, upper_index) = data.partition_mut();
+    /// // The pivot values are now located at `lower_index` and `upper_index`.
+    /// assert_eq!(data[lower_index], lower_value);
+    /// assert_eq!(data[upper_index], upper_value);
+    /// // Elements lower than the lower pivot value are moved to its left.
+    /// for i in 0..lower_index {
+    ///     assert!(data[i] < lower_value);
+    /// }
+    /// // Elements greater than or equal the lower pivot value and less than or equal the upper
+    /// // pivot value are moved between the two pivot indexes.
+    /// for i in lower_index + 1..upper_index {
+    ///     assert!(lower_value <= data[i]);
+    ///     assert!(data[i] <= upper_value);
     /// }
-    /// // Elements greater than or equal to that value are moved to the right.
-    /// for i in (new_index + 1)..data.len() {
-    ///      assert!(data[i] >= pivot_value);
+    /// // Elements greater than or equal the upper pivot value are moved to its right.
+    /// for i in upper_index + 1..data.len() {
+    ///     assert!(upper_value <= data[i]);
     /// }
     /// ```
-    fn partition_mut(&mut self, pivot_index: usize) -> usize
+    fn partition_mut(&mut self) -> (usize, usize)
     where
         A: Ord + Clone,
         S: DataMut;
@@ -115,17 +115,20 @@ where
         if n == 1 {
             self[0].clone()
         } else {
-            let mut rng = thread_rng();
-            let pivot_index = rng.gen_range(0..n);
-            let partition_index = self.partition_mut(pivot_index);
-            if i < partition_index {
-                self.slice_axis_mut(Axis(0), Slice::from(..partition_index))
+            let (lower_index, upper_index) = self.partition_mut();
+            if i < lower_index {
+                self.slice_axis_mut(Axis(0), Slice::from(..lower_index))
                     .get_from_sorted_mut(i)
-            } else if i == partition_index {
+            } else if i == lower_index {
+                self[i].clone()
+            } else if i < upper_index {
+                self.slice_axis_mut(Axis(0), Slice::from(lower_index + 1..upper_index))
+                    .get_from_sorted_mut(i - (lower_index + 1))
+            } else if i == upper_index {
                 self[i].clone()
             } else {
-                self.slice_axis_mut(Axis(0), Slice::from(partition_index + 1..))
-                    .get_from_sorted_mut(i - (partition_index + 1))
+                self.slice_axis_mut(Axis(0), Slice::from(upper_index + 1..))
+                    .get_from_sorted_mut(i - (upper_index + 1))
             }
         }
     }
@@ -143,42 +146,51 @@ where
         get_many_from_sorted_mut_unchecked(self, &deduped_indexes)
     }
 
-    fn partition_mut(&mut self, pivot_index: usize) -> usize
+    fn partition_mut(&mut self) -> (usize, usize)
     where
         A: Ord + Clone,
         S: DataMut,
     {
-        let pivot_value = self[pivot_index].clone();
-        self.swap(pivot_index, 0);
-        let n = self.len();
-        let mut i = 1;
-        let mut j = n - 1;
-        loop {
-            loop {
-                if i > j {
-                    break;
-                }
-                if self[i] >= pivot_value {
-                    break;
+        // Sort `lowermost` and `uppermost` elements and use them as dual pivot.
+        let lowermost = 0;
+        let uppermost = self.len() - 1;
+        if self[lowermost] > self[uppermost] {
+            self.swap(lowermost, uppermost);
+        }
+        // Increasing running and partition index starting after lower pivot.
+        let mut index = lowermost + 1;
+        let mut lower = lowermost + 1;
+        // Decreasing partition index starting before upper pivot.
+        let mut upper = uppermost - 1;
+        // Swap elements at `index` into their partitions.
+        while index <= upper {
+            if self[index] < self[lowermost] {
+                // Swap elements into lower partition.
+                self.swap(index, lower);
+                lower += 1;
+            } else if self[index] >= self[uppermost] {
+                // Search first element of upper partition.
+                while self[upper] > self[uppermost] && index < upper {
+                    upper -= 1;
                 }
-                i += 1;
-            }
-            while pivot_value <= self[j] {
-                if j == 1 {
-                    break;
+                // Swap elements into upper partition.
+                self.swap(index, upper);
+                if self[index] < self[lowermost] {
+                    // Swap swapped elements into lower partition.
+                    self.swap(index, lower);
+                    lower += 1;
                 }
-                j -= 1;
-            }
-            if i >= j {
-                break;
-            } else {
-                self.swap(i, j);
-                i += 1;
-                j -= 1;
+                upper -= 1;
             }
+            index += 1;
         }
-        self.swap(0, i - 1);
-        i - 1
+        lower -= 1;
+        upper += 1;
+        // Swap pivots to their new indexes.
+        self.swap(lowermost, lower);
+        self.swap(uppermost, upper);
+        // Lower and upper pivot index.
+        (lower, upper)
     }
 
     private_impl! {}
@@ -249,50 +261,72 @@ fn _get_many_from_sorted_mut_unchecked<A>(
         return;
     }
 
-    // We pick a random pivot index: the corresponding element is the pivot value
-    let mut rng = thread_rng();
-    let pivot_index = rng.gen_range(0..n);
+    // We partition the array with respect to the two pivot values. The pivot values move to
+    // `lower_index` and `upper_index`.
+    //
+    // Elements strictly less than the lower pivot value have indexes < `lower_index`.
+    //
+    // Elements greater than or equal the lower pivot value and less than or equal the upper pivot
+    // value have indexes > `lower_index` and < `upper_index`.
+    //
+    // Elements less than or equal the upper pivot value have indexes > `upper_index`.
+    let (lower_index, upper_index) = array.partition_mut();
 
-    // We partition the array with respect to the pivot value.
-    // The pivot value moves to `array_partition_index`.
-    // Elements strictly smaller than the pivot value have indexes < `array_partition_index`.
-    // Elements greater or equal to the pivot value have indexes > `array_partition_index`.
-    let array_partition_index = array.partition_mut(pivot_index);
+    // We use a divide-and-conquer strategy, splitting the indexes we are searching for (`indexes`)
+    // and the corresponding portions of the output slice (`values`) into partitions with respect to
+    // `lower_index` and `upper_index`.
+    let (found_exact, split_index) = match indexes.binary_search(&lower_index) {
+        Ok(index) => (true, index),
+        Err(index) => (false, index),
+    };
+    let (lower_indexes, inner_indexes) = indexes.split_at_mut(split_index);
+    let (lower_values, inner_values) = values.split_at_mut(split_index);
+    let (upper_indexes, upper_values) = if found_exact {
+        inner_values[0] = array[lower_index].clone(); // Write exactly found value.
+        (&mut inner_indexes[1..], &mut inner_values[1..])
+    } else {
+        (inner_indexes, inner_values)
+    };
 
-    // We use a divide-and-conquer strategy, splitting the indexes we are
-    // searching for (`indexes`) and the corresponding portions of the output
-    // slice (`values`) into pieces with respect to `array_partition_index`.
-    let (found_exact, index_split) = match indexes.binary_search(&array_partition_index) {
+    let (found_exact, split_index) = match upper_indexes.binary_search(&upper_index) {
         Ok(index) => (true, index),
         Err(index) => (false, index),
     };
-    let (smaller_indexes, other_indexes) = indexes.split_at_mut(index_split);
-    let (smaller_values, other_values) = values.split_at_mut(index_split);
-    let (bigger_indexes, bigger_values) = if found_exact {
-        other_values[0] = array[array_partition_index].clone(); // Write exactly found value.
-        (&mut other_indexes[1..], &mut other_values[1..])
+    let (inner_indexes, upper_indexes) = upper_indexes.split_at_mut(split_index);
+    let (inner_values, upper_values) = upper_values.split_at_mut(split_index);
+    let (upper_indexes, upper_values) = if found_exact {
+        upper_values[0] = array[upper_index].clone(); // Write exactly found value.
+        (&mut upper_indexes[1..], &mut upper_values[1..])
     } else {
-        (other_indexes, other_values)
+        (upper_indexes, upper_values)
     };
 
-    // We search recursively for the values corresponding to strictly smaller
-    // indexes to the left of `partition_index`.
+    // We search recursively for the values corresponding to indexes strictly less than
+    // `lower_index` in the lower partition.
+    _get_many_from_sorted_mut_unchecked(
+        array.slice_axis_mut(Axis(0), Slice::from(..lower_index)),
+        lower_indexes,
+        lower_values,
+    );
+
+    // We search recursively for the values corresponding to indexes greater than or equal
+    // `lower_index` in the inner partition, that is between the lower and upper partition. Since
+    // only the inner partition of the array is passed in, the indexes need to be shifted by length
+    // of the lower partition.
+    inner_indexes.iter_mut().for_each(|x| *x -= lower_index + 1);
     _get_many_from_sorted_mut_unchecked(
-        array.slice_axis_mut(Axis(0), Slice::from(..array_partition_index)),
-        smaller_indexes,
-        smaller_values,
+        array.slice_axis_mut(Axis(0), Slice::from(lower_index + 1..upper_index)),
+        inner_indexes,
+        inner_values,
     );
 
-    // We search recursively for the values corresponding to strictly bigger
-    // indexes to the right of `partition_index`. Since only the right portion
-    // of the array is passed in, the indexes need to be shifted by length of
-    // the removed portion.
-    bigger_indexes
-        .iter_mut()
-        .for_each(|x| *x -= array_partition_index + 1);
+    // We search recursively for the values corresponding to indexes greater than or equal
+    // `upper_index` in the upper partition. Since only the upper partition of the array is passed
+    // in, the indexes need to be shifted by the combined length of the lower and inner partition.
+    upper_indexes.iter_mut().for_each(|x| *x -= upper_index + 1);
     _get_many_from_sorted_mut_unchecked(
-        array.slice_axis_mut(Axis(0), Slice::from(array_partition_index + 1..)),
-        bigger_indexes,
-        bigger_values,
+        array.slice_axis_mut(Axis(0), Slice::from(upper_index + 1..)),
+        upper_indexes,
+        upper_values,
     );
 }
diff --git a/tests/sort.rs b/tests/sort.rs
index b2bd12f1..ff644356 100644
--- a/tests/sort.rs
+++ b/tests/sort.rs
@@ -19,15 +19,22 @@ fn test_partition_mut() {
     ];
     for a in l.iter_mut() {
         let n = a.len();
-        let pivot_index = n - 1;
-        let pivot_value = a[pivot_index].clone();
-        let partition_index = a.partition_mut(pivot_index);
-        for i in 0..partition_index {
-            assert!(a[i] < pivot_value);
+        let (mut lower_value, mut upper_value) = (a[0].clone(), a[n - 1].clone());
+        if lower_value > upper_value {
+            std::mem::swap(&mut lower_value, &mut upper_value);
         }
-        assert_eq!(a[partition_index], pivot_value);
-        for j in (partition_index + 1)..n {
-            assert!(pivot_value <= a[j]);
+        let (lower_index, upper_index) = a.partition_mut();
+        for i in 0..lower_index {
+            assert!(a[i] < lower_value);
+        }
+        assert_eq!(a[lower_index], lower_value);
+        for i in lower_index + 1..upper_index {
+            assert!(lower_value <= a[i]);
+            assert!(a[i] <= upper_value);
+        }
+        assert_eq!(a[upper_index], upper_value);
+        for i in (upper_index + 1)..n {
+            assert!(upper_value <= a[i]);
         }
     }
 }

From e6b7b02ba9e200f3bdfdf88daabc7421fcb3eee7 Mon Sep 17 00:00:00 2001
From: Rouven Spreckels <rs@qu1x.dev>
Date: Tue, 22 Jun 2021 22:29:26 +0200
Subject: [PATCH 02/13] Implement skewed pivot sampling.

---
 src/sort.rs | 48 +++++++++++++++++++++++++++++++++++-------------
 1 file changed, 35 insertions(+), 13 deletions(-)

diff --git a/src/sort.rs b/src/sort.rs
index 0e4631e4..3e19f2ea 100644
--- a/src/sort.rs
+++ b/src/sort.rs
@@ -48,21 +48,21 @@ where
         S: DataMut,
         S2: Data<Elem = usize>;
 
-    /// Partitions the array in increasing order based on the values initially located at `0` and
-    /// `n - 1` where `n` is the number of elements in the array and returns the new indexes of the
-    /// values.
+    /// Partitions the array in increasing order at two skewed pivot values as 1st and 3rd element
+    /// of a sorted sample of 5 equally spaced elements around the center and returns their indexes.
+    /// For arrays of less than 42 elements the outermost elements serve as sample for pivot values.
     ///
-    /// The elements are rearranged in such a way that the values initially located at `0` and
-    /// `n - 1` are moved to the position it would be in an array sorted in increasing order. The
-    /// return values are the new indexes of the values after rearrangement. All elements less than
-    /// the values are moved to their left and all elements equal or greater than the values are
-    /// moved to their right. The ordering of the elements in the three partitions is undefined.
+    /// The elements are rearranged in such a way that the two pivot values are moved to the indexes
+    /// they would be in an array sorted in increasing order. The return values are the new indexes
+    /// of the values after rearrangement. All elements less than the values are moved to their left
+    /// and all elements equal or greater than the values are moved to their right. The ordering of
+    /// the elements in the three partitions is undefined.
     ///
     /// The array is shuffled **in place**, no copy of the array is allocated.
     ///
-    /// This method implements the partitioning scheme of [Yaroslavskiy-Bentley-Bloch Quicksort].
+    /// This method performs [dual-pivot partitioning] with skewed pivot sampling.
     ///
-    /// [Yaroslavskiy-Bentley-Bloch Quicksort]: https://api.semanticscholar.org/CorpusID:51871084
+    /// [dual-pivot partitioning]: https://www.wild-inter.net/publications/wild-2018b.pdf
     ///
     /// # Example
     ///
@@ -151,11 +151,33 @@ where
         A: Ord + Clone,
         S: DataMut,
     {
-        // Sort `lowermost` and `uppermost` elements and use them as dual pivot.
         let lowermost = 0;
         let uppermost = self.len() - 1;
-        if self[lowermost] > self[uppermost] {
-            self.swap(lowermost, uppermost);
+        if self.len() < 42 {
+            // Sort outermost elements and use them as pivots.
+            if self[lowermost] > self[uppermost] {
+                self.swap(lowermost, uppermost);
+            }
+        } else {
+            // Sample indexes of 5 evenly spaced elements around the center element.
+            let mut samples = [0; 5];
+            // Assume array of at least 7 elements.
+            let seventh = self.len() / (samples.len() + 2);
+            samples[2] = self.len() / 2;
+            samples[1] = samples[2] - seventh;
+            samples[0] = samples[1] - seventh;
+            samples[3] = samples[2] + seventh;
+            samples[4] = samples[3] + seventh;
+            // Use insertion sort for sample elements by looking up their indexes.
+            for mut index in 1..samples.len() {
+                while index > 0 && self[samples[index - 1]] > self[samples[index]] {
+                    self.swap(samples[index - 1], samples[index]);
+                    index -= 1;
+                }
+            }
+            // Use 1st and 3rd element of sorted sample as skewed pivots.
+            self.swap(lowermost, samples[0]);
+            self.swap(uppermost, samples[2]);
         }
         // Increasing running and partition index starting after lower pivot.
         let mut index = lowermost + 1;

From 0919c2bc57bf41124134f9229d9be300c678c09a Mon Sep 17 00:00:00 2001
From: Rouven Spreckels <rs@qu1x.dev>
Date: Mon, 28 Jun 2021 21:12:36 +0200
Subject: [PATCH 03/13] Adapt pivot sampling to relative sought rank.

---
 src/sort.rs   | 276 +++++++++++++++++++++++++++++++++++++++-----------
 tests/sort.rs |  44 ++++++--
 2 files changed, 249 insertions(+), 71 deletions(-)

diff --git a/src/sort.rs b/src/sort.rs
index 3e19f2ea..b695b691 100644
--- a/src/sort.rs
+++ b/src/sort.rs
@@ -20,10 +20,10 @@ where
     /// No other assumptions should be made on the ordering of the
     /// elements after this computation.
     ///
-    /// Complexity ([quickselect](https://en.wikipedia.org/wiki/Quickselect)):
-    /// - average case: O(`n`);
-    /// - worst case: O(`n`^2);
-    /// where n is the number of elements in the array.
+    /// This method performs a variant of [Sesquickselect] with pivot samples of 5 equally spaced
+    /// elements around the center.
+    ///
+    /// [Sesquickselect]: https://www.wild-inter.net/publications/martinez-nebel-wild-2019.pdf
     ///
     /// **Panics** if `i` is greater than or equal to `n`.
     fn get_from_sorted_mut(&mut self, i: usize) -> A
@@ -48,9 +48,65 @@ where
         S: DataMut,
         S2: Data<Elem = usize>;
 
-    /// Partitions the array in increasing order at two skewed pivot values as 1st and 3rd element
-    /// of a sorted sample of 5 equally spaced elements around the center and returns their indexes.
-    /// For arrays of less than 42 elements the outermost elements serve as sample for pivot values.
+    /// Sorts a sample of 5 equally spaced elements around the center and returns their indexes.
+    ///
+    /// **Panics** for arrays of less than 7 elements.
+    fn sample_mut(&mut self) -> [usize; 5]
+    where
+        A: Ord + Clone,
+        S: DataMut;
+
+    /// Partitions the array in increasing order based on the value initially
+    /// located at `pivot_index` and returns the new index of the value.
+    ///
+    /// The elements are rearranged in such a way that the value initially
+    /// located at `pivot_index` is moved to the position it would be in an
+    /// array sorted in increasing order. The return value is the new index of
+    /// the value after rearrangement. All elements smaller than the value are
+    /// moved to its left and all elements equal or greater than the value are
+    /// moved to its right. The ordering of the elements in the two partitions
+    /// is undefined.
+    ///
+    /// `self` is shuffled **in place** to operate the desired partition:
+    /// no copy of the array is allocated.
+    ///
+    /// The method uses Hoare's partition algorithm.
+    /// Complexity: O(`n`), where `n` is the number of elements in the array.
+    /// Average number of element swaps: n/6 - 1/3 (see
+    /// [link](https://cs.stackexchange.com/questions/11458/quicksort-partitioning-hoare-vs-lomuto/11550))
+    ///
+    /// **Panics** if `pivot_index` is greater than or equal to `n`.
+    ///
+    /// # Example
+    ///
+    /// ```
+    /// use ndarray::array;
+    /// use ndarray_stats::Sort1dExt;
+    ///
+    /// let mut data = array![3, 1, 4, 5, 2];
+    /// let pivot_index = 2;
+    /// let pivot_value = data[pivot_index];
+    ///
+    /// // Partition by the value located at `pivot_index`.
+    /// let new_index = data.partition_mut(pivot_index);
+    /// // The pivot value is now located at `new_index`.
+    /// assert_eq!(data[new_index], pivot_value);
+    /// // Elements less than that value are moved to the left.
+    /// for i in 0..new_index {
+    ///     assert!(data[i] < pivot_value);
+    /// }
+    /// // Elements greater than or equal to that value are moved to the right.
+    /// for i in (new_index + 1)..data.len() {
+    ///      assert!(data[i] >= pivot_value);
+    /// }
+    /// ```
+    fn partition_mut(&mut self, pivot_index: usize) -> usize
+    where
+        A: Ord + Clone,
+        S: DataMut;
+
+    /// Partitions the array in increasing order based on the values initially located at the two
+    /// pivot indexes `lower` and `upper` and returns the new indexes of their values.
     ///
     /// The elements are rearranged in such a way that the two pivot values are moved to the indexes
     /// they would be in an array sorted in increasing order. The return values are the new indexes
@@ -60,10 +116,12 @@ where
     ///
     /// The array is shuffled **in place**, no copy of the array is allocated.
     ///
-    /// This method performs [dual-pivot partitioning] with skewed pivot sampling.
+    /// This method performs [dual-pivot partitioning].
     ///
     /// [dual-pivot partitioning]: https://www.wild-inter.net/publications/wild-2018b.pdf
     ///
+    /// **Panics** if `lower` or `upper` is out of bound.
+    ///
     /// # Example
     ///
     /// ```
@@ -71,11 +129,11 @@ where
     /// use ndarray_stats::Sort1dExt;
     ///
     /// let mut data = array![3, 1, 4, 5, 2];
-    /// // Sorted pivot values.
-    /// let (lower_value, upper_value) = (data[data.len() - 1], data[0]);
+    /// // Skewed pivot values.
+    /// let (lower_value, upper_value) = (1, 5);
     ///
-    /// // Partitions by the values located at `0` and `data.len() - 1`.
-    /// let (lower_index, upper_index) = data.partition_mut();
+    /// // Partitions by the values located at `1` and `3`.
+    /// let (lower_index, upper_index) = data.dual_partition_mut(1, 3);
     /// // The pivot values are now located at `lower_index` and `upper_index`.
     /// assert_eq!(data[lower_index], lower_value);
     /// assert_eq!(data[upper_index], upper_value);
@@ -94,7 +152,7 @@ where
     ///     assert!(upper_value <= data[i]);
     /// }
     /// ```
-    fn partition_mut(&mut self) -> (usize, usize)
+    fn dual_partition_mut(&mut self, lower: usize, upper: usize) -> (usize, usize)
     where
         A: Ord + Clone,
         S: DataMut;
@@ -112,23 +170,67 @@ where
         S: DataMut,
     {
         let n = self.len();
-        if n == 1 {
-            self[0].clone()
+        // Recursion cutoff at integer multiple of sample space divider of 7 elements.
+        if n < 21 {
+            for mut index in 1..n {
+                while index > 0 && self[index - 1] > self[index] {
+                    self.swap(index - 1, index);
+                    index -= 1;
+                }
+            }
+            self[i].clone()
         } else {
-            let (lower_index, upper_index) = self.partition_mut();
-            if i < lower_index {
-                self.slice_axis_mut(Axis(0), Slice::from(..lower_index))
-                    .get_from_sorted_mut(i)
-            } else if i == lower_index {
-                self[i].clone()
-            } else if i < upper_index {
-                self.slice_axis_mut(Axis(0), Slice::from(lower_index + 1..upper_index))
-                    .get_from_sorted_mut(i - (lower_index + 1))
-            } else if i == upper_index {
-                self[i].clone()
+            // Adapt pivot sampling to relative sought rank and switch from dual-pivot to
+            // single-pivot partitioning for extreme sought ranks.
+            let sought_rank = i as f64 / n as f64;
+            if (0.036..=0.964).contains(&sought_rank) {
+                let (lower_index, upper_index) = if sought_rank <= 0.5 {
+                    if sought_rank <= 0.153 {
+                        (0, 1) // (0, 0, 3)
+                    } else {
+                        (0, 2) // (0, 1, 2)
+                    }
+                } else {
+                    if sought_rank <= 0.847 {
+                        (2, 4) // (2, 1, 0)
+                    } else {
+                        (3, 4) // (3, 0, 0)
+                    }
+                };
+                let sample = self.sample_mut();
+                let (lower_index, upper_index) =
+                    self.dual_partition_mut(sample[lower_index], sample[upper_index]);
+                if i < lower_index {
+                    self.slice_axis_mut(Axis(0), Slice::from(..lower_index))
+                        .get_from_sorted_mut(i)
+                } else if i == lower_index {
+                    self[i].clone()
+                } else if i < upper_index {
+                    self.slice_axis_mut(Axis(0), Slice::from(lower_index + 1..upper_index))
+                        .get_from_sorted_mut(i - (lower_index + 1))
+                } else if i == upper_index {
+                    self[i].clone()
+                } else {
+                    self.slice_axis_mut(Axis(0), Slice::from(upper_index + 1..))
+                        .get_from_sorted_mut(i - (upper_index + 1))
+                }
             } else {
-                self.slice_axis_mut(Axis(0), Slice::from(upper_index + 1..))
-                    .get_from_sorted_mut(i - (upper_index + 1))
+                let sample = self.sample_mut();
+                let pivot_index = if sought_rank <= 0.5 {
+                    0 // (0, 4)
+                } else {
+                    4 // (4, 0)
+                };
+                let pivot_index = self.partition_mut(sample[pivot_index]);
+                if i < pivot_index {
+                    self.slice_axis_mut(Axis(0), Slice::from(..pivot_index))
+                        .get_from_sorted_mut(i)
+                } else if i == pivot_index {
+                    self[i].clone()
+                } else {
+                    self.slice_axis_mut(Axis(0), Slice::from(pivot_index + 1..))
+                        .get_from_sorted_mut(i - (pivot_index + 1))
+                }
             }
         }
     }
@@ -146,38 +248,80 @@ where
         get_many_from_sorted_mut_unchecked(self, &deduped_indexes)
     }
 
-    fn partition_mut(&mut self) -> (usize, usize)
+    fn sample_mut(&mut self) -> [usize; 5]
     where
         A: Ord + Clone,
         S: DataMut,
     {
-        let lowermost = 0;
-        let uppermost = self.len() - 1;
-        if self.len() < 42 {
-            // Sort outermost elements and use them as pivots.
-            if self[lowermost] > self[uppermost] {
-                self.swap(lowermost, uppermost);
+        // Space between sample indexes.
+        let space = self.len() / 7;
+        assert!(space > 0, "Cannot sample array of less then 7 elements");
+        // Initialize sample indexes with their lowermost index.
+        let mut samples = [self.len() / 2 - 2 * space; 5];
+        // Equally space sample indexes and sort by their values by looking up their indexes.
+        for mut index in 1..samples.len() {
+            // Equally space sample indexes based on their lowermost index.
+            samples[index] += index * space;
+            // Insertion sort looking up only the already equally spaced sample indexes.
+            while index > 0 && self[samples[index - 1]] > self[samples[index]] {
+                self.swap(samples[index - 1], samples[index]);
+                index -= 1;
             }
-        } else {
-            // Sample indexes of 5 evenly spaced elements around the center element.
-            let mut samples = [0; 5];
-            // Assume array of at least 7 elements.
-            let seventh = self.len() / (samples.len() + 2);
-            samples[2] = self.len() / 2;
-            samples[1] = samples[2] - seventh;
-            samples[0] = samples[1] - seventh;
-            samples[3] = samples[2] + seventh;
-            samples[4] = samples[3] + seventh;
-            // Use insertion sort for sample elements by looking up their indexes.
-            for mut index in 1..samples.len() {
-                while index > 0 && self[samples[index - 1]] > self[samples[index]] {
-                    self.swap(samples[index - 1], samples[index]);
-                    index -= 1;
+        }
+        samples
+    }
+
+    fn partition_mut(&mut self, pivot_index: usize) -> usize
+    where
+        A: Ord + Clone,
+        S: DataMut,
+    {
+        let pivot_value = self[pivot_index].clone();
+        self.swap(pivot_index, 0);
+        let n = self.len();
+        let mut i = 1;
+        let mut j = n - 1;
+        loop {
+            loop {
+                if i > j {
+                    break;
                 }
+                if self[i] >= pivot_value {
+                    break;
+                }
+                i += 1;
+            }
+            while pivot_value <= self[j] {
+                if j == 1 {
+                    break;
+                }
+                j -= 1;
+            }
+            if i >= j {
+                break;
+            } else {
+                self.swap(i, j);
+                i += 1;
+                j -= 1;
             }
-            // Use 1st and 3rd element of sorted sample as skewed pivots.
-            self.swap(lowermost, samples[0]);
-            self.swap(uppermost, samples[2]);
+        }
+        self.swap(0, i - 1);
+        i - 1
+    }
+
+    fn dual_partition_mut(&mut self, lower: usize, upper: usize) -> (usize, usize)
+    where
+        A: Ord + Clone,
+        S: DataMut,
+    {
+        let lowermost = 0;
+        let uppermost = self.len() - 1;
+        // Swap pivots with outermost elements.
+        self.swap(lowermost, lower);
+        self.swap(uppermost, upper);
+        if self[lowermost] > self[uppermost] {
+            // Sort pivots instead of panicking via assertion.
+            self.swap(lowermost, uppermost);
         }
         // Increasing running and partition index starting after lower pivot.
         let mut index = lowermost + 1;
@@ -274,16 +418,26 @@ fn _get_many_from_sorted_mut_unchecked<A>(
         return;
     }
 
-    // At this point, `n >= 1` since `indexes.len() >= 1`.
-    if n == 1 {
-        // We can only reach this point if `indexes.len() == 1`, so we only
-        // need to assign the single value, and then we're done.
-        debug_assert_eq!(indexes.len(), 1);
-        values[0] = array[0].clone();
+    // Recursion cutoff at integer multiple of sample space divider of 7 elements.
+    if n < 21 {
+        for mut index in 1..n {
+            while index > 0 && array[index - 1] > array[index] {
+                array.swap(index - 1, index);
+                index -= 1;
+            }
+        }
+        for (value, index) in values.iter_mut().zip(indexes.iter()) {
+            *value = array[*index].clone();
+        }
         return;
     }
 
-    // We partition the array with respect to the two pivot values. The pivot values move to
+    // Since there is no single sought rank to adapt pivot sampling to, the recommended skewed pivot
+    // sampling of dual-pivot Quicksort is used.
+    let sample = array.sample_mut();
+    let (lower_index, upper_index) = (sample[0], sample[2]); // (0, 1, 2)
+
+    // We partition the array with respect to the two pivot values. The pivot values move to the new
     // `lower_index` and `upper_index`.
     //
     // Elements strictly less than the lower pivot value have indexes < `lower_index`.
@@ -292,7 +446,7 @@ fn _get_many_from_sorted_mut_unchecked<A>(
     // value have indexes > `lower_index` and < `upper_index`.
     //
     // Elements less than or equal the upper pivot value have indexes > `upper_index`.
-    let (lower_index, upper_index) = array.partition_mut();
+    let (lower_index, upper_index) = array.dual_partition_mut(lower_index, upper_index);
 
     // We use a divide-and-conquer strategy, splitting the indexes we are searching for (`indexes`)
     // and the corresponding portions of the output slice (`values`) into partitions with respect to
diff --git a/tests/sort.rs b/tests/sort.rs
index ff644356..b298e88d 100644
--- a/tests/sort.rs
+++ b/tests/sort.rs
@@ -19,22 +19,46 @@ fn test_partition_mut() {
     ];
     for a in l.iter_mut() {
         let n = a.len();
-        let (mut lower_value, mut upper_value) = (a[0].clone(), a[n - 1].clone());
-        if lower_value > upper_value {
-            std::mem::swap(&mut lower_value, &mut upper_value);
+        let pivot_index = n - 1;
+        let pivot_value = a[pivot_index].clone();
+        let partition_index = a.partition_mut(pivot_index);
+        for i in 0..partition_index {
+            assert!(a[i] < pivot_value);
         }
-        let (lower_index, upper_index) = a.partition_mut();
+        assert_eq!(a[partition_index], pivot_value);
+        for j in (partition_index + 1)..n {
+            assert!(pivot_value <= a[j]);
+        }
+    }
+}
+
+#[test]
+fn test_dual_partition_mut() {
+    let mut l = vec![
+        arr1(&[1, 1, 1, 1, 1]),
+        arr1(&[1, 3, 2, 10, 10]),
+        arr1(&[2, 3, 4, 1]),
+        arr1(&[
+            355, 453, 452, 391, 289, 343, 44, 154, 271, 44, 314, 276, 160, 469, 191, 138, 163, 308,
+            395, 3, 416, 391, 210, 354, 200,
+        ]),
+        arr1(&[
+            84, 192, 216, 159, 89, 296, 35, 213, 456, 278, 98, 52, 308, 418, 329, 173, 286, 106,
+            366, 129, 125, 450, 23, 463, 151,
+        ]),
+    ];
+    for a in l.iter_mut() {
+        let n = a.len();
+        let (lower_index, upper_index) = a.dual_partition_mut(1, a.len() / 2);
         for i in 0..lower_index {
-            assert!(a[i] < lower_value);
+            assert!(a[i] < a[lower_index]);
         }
-        assert_eq!(a[lower_index], lower_value);
         for i in lower_index + 1..upper_index {
-            assert!(lower_value <= a[i]);
-            assert!(a[i] <= upper_value);
+            assert!(a[lower_index] <= a[i]);
+            assert!(a[i] <= a[upper_index]);
         }
-        assert_eq!(a[upper_index], upper_value);
         for i in (upper_index + 1)..n {
-            assert!(upper_value <= a[i]);
+            assert!(a[upper_index] <= a[i]);
         }
     }
 }

From fc13897002acf85a3b692fd88967333046cfb1eb Mon Sep 17 00:00:00 2001
From: Rouven Spreckels <rs@qu1x.dev>
Date: Wed, 30 Jun 2021 16:03:05 +0200
Subject: [PATCH 04/13] Make sampling an implementation detail.

---
 src/sort.rs | 67 +++++++++++++++++++++++++----------------------------
 1 file changed, 31 insertions(+), 36 deletions(-)

diff --git a/src/sort.rs b/src/sort.rs
index b695b691..54f9d4a5 100644
--- a/src/sort.rs
+++ b/src/sort.rs
@@ -20,8 +20,7 @@ where
     /// No other assumptions should be made on the ordering of the
     /// elements after this computation.
     ///
-    /// This method performs a variant of [Sesquickselect] with pivot samples of 5 equally spaced
-    /// elements around the center.
+    /// This method performs [Sesquickselect].
     ///
     /// [Sesquickselect]: https://www.wild-inter.net/publications/martinez-nebel-wild-2019.pdf
     ///
@@ -48,14 +47,6 @@ where
         S: DataMut,
         S2: Data<Elem = usize>;
 
-    /// Sorts a sample of 5 equally spaced elements around the center and returns their indexes.
-    ///
-    /// **Panics** for arrays of less than 7 elements.
-    fn sample_mut(&mut self) -> [usize; 5]
-    where
-        A: Ord + Clone,
-        S: DataMut;
-
     /// Partitions the array in increasing order based on the value initially
     /// located at `pivot_index` and returns the new index of the value.
     ///
@@ -182,6 +173,8 @@ where
         } else {
             // Adapt pivot sampling to relative sought rank and switch from dual-pivot to
             // single-pivot partitioning for extreme sought ranks.
+            let mut sample = [0; 5];
+            sample_mut(self, &mut sample);
             let sought_rank = i as f64 / n as f64;
             if (0.036..=0.964).contains(&sought_rank) {
                 let (lower_index, upper_index) = if sought_rank <= 0.5 {
@@ -197,7 +190,6 @@ where
                         (3, 4) // (3, 0, 0)
                     }
                 };
-                let sample = self.sample_mut();
                 let (lower_index, upper_index) =
                     self.dual_partition_mut(sample[lower_index], sample[upper_index]);
                 if i < lower_index {
@@ -215,7 +207,6 @@ where
                         .get_from_sorted_mut(i - (upper_index + 1))
                 }
             } else {
-                let sample = self.sample_mut();
                 let pivot_index = if sought_rank <= 0.5 {
                     0 // (0, 4)
                 } else {
@@ -248,29 +239,6 @@ where
         get_many_from_sorted_mut_unchecked(self, &deduped_indexes)
     }
 
-    fn sample_mut(&mut self) -> [usize; 5]
-    where
-        A: Ord + Clone,
-        S: DataMut,
-    {
-        // Space between sample indexes.
-        let space = self.len() / 7;
-        assert!(space > 0, "Cannot sample array of less then 7 elements");
-        // Initialize sample indexes with their lowermost index.
-        let mut samples = [self.len() / 2 - 2 * space; 5];
-        // Equally space sample indexes and sort by their values by looking up their indexes.
-        for mut index in 1..samples.len() {
-            // Equally space sample indexes based on their lowermost index.
-            samples[index] += index * space;
-            // Insertion sort looking up only the already equally spaced sample indexes.
-            while index > 0 && self[samples[index - 1]] > self[samples[index]] {
-                self.swap(samples[index - 1], samples[index]);
-                index -= 1;
-            }
-        }
-        samples
-    }
-
     fn partition_mut(&mut self, pivot_index: usize) -> usize
     where
         A: Ord + Clone,
@@ -434,7 +402,8 @@ fn _get_many_from_sorted_mut_unchecked<A>(
 
     // Since there is no single sought rank to adapt pivot sampling to, the recommended skewed pivot
     // sampling of dual-pivot Quicksort is used.
-    let sample = array.sample_mut();
+    let mut sample = [0; 5];
+    sample_mut(&mut array, &mut sample);
     let (lower_index, upper_index) = (sample[0], sample[2]); // (0, 1, 2)
 
     // We partition the array with respect to the two pivot values. The pivot values move to the new
@@ -506,3 +475,29 @@ fn _get_many_from_sorted_mut_unchecked<A>(
         upper_values,
     );
 }
+
+/// Equally space `sample` indexes around the center of `array` and sort them by their values.
+///
+/// `sample` content is ignored but its length defines the sample size and the sample space divider.
+///
+/// Assumes arrays of at least `sample.len() + 2` elements.
+pub(crate) fn sample_mut<A, S>(array: &mut ArrayBase<S, Ix1>, sample: &mut [usize])
+where
+    A: Ord + Clone,
+    S: DataMut<Elem = A>,
+{
+    // Space between sample indexes.
+    let space = array.len() / (sample.len() + 2);
+    // Lowermost sample index.
+    let lowermost = array.len() / 2 - (sample.len() / 2) * space;
+    // Equally space sample indexes and sort them by their values by looking up their indexes.
+    for mut index in 1..sample.len() {
+        // Equally space sample indexes based on their lowermost index.
+        sample[index] = lowermost + index * space;
+        // Insertion sort looking up only the already equally spaced sample indexes.
+        while index > 0 && array[sample[index - 1]] > array[sample[index]] {
+            array.swap(sample[index - 1], sample[index]);
+            index -= 1;
+        }
+    }
+}

From 2b1b502bd0d0c2d9971832ecf6b4f2c989c91b01 Mon Sep 17 00:00:00 2001
From: Rouven Spreckels <rs@qu1x.dev>
Date: Wed, 30 Jun 2021 16:04:22 +0200
Subject: [PATCH 05/13] Adapt pivot sampling for bulk version.

---
 src/sort.rs | 85 ++++++++++++++++++++++++++++++++++++++++++++++++-----
 1 file changed, 78 insertions(+), 7 deletions(-)

diff --git a/src/sort.rs b/src/sort.rs
index 54f9d4a5..0796ffee 100644
--- a/src/sort.rs
+++ b/src/sort.rs
@@ -171,10 +171,11 @@ where
             }
             self[i].clone()
         } else {
-            // Adapt pivot sampling to relative sought rank and switch from dual-pivot to
-            // single-pivot partitioning for extreme sought ranks.
+            // Sorted sample of 5 equally spaced elements around the center.
             let mut sample = [0; 5];
             sample_mut(self, &mut sample);
+            // Adapt pivot sampling to relative sought rank and switch from dual-pivot to
+            // single-pivot partitioning for extreme sought ranks.
             let sought_rank = i as f64 / n as f64;
             if (0.036..=0.964).contains(&sought_rank) {
                 let (lower_index, upper_index) = if sought_rank <= 0.5 {
@@ -400,11 +401,80 @@ fn _get_many_from_sorted_mut_unchecked<A>(
         return;
     }
 
-    // Since there is no single sought rank to adapt pivot sampling to, the recommended skewed pivot
-    // sampling of dual-pivot Quicksort is used.
+    // Sorted sample of 5 equally spaced elements around the center.
     let mut sample = [0; 5];
     sample_mut(&mut array, &mut sample);
-    let (lower_index, upper_index) = (sample[0], sample[2]); // (0, 1, 2)
+    let (lower_index, upper_index) = if indexes.len() == 1 {
+        // Adapt pivot sampling to relative sought rank and switch from dual-pivot to single-pivot
+        // partitioning for extreme sought ranks.
+        let sought_rank = indexes[0] as f64 / n as f64;
+        if (0.036..=0.964).contains(&sought_rank) {
+            if sought_rank <= 0.5 {
+                if sought_rank <= 0.153 {
+                    (0, 1) // (0, 0, 3)
+                } else {
+                    (0, 2) // (0, 1, 2)
+                }
+            } else {
+                if sought_rank <= 0.847 {
+                    (2, 4) // (2, 1, 0)
+                } else {
+                    (3, 4) // (3, 0, 0)
+                }
+            }
+        } else {
+            let pivot_index = if sought_rank <= 0.5 {
+                0 // (0, 4)
+            } else {
+                4 // (4, 0)
+            };
+
+            // We partition the array with respect to the pivot value. The pivot value moves to the
+            // new `pivot_index`.
+            //
+            // Elements strictly less than the pivot value have indexes < `pivot_index`.
+            //
+            // Elements greater than or equal the pivot value have indexes > `pivot_index`.
+            let pivot_index = array.partition_mut(sample[pivot_index]);
+            let (found_exact, split_index) = match indexes.binary_search(&pivot_index) {
+                Ok(index) => (true, index),
+                Err(index) => (false, index),
+            };
+            let (lower_indexes, upper_indexes) = indexes.split_at_mut(split_index);
+            let (lower_values, upper_values) = values.split_at_mut(split_index);
+            let (upper_indexes, upper_values) = if found_exact {
+                upper_values[0] = array[pivot_index].clone(); // Write exactly found value.
+                (&mut upper_indexes[1..], &mut upper_values[1..])
+            } else {
+                (upper_indexes, upper_values)
+            };
+
+            // We search recursively for the values corresponding to indexes strictly less than
+            // `pivot_index` in the lower partition.
+            _get_many_from_sorted_mut_unchecked(
+                array.slice_axis_mut(Axis(0), Slice::from(..pivot_index)),
+                lower_indexes,
+                lower_values,
+            );
+
+            // We search recursively for the values corresponding to indexes greater than or equal
+            // `pivot_index` in the upper partition. Since only the upper partition of the array is
+            // passed in, the indexes need to be shifted by length of the lower partition.
+            upper_indexes.iter_mut().for_each(|x| *x -= pivot_index + 1);
+            _get_many_from_sorted_mut_unchecked(
+                array.slice_axis_mut(Axis(0), Slice::from(pivot_index + 1..)),
+                upper_indexes,
+                upper_values,
+            );
+
+            return;
+        }
+    } else {
+        // Since there is no single sought rank to adapt pivot sampling to, the recommended skewed
+        // pivot sampling of dual-pivot Quicksort is used in the assumption that multiple indexes
+        // change characteristics from Quickselect towards Quicksort.
+        (0, 2) // (0, 1, 2)
+    };
 
     // We partition the array with respect to the two pivot values. The pivot values move to the new
     // `lower_index` and `upper_index`.
@@ -414,8 +484,9 @@ fn _get_many_from_sorted_mut_unchecked<A>(
     // Elements greater than or equal the lower pivot value and less than or equal the upper pivot
     // value have indexes > `lower_index` and < `upper_index`.
     //
-    // Elements less than or equal the upper pivot value have indexes > `upper_index`.
-    let (lower_index, upper_index) = array.dual_partition_mut(lower_index, upper_index);
+    // Elements greater than or equal the upper pivot value have indexes > `upper_index`.
+    let (lower_index, upper_index) =
+        array.dual_partition_mut(sample[lower_index], sample[upper_index]);
 
     // We use a divide-and-conquer strategy, splitting the indexes we are searching for (`indexes`)
     // and the corresponding portions of the output slice (`values`) into partitions with respect to

From 78d334bb16ac26a59f2f646a022f783cd43cea4e Mon Sep 17 00:00:00 2001
From: Rouven Spreckels <rs@qu1x.dev>
Date: Sat, 1 Apr 2023 11:37:56 +0200
Subject: [PATCH 06/13] Recursively grow stack on heap whenever necessary.

  * Add test with large array of equal floats.
  * Enable optimization for test profile to reduce execution time.
---
 Cargo.toml    |  4 +++
 src/sort.rs   | 70 ++++++++++++++++++++++++++++++++-------------------
 tests/sort.rs |  8 ++++++
 3 files changed, 56 insertions(+), 26 deletions(-)

diff --git a/Cargo.toml b/Cargo.toml
index 892533e5..7587537e 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -23,6 +23,7 @@ num-traits = "0.2"
 rand = "0.8.3"
 itertools = { version = "0.10.0", default-features = false }
 indexmap = "1.6.2"
+stacker = "0.1.15"
 
 [dev-dependencies]
 ndarray = { version = "0.15.0", features = ["approx"] }
@@ -44,3 +45,6 @@ harness = false
 [[bench]]
 name = "deviation"
 harness = false
+
+[profile.test]
+opt-level = 2
diff --git a/src/sort.rs b/src/sort.rs
index 0796ffee..6d667302 100644
--- a/src/sort.rs
+++ b/src/sort.rs
@@ -1,6 +1,12 @@
 use indexmap::IndexMap;
 use ndarray::prelude::*;
 use ndarray::{Data, DataMut, Slice};
+use stacker::maybe_grow;
+
+/// Guaranteed stack size per recursion step of 1 MiB.
+const RED_ZONE: usize = 1_024 * 1_024;
+/// New stack space of 8 MiB to allocate if within [`RED_ZONE`].
+const STACK_SIZE: usize = 8 * RED_ZONE;
 
 /// Methods for sorting and partitioning 1-D arrays.
 pub trait Sort1dExt<A, S>
@@ -357,7 +363,9 @@ where
     // Since `!indexes.is_empty()` and indexes must be in-bounds, `array` must
     // be non-empty.
     let mut values = vec![array[0].clone(); indexes.len()];
-    _get_many_from_sorted_mut_unchecked(array.view_mut(), &mut indexes.to_owned(), &mut values);
+    maybe_grow(RED_ZONE, STACK_SIZE, || {
+        _get_many_from_sorted_mut_unchecked(array.view_mut(), &mut indexes.to_owned(), &mut values);
+    });
 
     // We convert the vector to a more search-friendly `IndexMap`.
     indexes.iter().cloned().zip(values.into_iter()).collect()
@@ -451,21 +459,25 @@ fn _get_many_from_sorted_mut_unchecked<A>(
 
             // We search recursively for the values corresponding to indexes strictly less than
             // `pivot_index` in the lower partition.
-            _get_many_from_sorted_mut_unchecked(
-                array.slice_axis_mut(Axis(0), Slice::from(..pivot_index)),
-                lower_indexes,
-                lower_values,
-            );
+            maybe_grow(RED_ZONE, STACK_SIZE, || {
+                _get_many_from_sorted_mut_unchecked(
+                    array.slice_axis_mut(Axis(0), Slice::from(..pivot_index)),
+                    lower_indexes,
+                    lower_values,
+                );
+            });
 
             // We search recursively for the values corresponding to indexes greater than or equal
             // `pivot_index` in the upper partition. Since only the upper partition of the array is
             // passed in, the indexes need to be shifted by length of the lower partition.
             upper_indexes.iter_mut().for_each(|x| *x -= pivot_index + 1);
-            _get_many_from_sorted_mut_unchecked(
-                array.slice_axis_mut(Axis(0), Slice::from(pivot_index + 1..)),
-                upper_indexes,
-                upper_values,
-            );
+            maybe_grow(RED_ZONE, STACK_SIZE, || {
+                _get_many_from_sorted_mut_unchecked(
+                    array.slice_axis_mut(Axis(0), Slice::from(pivot_index + 1..)),
+                    upper_indexes,
+                    upper_values,
+                );
+            });
 
             return;
         }
@@ -519,32 +531,38 @@ fn _get_many_from_sorted_mut_unchecked<A>(
 
     // We search recursively for the values corresponding to indexes strictly less than
     // `lower_index` in the lower partition.
-    _get_many_from_sorted_mut_unchecked(
-        array.slice_axis_mut(Axis(0), Slice::from(..lower_index)),
-        lower_indexes,
-        lower_values,
-    );
+    maybe_grow(RED_ZONE, STACK_SIZE, || {
+        _get_many_from_sorted_mut_unchecked(
+            array.slice_axis_mut(Axis(0), Slice::from(..lower_index)),
+            lower_indexes,
+            lower_values,
+        );
+    });
 
     // We search recursively for the values corresponding to indexes greater than or equal
     // `lower_index` in the inner partition, that is between the lower and upper partition. Since
     // only the inner partition of the array is passed in, the indexes need to be shifted by length
     // of the lower partition.
     inner_indexes.iter_mut().for_each(|x| *x -= lower_index + 1);
-    _get_many_from_sorted_mut_unchecked(
-        array.slice_axis_mut(Axis(0), Slice::from(lower_index + 1..upper_index)),
-        inner_indexes,
-        inner_values,
-    );
+    maybe_grow(RED_ZONE, STACK_SIZE, || {
+        _get_many_from_sorted_mut_unchecked(
+            array.slice_axis_mut(Axis(0), Slice::from(lower_index + 1..upper_index)),
+            inner_indexes,
+            inner_values,
+        );
+    });
 
     // We search recursively for the values corresponding to indexes greater than or equal
     // `upper_index` in the upper partition. Since only the upper partition of the array is passed
     // in, the indexes need to be shifted by the combined length of the lower and inner partition.
     upper_indexes.iter_mut().for_each(|x| *x -= upper_index + 1);
-    _get_many_from_sorted_mut_unchecked(
-        array.slice_axis_mut(Axis(0), Slice::from(upper_index + 1..)),
-        upper_indexes,
-        upper_values,
-    );
+    maybe_grow(RED_ZONE, STACK_SIZE, || {
+        _get_many_from_sorted_mut_unchecked(
+            array.slice_axis_mut(Axis(0), Slice::from(upper_index + 1..)),
+            upper_indexes,
+            upper_values,
+        );
+    });
 }
 
 /// Equally space `sample` indexes around the center of `array` and sort them by their values.
diff --git a/tests/sort.rs b/tests/sort.rs
index b298e88d..56c14d66 100644
--- a/tests/sort.rs
+++ b/tests/sort.rs
@@ -1,5 +1,7 @@
 use ndarray::prelude::*;
 use ndarray_stats::Sort1dExt;
+use ndarray_stats::{interpolate::Linear, Quantile1dExt};
+use noisy_float::types::{n64, N64};
 use quickcheck_macros::quickcheck;
 
 #[test]
@@ -63,6 +65,12 @@ fn test_dual_partition_mut() {
     }
 }
 
+#[test]
+fn test_quantile_mut_with_large_array_of_equal_floats() {
+    let mut array: Array1<N64> = Array1::ones(100000);
+    array.quantile_mut(n64(0.5), &Linear).unwrap();
+}
+
 #[test]
 fn test_sorted_get_mut() {
     let a = arr1(&[1, 3, 2, 10]);

From 23782e775a0e47eb973daff423b1ab7ce1d0a468 Mon Sep 17 00:00:00 2001
From: Rouven Spreckels <rs@qu1x.dev>
Date: Sun, 2 Apr 2023 11:04:06 +0200
Subject: [PATCH 07/13] Also protect non-bulk version from stack overflow.

---
 src/sort.rs   | 34 +++++++++++++++++++++-------------
 tests/sort.rs |  2 +-
 2 files changed, 22 insertions(+), 14 deletions(-)

diff --git a/src/sort.rs b/src/sort.rs
index 6d667302..d3d3440a 100644
--- a/src/sort.rs
+++ b/src/sort.rs
@@ -200,18 +200,24 @@ where
                 let (lower_index, upper_index) =
                     self.dual_partition_mut(sample[lower_index], sample[upper_index]);
                 if i < lower_index {
-                    self.slice_axis_mut(Axis(0), Slice::from(..lower_index))
-                        .get_from_sorted_mut(i)
+                    maybe_grow(RED_ZONE, STACK_SIZE, || {
+                        self.slice_axis_mut(Axis(0), Slice::from(..lower_index))
+                            .get_from_sorted_mut(i)
+                    })
                 } else if i == lower_index {
                     self[i].clone()
                 } else if i < upper_index {
-                    self.slice_axis_mut(Axis(0), Slice::from(lower_index + 1..upper_index))
-                        .get_from_sorted_mut(i - (lower_index + 1))
+                    maybe_grow(RED_ZONE, STACK_SIZE, || {
+                        self.slice_axis_mut(Axis(0), Slice::from(lower_index + 1..upper_index))
+                            .get_from_sorted_mut(i - (lower_index + 1))
+                    })
                 } else if i == upper_index {
                     self[i].clone()
                 } else {
-                    self.slice_axis_mut(Axis(0), Slice::from(upper_index + 1..))
-                        .get_from_sorted_mut(i - (upper_index + 1))
+                    maybe_grow(RED_ZONE, STACK_SIZE, || {
+                        self.slice_axis_mut(Axis(0), Slice::from(upper_index + 1..))
+                            .get_from_sorted_mut(i - (upper_index + 1))
+                    })
                 }
             } else {
                 let pivot_index = if sought_rank <= 0.5 {
@@ -221,13 +227,17 @@ where
                 };
                 let pivot_index = self.partition_mut(sample[pivot_index]);
                 if i < pivot_index {
-                    self.slice_axis_mut(Axis(0), Slice::from(..pivot_index))
-                        .get_from_sorted_mut(i)
+                    maybe_grow(RED_ZONE, STACK_SIZE, || {
+                        self.slice_axis_mut(Axis(0), Slice::from(..pivot_index))
+                            .get_from_sorted_mut(i)
+                    })
                 } else if i == pivot_index {
                     self[i].clone()
                 } else {
-                    self.slice_axis_mut(Axis(0), Slice::from(pivot_index + 1..))
-                        .get_from_sorted_mut(i - (pivot_index + 1))
+                    maybe_grow(RED_ZONE, STACK_SIZE, || {
+                        self.slice_axis_mut(Axis(0), Slice::from(pivot_index + 1..))
+                            .get_from_sorted_mut(i - (pivot_index + 1))
+                    })
                 }
             }
         }
@@ -363,9 +373,7 @@ where
     // Since `!indexes.is_empty()` and indexes must be in-bounds, `array` must
     // be non-empty.
     let mut values = vec![array[0].clone(); indexes.len()];
-    maybe_grow(RED_ZONE, STACK_SIZE, || {
-        _get_many_from_sorted_mut_unchecked(array.view_mut(), &mut indexes.to_owned(), &mut values);
-    });
+    _get_many_from_sorted_mut_unchecked(array.view_mut(), &mut indexes.to_owned(), &mut values);
 
     // We convert the vector to a more search-friendly `IndexMap`.
     indexes.iter().cloned().zip(values.into_iter()).collect()
diff --git a/tests/sort.rs b/tests/sort.rs
index 56c14d66..ad4bc5ec 100644
--- a/tests/sort.rs
+++ b/tests/sort.rs
@@ -67,7 +67,7 @@ fn test_dual_partition_mut() {
 
 #[test]
 fn test_quantile_mut_with_large_array_of_equal_floats() {
-    let mut array: Array1<N64> = Array1::ones(100000);
+    let mut array: Array1<N64> = Array1::ones(100_000);
     array.quantile_mut(n64(0.5), &Linear).unwrap();
 }
 

From 359dc0c4cf1f40ee3292db31f8868a64f77fb379 Mon Sep 17 00:00:00 2001
From: Rouven Spreckels <rs@qu1x.dev>
Date: Sun, 2 Apr 2023 12:03:31 +0200
Subject: [PATCH 08/13] Avoid worst case complexity for equal elements.

  * Delegate single-index invocation to non-bulk implementation.
  * Non-bulk implementation skips recursion if `index == pivot`.
---
 src/sort.rs   | 105 ++++++++++----------------------------------------
 tests/sort.rs |   2 +-
 2 files changed, 21 insertions(+), 86 deletions(-)

diff --git a/src/sort.rs b/src/sort.rs
index d3d3440a..d5351fe5 100644
--- a/src/sort.rs
+++ b/src/sort.rs
@@ -366,17 +366,23 @@ where
     A: Ord + Clone,
     S: DataMut<Elem = A>,
 {
-    if indexes.is_empty() {
-        return IndexMap::new();
-    }
-
-    // Since `!indexes.is_empty()` and indexes must be in-bounds, `array` must
-    // be non-empty.
-    let mut values = vec![array[0].clone(); indexes.len()];
-    _get_many_from_sorted_mut_unchecked(array.view_mut(), &mut indexes.to_owned(), &mut values);
+    match indexes.len() {
+        0 => IndexMap::new(),
+        1 => IndexMap::from([(indexes[0], array.get_from_sorted_mut(indexes[0]))]),
+        _ => {
+            // Since `!indexes.is_empty()` and indexes must be in-bounds, `array` must
+            // be non-empty.
+            let mut values = vec![array[0].clone(); indexes.len()];
+            _get_many_from_sorted_mut_unchecked(
+                array.view_mut(),
+                &mut indexes.to_owned(),
+                &mut values,
+            );
 
-    // We convert the vector to a more search-friendly `IndexMap`.
-    indexes.iter().cloned().zip(values.into_iter()).collect()
+            // We convert the vector to a more search-friendly `IndexMap`.
+            indexes.iter().cloned().zip(values.into_iter()).collect()
+        }
+    }
 }
 
 /// This is the recursive portion of `get_many_from_sorted_mut_unchecked`.
@@ -420,81 +426,10 @@ fn _get_many_from_sorted_mut_unchecked<A>(
     // Sorted sample of 5 equally spaced elements around the center.
     let mut sample = [0; 5];
     sample_mut(&mut array, &mut sample);
-    let (lower_index, upper_index) = if indexes.len() == 1 {
-        // Adapt pivot sampling to relative sought rank and switch from dual-pivot to single-pivot
-        // partitioning for extreme sought ranks.
-        let sought_rank = indexes[0] as f64 / n as f64;
-        if (0.036..=0.964).contains(&sought_rank) {
-            if sought_rank <= 0.5 {
-                if sought_rank <= 0.153 {
-                    (0, 1) // (0, 0, 3)
-                } else {
-                    (0, 2) // (0, 1, 2)
-                }
-            } else {
-                if sought_rank <= 0.847 {
-                    (2, 4) // (2, 1, 0)
-                } else {
-                    (3, 4) // (3, 0, 0)
-                }
-            }
-        } else {
-            let pivot_index = if sought_rank <= 0.5 {
-                0 // (0, 4)
-            } else {
-                4 // (4, 0)
-            };
-
-            // We partition the array with respect to the pivot value. The pivot value moves to the
-            // new `pivot_index`.
-            //
-            // Elements strictly less than the pivot value have indexes < `pivot_index`.
-            //
-            // Elements greater than or equal the pivot value have indexes > `pivot_index`.
-            let pivot_index = array.partition_mut(sample[pivot_index]);
-            let (found_exact, split_index) = match indexes.binary_search(&pivot_index) {
-                Ok(index) => (true, index),
-                Err(index) => (false, index),
-            };
-            let (lower_indexes, upper_indexes) = indexes.split_at_mut(split_index);
-            let (lower_values, upper_values) = values.split_at_mut(split_index);
-            let (upper_indexes, upper_values) = if found_exact {
-                upper_values[0] = array[pivot_index].clone(); // Write exactly found value.
-                (&mut upper_indexes[1..], &mut upper_values[1..])
-            } else {
-                (upper_indexes, upper_values)
-            };
-
-            // We search recursively for the values corresponding to indexes strictly less than
-            // `pivot_index` in the lower partition.
-            maybe_grow(RED_ZONE, STACK_SIZE, || {
-                _get_many_from_sorted_mut_unchecked(
-                    array.slice_axis_mut(Axis(0), Slice::from(..pivot_index)),
-                    lower_indexes,
-                    lower_values,
-                );
-            });
-
-            // We search recursively for the values corresponding to indexes greater than or equal
-            // `pivot_index` in the upper partition. Since only the upper partition of the array is
-            // passed in, the indexes need to be shifted by length of the lower partition.
-            upper_indexes.iter_mut().for_each(|x| *x -= pivot_index + 1);
-            maybe_grow(RED_ZONE, STACK_SIZE, || {
-                _get_many_from_sorted_mut_unchecked(
-                    array.slice_axis_mut(Axis(0), Slice::from(pivot_index + 1..)),
-                    upper_indexes,
-                    upper_values,
-                );
-            });
-
-            return;
-        }
-    } else {
-        // Since there is no single sought rank to adapt pivot sampling to, the recommended skewed
-        // pivot sampling of dual-pivot Quicksort is used in the assumption that multiple indexes
-        // change characteristics from Quickselect towards Quicksort.
-        (0, 2) // (0, 1, 2)
-    };
+    // Since there is no single sought rank to adapt pivot sampling to, the recommended skewed
+    // pivot sampling of dual-pivot Quicksort is used in the assumption that multiple indexes
+    // change characteristics from Quickselect towards Quicksort.
+    let (lower_index, upper_index) = (0, 2); // (0, 1, 2)
 
     // We partition the array with respect to the two pivot values. The pivot values move to the new
     // `lower_index` and `upper_index`.
diff --git a/tests/sort.rs b/tests/sort.rs
index ad4bc5ec..093fb0f6 100644
--- a/tests/sort.rs
+++ b/tests/sort.rs
@@ -67,7 +67,7 @@ fn test_dual_partition_mut() {
 
 #[test]
 fn test_quantile_mut_with_large_array_of_equal_floats() {
-    let mut array: Array1<N64> = Array1::ones(100_000);
+    let mut array: Array1<N64> = Array1::ones(1_000_000);
     array.quantile_mut(n64(0.5), &Linear).unwrap();
 }
 

From 4c3ae67468a3cd5befc8ea3bebb69c25e0ee0dc7 Mon Sep 17 00:00:00 2001
From: Rouven Spreckels <rs@qu1x.dev>
Date: Sun, 2 Apr 2023 14:01:51 +0200
Subject: [PATCH 09/13] Fix sample, allow single-pivot in bulk mode again.

  * Move single-index delegation into recursion allowing single-pivoting
    in bulk mode.
  * Avoid worst case complexity in non-bulk and bulk mode by using dual-pivot
    if both sample choices of single-pivot are equal.
  * Fix first sample index from being skipped.
---
 src/sort.rs | 34 ++++++++++++++++------------------
 1 file changed, 16 insertions(+), 18 deletions(-)

diff --git a/src/sort.rs b/src/sort.rs
index d5351fe5..84bba799 100644
--- a/src/sort.rs
+++ b/src/sort.rs
@@ -183,7 +183,7 @@ where
             // Adapt pivot sampling to relative sought rank and switch from dual-pivot to
             // single-pivot partitioning for extreme sought ranks.
             let sought_rank = i as f64 / n as f64;
-            if (0.036..=0.964).contains(&sought_rank) {
+            if (0.036..=0.964).contains(&sought_rank) || self[sample[0]] == self[sample[4]] {
                 let (lower_index, upper_index) = if sought_rank <= 0.5 {
                     if sought_rank <= 0.153 {
                         (0, 1) // (0, 0, 3)
@@ -366,23 +366,17 @@ where
     A: Ord + Clone,
     S: DataMut<Elem = A>,
 {
-    match indexes.len() {
-        0 => IndexMap::new(),
-        1 => IndexMap::from([(indexes[0], array.get_from_sorted_mut(indexes[0]))]),
-        _ => {
-            // Since `!indexes.is_empty()` and indexes must be in-bounds, `array` must
-            // be non-empty.
-            let mut values = vec![array[0].clone(); indexes.len()];
-            _get_many_from_sorted_mut_unchecked(
-                array.view_mut(),
-                &mut indexes.to_owned(),
-                &mut values,
-            );
-
-            // We convert the vector to a more search-friendly `IndexMap`.
-            indexes.iter().cloned().zip(values.into_iter()).collect()
-        }
+    if indexes.is_empty() {
+        return IndexMap::new();
     }
+
+    // Since `!indexes.is_empty()` and indexes must be in-bounds, `array` must
+    // be non-empty.
+    let mut values = vec![array[0].clone(); indexes.len()];
+    _get_many_from_sorted_mut_unchecked(array.view_mut(), &mut indexes.to_owned(), &mut values);
+
+    // We convert the vector to a more search-friendly `IndexMap`.
+    indexes.iter().cloned().zip(values.into_iter()).collect()
 }
 
 /// This is the recursive portion of `get_many_from_sorted_mut_unchecked`.
@@ -408,6 +402,10 @@ fn _get_many_from_sorted_mut_unchecked<A>(
         // Nothing to do in this case.
         return;
     }
+    if indexes.len() == 1 {
+        values[0] = array.get_from_sorted_mut(indexes[0]);
+        return;
+    }
 
     // Recursion cutoff at integer multiple of sample space divider of 7 elements.
     if n < 21 {
@@ -523,7 +521,7 @@ where
     // Lowermost sample index.
     let lowermost = array.len() / 2 - (sample.len() / 2) * space;
     // Equally space sample indexes and sort them by their values by looking up their indexes.
-    for mut index in 1..sample.len() {
+    for mut index in 0..sample.len() {
         // Equally space sample indexes based on their lowermost index.
         sample[index] = lowermost + index * space;
         // Insertion sort looking up only the already equally spaced sample indexes.

From 05b87b4471f15e68a51380e0090898335661f782 Mon Sep 17 00:00:00 2001
From: Rouven Spreckels <rs@qu1x.dev>
Date: Mon, 3 Apr 2023 08:22:38 +0200
Subject: [PATCH 10/13] Make sampling public and test it.

---
 src/sort.rs   | 88 +++++++++++++++++++++++++++++++++++----------------
 tests/sort.rs | 47 +++++++++++++++++++++++++++
 2 files changed, 107 insertions(+), 28 deletions(-)

diff --git a/src/sort.rs b/src/sort.rs
index 84bba799..26918c49 100644
--- a/src/sort.rs
+++ b/src/sort.rs
@@ -154,6 +154,41 @@ where
         A: Ord + Clone,
         S: DataMut;
 
+    /// Equally spaces `sample` indexes around the center of `array` and sorts them by their values.
+    ///
+    /// `sample` content is ignored but its length defines the sample size and sample space divider.
+    ///
+    /// **Panics** if `self.len() < sample.len() + 2`.
+    ///
+    /// ```
+    /// use ndarray::Array1;
+    /// use ndarray_stats::Sort1dExt;
+    /// use noisy_float::types::n64;
+    ///
+    /// // Define sample size of 5.
+    /// let mut sample = [0; 5];
+    ///
+    /// // Create array from 100 down to 1;
+    /// let mut array = Array1::range(n64(100.0), n64(0.0), n64(-1.0));
+    /// assert_eq!(array.len(), 100);
+    ///
+    /// // Find `sample` indices and sort their values `array[sample[s]]`.
+    /// array.sample_mut(&mut sample);
+    ///
+    /// // Equally spaced by 13 non-sample elements around center of 50.
+    /// let assert = [22, 36, 50, 64, 78];
+    /// assert_eq!(sample, &assert[..]);
+    /// // Ensure `array[sample[s]]` is sorted.
+    /// assert_eq!(sample.map(|s| array[s]), assert.map(|s| s as f64));
+    /// // Ensure reverse order of non-sample elements is preserved.
+    /// assert_eq!(array[49], n64(51.0));
+    /// assert_eq!(array[51], n64(49.0));
+    /// ```
+    fn sample_mut(&mut self, sample: &mut [usize])
+    where
+        A: Ord + Clone,
+        S: DataMut;
+
     private_decl! {}
 }
 
@@ -179,7 +214,7 @@ where
         } else {
             // Sorted sample of 5 equally spaced elements around the center.
             let mut sample = [0; 5];
-            sample_mut(self, &mut sample);
+            self.sample_mut(&mut sample);
             // Adapt pivot sampling to relative sought rank and switch from dual-pivot to
             // single-pivot partitioning for extreme sought ranks.
             let sought_rank = i as f64 / n as f64;
@@ -344,6 +379,29 @@ where
         (lower, upper)
     }
 
+    fn sample_mut(&mut self, sample: &mut [usize])
+    where
+        A: Ord + Clone,
+        S: DataMut,
+    {
+        // Assumes arrays of at least `sample.len() + 2` elements.
+        assert!(self.len() >= sample.len() + 2);
+        // Space between sample indexes.
+        let space = self.len() / (sample.len() + 2);
+        // Lowermost sample index.
+        let lowermost = self.len() / 2 - (sample.len() / 2) * space;
+        // Equally space sample indexes and sort them by their values by looking up their indexes.
+        for mut index in 0..sample.len() {
+            // Equally space sample indexes based on their lowermost index.
+            sample[index] = lowermost + index * space;
+            // Insertion sort looking up only the already equally spaced sample indexes.
+            while index > 0 && self[sample[index - 1]] > self[sample[index]] {
+                self.swap(sample[index - 1], sample[index]);
+                index -= 1;
+            }
+        }
+    }
+
     private_impl! {}
 }
 
@@ -423,7 +481,7 @@ fn _get_many_from_sorted_mut_unchecked<A>(
 
     // Sorted sample of 5 equally spaced elements around the center.
     let mut sample = [0; 5];
-    sample_mut(&mut array, &mut sample);
+    array.sample_mut(&mut sample);
     // Since there is no single sought rank to adapt pivot sampling to, the recommended skewed
     // pivot sampling of dual-pivot Quicksort is used in the assumption that multiple indexes
     // change characteristics from Quickselect towards Quicksort.
@@ -505,29 +563,3 @@ fn _get_many_from_sorted_mut_unchecked<A>(
         );
     });
 }
-
-/// Equally space `sample` indexes around the center of `array` and sort them by their values.
-///
-/// `sample` content is ignored but its length defines the sample size and the sample space divider.
-///
-/// Assumes arrays of at least `sample.len() + 2` elements.
-pub(crate) fn sample_mut<A, S>(array: &mut ArrayBase<S, Ix1>, sample: &mut [usize])
-where
-    A: Ord + Clone,
-    S: DataMut<Elem = A>,
-{
-    // Space between sample indexes.
-    let space = array.len() / (sample.len() + 2);
-    // Lowermost sample index.
-    let lowermost = array.len() / 2 - (sample.len() / 2) * space;
-    // Equally space sample indexes and sort them by their values by looking up their indexes.
-    for mut index in 0..sample.len() {
-        // Equally space sample indexes based on their lowermost index.
-        sample[index] = lowermost + index * space;
-        // Insertion sort looking up only the already equally spaced sample indexes.
-        while index > 0 && array[sample[index - 1]] > array[sample[index]] {
-            array.swap(sample[index - 1], sample[index]);
-            index -= 1;
-        }
-    }
-}
diff --git a/tests/sort.rs b/tests/sort.rs
index 093fb0f6..09654abb 100644
--- a/tests/sort.rs
+++ b/tests/sort.rs
@@ -65,6 +65,53 @@ fn test_dual_partition_mut() {
     }
 }
 
+#[test]
+fn test_sample_mut() {
+    let mut sample = [0; 5];
+
+    let assert = [1, 2, 3, 4, 5];
+    let mut array = Array1::range(n64(0.0), n64(7.0), n64(1.0));
+    assert_eq!(array.len(), 7);
+    array.sample_mut(&mut sample);
+    assert_eq!(sample, &assert[..]);
+
+    let assert = [2, 3, 4, 5, 6];
+    let mut array = Array1::range(n64(0.0), n64(8.0), n64(1.0));
+    assert_eq!(array.len(), 8);
+    array.sample_mut(&mut sample);
+    assert_eq!(sample, &assert[..]);
+    let mut array = Array1::range(n64(0.0), n64(9.0), n64(1.0));
+    assert_eq!(array.len(), 9);
+    array.sample_mut(&mut sample);
+    assert_eq!(sample, &assert[..]);
+
+    let assert = [3, 4, 5, 6, 7];
+    let mut array = Array1::range(n64(0.0), n64(10.0), n64(1.0));
+    assert_eq!(array.len(), 10);
+    array.sample_mut(&mut sample);
+    assert_eq!(sample, &assert[..]);
+
+    let assert = [3, 5, 7, 9, 11];
+    let mut array = Array1::range(n64(0.0), n64(14.0), n64(1.0));
+    assert_eq!(array.len(), 14);
+    array.sample_mut(&mut sample);
+    assert_eq!(sample, &assert[..]);
+
+    let assert = [22, 36, 50, 64, 78];
+    let mut array = Array1::range(n64(0.0), n64(100.0), n64(1.0));
+    assert_eq!(array.len(), 100);
+    array.sample_mut(&mut sample);
+    assert_eq!(sample, &assert[..]);
+    assert_eq!(sample.map(|s| array[s]), assert.map(|s| s as f64));
+
+    let assert = [22, 36, 50, 64, 78];
+    let mut array = Array1::range(n64(100.0), n64(0.0), n64(-1.0));
+    assert_eq!(array.len(), 100);
+    array.sample_mut(&mut sample);
+    assert_eq!(sample, &assert[..]);
+    assert_eq!(sample.map(|s| array[s]), assert.map(|s| s as f64));
+}
+
 #[test]
 fn test_quantile_mut_with_large_array_of_equal_floats() {
     let mut array: Array1<N64> = Array1::ones(1_000_000);

From 34952ceddb767706f77815f7483e51bf48e5125a Mon Sep 17 00:00:00 2001
From: Rouven Spreckels <rs@qu1x.dev>
Date: Mon, 3 Apr 2023 08:53:24 +0200
Subject: [PATCH 11/13] Add test with large (rev) sorted array of floats.

---
 tests/sort.rs | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/tests/sort.rs b/tests/sort.rs
index 09654abb..7f13a78e 100644
--- a/tests/sort.rs
+++ b/tests/sort.rs
@@ -118,6 +118,18 @@ fn test_quantile_mut_with_large_array_of_equal_floats() {
     array.quantile_mut(n64(0.5), &Linear).unwrap();
 }
 
+#[test]
+fn test_quantile_mut_with_large_array_of_sorted_floats() {
+    let mut array: Array1<N64> = Array1::range(n64(0.0), n64(1_000_000.0), n64(1.0));
+    array.quantile_mut(n64(0.5), &Linear).unwrap();
+}
+
+#[test]
+fn test_quantile_mut_with_large_array_of_rev_sorted_floats() {
+    let mut array: Array1<N64> = Array1::range(n64(1_000_000.0), n64(0.0), n64(-1.0));
+    array.quantile_mut(n64(0.5), &Linear).unwrap();
+}
+
 #[test]
 fn test_sorted_get_mut() {
     let a = arr1(&[1, 3, 2, 10]);

From 4dc45a4806fd7e0c0feeb28ed7e01bfcb31df389 Mon Sep 17 00:00:00 2001
From: Rouven Spreckels <rs@qu1x.dev>
Date: Mon, 3 Apr 2023 09:08:05 +0200
Subject: [PATCH 12/13] Add tests and enlarge arrays.

---
 tests/sort.rs | 27 ++++++++++++++++++++++++---
 1 file changed, 24 insertions(+), 3 deletions(-)

diff --git a/tests/sort.rs b/tests/sort.rs
index 7f13a78e..7657704d 100644
--- a/tests/sort.rs
+++ b/tests/sort.rs
@@ -114,22 +114,43 @@ fn test_sample_mut() {
 
 #[test]
 fn test_quantile_mut_with_large_array_of_equal_floats() {
-    let mut array: Array1<N64> = Array1::ones(1_000_000);
+    let mut array: Array1<N64> = Array1::ones(10_000_000);
     array.quantile_mut(n64(0.5), &Linear).unwrap();
 }
 
 #[test]
 fn test_quantile_mut_with_large_array_of_sorted_floats() {
-    let mut array: Array1<N64> = Array1::range(n64(0.0), n64(1_000_000.0), n64(1.0));
+    let mut array: Array1<N64> = Array1::range(n64(0.0), n64(1e7), n64(1.0));
     array.quantile_mut(n64(0.5), &Linear).unwrap();
 }
 
 #[test]
 fn test_quantile_mut_with_large_array_of_rev_sorted_floats() {
-    let mut array: Array1<N64> = Array1::range(n64(1_000_000.0), n64(0.0), n64(-1.0));
+    let mut array: Array1<N64> = Array1::range(n64(1e7), n64(0.0), n64(-1.0));
     array.quantile_mut(n64(0.5), &Linear).unwrap();
 }
 
+#[test]
+fn test_quantiles_mut_with_large_array_of_equal_floats() {
+    let mut array: Array1<N64> = Array1::ones(10_000_000);
+    let quantiles: Array1<N64> = Array1::range(n64(0.0), n64(1.0), n64(1e-6));
+    array.quantiles_mut(&quantiles, &Linear).unwrap();
+}
+
+#[test]
+fn test_quantiles_mut_with_large_array_of_sorted_floats() {
+    let mut array: Array1<N64> = Array1::range(n64(0.0), n64(1e7), n64(1.0));
+    let quantiles: Array1<N64> = Array1::range(n64(0.0), n64(1.0), n64(1e-6));
+    array.quantiles_mut(&quantiles, &Linear).unwrap();
+}
+
+#[test]
+fn test_quantiles_mut_with_large_array_of_rev_sorted_floats() {
+    let mut array: Array1<N64> = Array1::range(n64(1e7), n64(0.0), n64(-1.0));
+    let quantiles: Array1<N64> = Array1::range(n64(0.0), n64(1.0), n64(1e-6));
+    array.quantiles_mut(&quantiles, &Linear).unwrap();
+}
+
 #[test]
 fn test_sorted_get_mut() {
     let a = arr1(&[1, 3, 2, 10]);

From 0821516f551b3d2f02d83de8c4634a46ebd36241 Mon Sep 17 00:00:00 2001
From: Rouven Spreckels <rs@qu1x.dev>
Date: Mon, 3 Apr 2023 10:11:47 +0200
Subject: [PATCH 13/13] Bump MSRV to 1.56 for `array_map`/`rust-version`.

---
 .github/workflows/ci.yml | 2 +-
 Cargo.toml               | 1 +
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index c651a39b..2db250fa 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -20,7 +20,7 @@ jobs:
           - stable
           - beta
           - nightly
-          - 1.49.0  # MSRV
+          - 1.56.0  # MSRV
     steps:
       - uses: actions/checkout@v2
       - uses: actions-rs/toolchain@v1
diff --git a/Cargo.toml b/Cargo.toml
index 7587537e..f0537f58 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -3,6 +3,7 @@ name = "ndarray-stats"
 version = "0.5.1"
 authors = ["Jim Turner <ndarray-stats@turner.link>", "LukeMathWalker <rust@lpalmieri.com>"]
 edition = "2018"
+rust-version = "1.56"
 
 license = "MIT/Apache-2.0"