rust-lang-ci · Oct 12, 2024
diff --git a/‎library/alloc/src/slice.rs
+10-14 b/‎library/alloc/src/slice.rs
+10-14
diff --git a/‎library/alloc/src/slice/tests.rs
-369 b/‎library/alloc/src/slice/tests.rs
-369
diff --git a/‎library/alloc/tests/lib.rs
+2 b/‎library/alloc/tests/lib.rs
+2
diff --git a/‎library/alloc/tests/sort/ffi_types.rs
+82 b/‎library/alloc/tests/sort/ffi_types.rs
+82
diff --git a/‎library/alloc/tests/sort/known_good_stable_sort.rs
+192 b/‎library/alloc/tests/sort/known_good_stable_sort.rs
+192
diff --git a/‎library/alloc/tests/sort/mod.rs
+17 b/‎library/alloc/tests/sort/mod.rs
+17
diff --git a/‎library/alloc/tests/sort/patterns.rs
+211 b/‎library/alloc/tests/sort/patterns.rs
+211
diff --git a/‎library/alloc/tests/sort/tests.rs
+1,233 b/‎library/alloc/tests/sort/tests.rs
+1,233
diff --git a/‎library/alloc/tests/sort/zipf.rs
+208 b/‎library/alloc/tests/sort/zipf.rs
+208
diff --git a/‎library/core/tests/slice.rs
-51 b/‎library/core/tests/slice.rs
-51
@@ -19,20 +19,6 @@ use core::cmp::Ordering::{self, Less};
 use core::mem::{self, MaybeUninit};
 #[cfg(not(no_global_oom_handling))]
 use core::ptr;
-#[cfg(not(no_global_oom_handling))]
-use core::slice::sort;
-
-use crate::alloc::Allocator;
-#[cfg(not(no_global_oom_handling))]
-use crate::alloc::Global;
-#[cfg(not(no_global_oom_handling))]
-use crate::borrow::ToOwned;
-use crate::boxed::Box;
-use crate::vec::Vec;
-
-#[cfg(test)]
-mod tests;
-
 #[unstable(feature = "array_chunks", issue = "74985")]
 pub use core::slice::ArrayChunks;
 #[unstable(feature = "array_chunks", issue = "74985")]
@@ -43,6 +29,8 @@ pub use core::slice::ArrayWindows;
 pub use core::slice::EscapeAscii;
 #[stable(feature = "slice_get_slice", since = "1.28.0")]
 pub use core::slice::SliceIndex;
+#[cfg(not(no_global_oom_handling))]
+use core::slice::sort;
 #[stable(feature = "slice_group_by", since = "1.77.0")]
 pub use core::slice::{ChunkBy, ChunkByMut};
 #[stable(feature = "rust1", since = "1.0.0")]
@@ -83,6 +71,14 @@ pub use hack::into_vec;
 #[cfg(test)]
 pub use hack::to_vec;
 
+use crate::alloc::Allocator;
+#[cfg(not(no_global_oom_handling))]
+use crate::alloc::Global;
+#[cfg(not(no_global_oom_handling))]
+use crate::borrow::ToOwned;
+use crate::boxed::Box;
+use crate::vec::Vec;
+
 // HACK(japaric): With cfg(test) `impl [T]` is not available, these three
 // functions are actually methods that are in `impl [T]` but not in
 // `core::slice::SliceExt` - we need to supply these functions for the
 
@@ -40,6 +40,7 @@
 #![feature(local_waker)]
 #![feature(vec_pop_if)]
 #![feature(unique_rc_arc)]
+#![feature(macro_metavar_expr_concat)]
 #![allow(internal_features)]
 #![deny(fuzzy_provenance_casts)]
 #![deny(unsafe_op_in_unsafe_fn)]
@@ -59,6 +60,7 @@ mod heap;
 mod linked_list;
 mod rc;
 mod slice;
+mod sort;
 mod str;
 mod string;
 mod task;
 
@@ -0,0 +1,82 @@
+use std::cmp::Ordering;
+
+// Very large stack value.
+#[repr(C)]
+#[derive(PartialEq, Eq, Debug, Clone)]
+pub struct FFIOneKibiByte {
+    values: [i64; 128],
+}
+
+impl FFIOneKibiByte {
+    pub fn new(val: i32) -> Self {
+        let mut values = [0i64; 128];
+        let mut val_i64 = val as i64;
+
+        for elem in &mut values {
+            *elem = val_i64;
+            val_i64 = std::hint::black_box(val_i64 + 1);
+        }
+        Self { values }
+    }
+
+    fn as_i64(&self) -> i64 {
+        self.values[11] + self.values[55] + self.values[77]
+    }
+}
+
+impl PartialOrd for FFIOneKibiByte {
+    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
+        Some(self.cmp(other))
+    }
+}
+
+impl Ord for FFIOneKibiByte {
+    fn cmp(&self, other: &Self) -> Ordering {
+        self.as_i64().cmp(&other.as_i64())
+    }
+}
+
+// 16 byte stack value, with more expensive comparison.
+#[repr(C)]
+#[derive(PartialEq, Debug, Clone, Copy)]
+pub struct F128 {
+    x: f64,
+    y: f64,
+}
+
+impl F128 {
+    pub fn new(val: i32) -> Self {
+        let val_f = (val as f64) + (i32::MAX as f64) + 10.0;
+
+        let x = val_f + 0.1;
+        let y = val_f.log(4.1);
+
+        assert!(y < x);
+        assert!(x.is_normal() && y.is_normal());
+
+        Self { x, y }
+    }
+}
+
+// This is kind of hacky, but we know we only have normal comparable floats in there.
+impl Eq for F128 {}
+
+impl PartialOrd for F128 {
+    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
+        Some(self.cmp(other))
+    }
+}
+
+// Goal is similar code-gen between Rust and C++
+// - Rust https://godbolt.org/z/3YM3xenPP
+// - C++ https://godbolt.org/z/178M6j1zz
+impl Ord for F128 {
+    fn cmp(&self, other: &Self) -> Ordering {
+        // Simulate expensive comparison function.
+        let this_div = self.x / self.y;
+        let other_div = other.x / other.y;
+
+        // SAFETY: We checked in the ctor that both are normal.
+        unsafe { this_div.partial_cmp(&other_div).unwrap_unchecked() }
+    }
+}
@@ -0,0 +1,192 @@
+// This module implements a known good stable sort implementation that helps provide better error
+// messages when the correctness tests fail, we can't use the stdlib sort functions because we are
+// testing them for correctness.
+//
+// Based on https://github.com/voultapher/tiny-sort-rs.
+
+use alloc::alloc::{Layout, alloc, dealloc};
+use std::{mem, ptr};
+
+/// Sort `v` preserving initial order of equal elements.
+///
+/// - Guaranteed O(N * log(N)) worst case perf
+/// - No adaptiveness
+/// - Branch miss-prediction not affected by outcome of comparison function
+/// - Uses `v.len()` auxiliary memory.
+///
+/// If `T: Ord` does not implement a total order the resulting order is
+/// unspecified. All original elements will remain in `v` and any possible modifications via
+/// interior mutability will be observable. Same is true if `T: Ord` panics.
+///
+/// Panics if allocating the auxiliary memory fails.
+#[inline(always)]
+pub fn sort<T: Ord>(v: &mut [T]) {
+    stable_sort(v, |a, b| a.lt(b))
+}
+
+#[inline(always)]
+fn stable_sort<T, F: FnMut(&T, &T) -> bool>(v: &mut [T], mut is_less: F) {
+    if mem::size_of::<T>() == 0 {
+        return;
+    }
+
+    let len = v.len();
+
+    // Inline the check for len < 2. This happens a lot, instrumenting the Rust compiler suggests
+    // len < 2 accounts for 94% of its calls to `slice::sort`.
+    if len < 2 {
+        return;
+    }
+
+    // SAFETY: We checked that len is > 0 and that T is not a ZST.
+    unsafe {
+        mergesort_main(v, &mut is_less);
+    }
+}
+
+/// The core logic should not be inlined.
+///
+/// SAFETY: The caller has to ensure that len is > 0 and that T is not a ZST.
+#[inline(never)]
+unsafe fn mergesort_main<T, F: FnMut(&T, &T) -> bool>(v: &mut [T], is_less: &mut F) {
+    // While it would be nice to have a merge implementation that only requires N / 2 auxiliary
+    // memory. Doing so would make the merge implementation significantly more complex and
+
+    // SAFETY: See function safety description.
+    let buf = unsafe { BufGuard::new(v.len()) };
+
+    // SAFETY: `scratch` has space for `v.len()` writes. And does not alias `v`.
+    unsafe {
+        mergesort_core(v, buf.buf_ptr.as_ptr(), is_less);
+    }
+}
+
+/// Tiny recursive top-down merge sort optimized for binary size. It has no adaptiveness whatsoever,
+/// no run detection, etc.
+///
+/// Buffer as pointed to by `scratch` must have space for `v.len()` writes. And must not alias `v`.
+#[inline(always)]
+unsafe fn mergesort_core<T, F: FnMut(&T, &T) -> bool>(
+    v: &mut [T],
+    scratch_ptr: *mut T,
+    is_less: &mut F,
+) {
+    let len = v.len();
+
+    if len > 2 {
+        // SAFETY: `mid` is guaranteed in-bounds. And caller has to ensure that `scratch_ptr` can
+        // hold `v.len()` values.
+        unsafe {
+            let mid = len / 2;
+            // Sort the left half recursively.
+            mergesort_core(v.get_unchecked_mut(..mid), scratch_ptr, is_less);
+            // Sort the right half recursively.
+            mergesort_core(v.get_unchecked_mut(mid..), scratch_ptr, is_less);
+            // Combine the two halves.
+            merge(v, scratch_ptr, is_less, mid);
+        }
+    } else if len == 2 {
+        if is_less(&v[1], &v[0]) {
+            v.swap(0, 1);
+        }
+    }
+}
+
+/// Branchless merge function.
+///
+/// SAFETY: The caller must ensure that `scratch_ptr` is valid for `v.len()` writes. And that mid is
+/// in-bounds.
+#[inline(always)]
+unsafe fn merge<T, F>(v: &mut [T], scratch_ptr: *mut T, is_less: &mut F, mid: usize)
+where
+    F: FnMut(&T, &T) -> bool,
+{
+    let len = v.len();
+    debug_assert!(mid > 0 && mid < len);
+
+    let len = v.len();
+
+    // Indexes to track the positions while merging.
+    let mut l = 0;
+    let mut r = mid;
+
+    // SAFETY: No matter what the result of is_less is we check that l and r remain in-bounds and if
+    // is_less panics the original elements remain in `v`.
+    unsafe {
+        let arr_ptr = v.as_ptr();
+
+        for i in 0..len {
+            let left_ptr = arr_ptr.add(l);
+            let right_ptr = arr_ptr.add(r);
+
+            let is_lt = !is_less(&*right_ptr, &*left_ptr);
+            let copy_ptr = if is_lt { left_ptr } else { right_ptr };
+            ptr::copy_nonoverlapping(copy_ptr, scratch_ptr.add(i), 1);
+
+            l += is_lt as usize;
+            r += !is_lt as usize;
+
+            // As long as neither side is exhausted merge left and right elements.
+            if ((l == mid) as u8 + (r == len) as u8) != 0 {
+                break;
+            }
+        }
+
+        // The left or right side is exhausted, drain the right side in one go.
+        let copy_ptr = if l == mid { arr_ptr.add(r) } else { arr_ptr.add(l) };
+        let i = l + (r - mid);
+        ptr::copy_nonoverlapping(copy_ptr, scratch_ptr.add(i), len - i);
+
+        // Now that scratch_ptr holds the full merged content, write it back on-top of v.
+        ptr::copy_nonoverlapping(scratch_ptr, v.as_mut_ptr(), len);
+    }
+}
+
+// SAFETY: The caller has to ensure that Option is Some, UB otherwise.
+unsafe fn unwrap_unchecked<T>(opt_val: Option<T>) -> T {
+    match opt_val {
+        Some(val) => val,
+        None => {
+            // SAFETY: See function safety description.
+            unsafe {
+                core::hint::unreachable_unchecked();
+            }
+        }
+    }
+}
+
+// Extremely basic versions of Vec.
+// Their use is super limited and by having the code here, it allows reuse between the sort
+// implementations.
+struct BufGuard<T> {
+    buf_ptr: ptr::NonNull<T>,
+    capacity: usize,
+}
+
+impl<T> BufGuard<T> {
+    // SAFETY: The caller has to ensure that len is not 0 and that T is not a ZST.
+    unsafe fn new(len: usize) -> Self {
+        debug_assert!(len > 0 && mem::size_of::<T>() > 0);
+
+        // SAFETY: See function safety description.
+        let layout = unsafe { unwrap_unchecked(Layout::array::<T>(len).ok()) };
+
+        // SAFETY: We checked that T is not a ZST.
+        let buf_ptr = unsafe { alloc(layout) as *mut T };
+
+        if buf_ptr.is_null() {
+            panic!("allocation failure");
+        }
+
+        Self { buf_ptr: ptr::NonNull::new(buf_ptr).unwrap(), capacity: len }
+    }
+}
+
+impl<T> Drop for BufGuard<T> {
+    fn drop(&mut self) {
+        // SAFETY: We checked that T is not a ZST.
+        unsafe {
+            dealloc(self.buf_ptr.as_ptr() as *mut u8, Layout::array::<T>(self.capacity).unwrap());
+        }
+    }
+}
@@ -0,0 +1,17 @@
+pub trait Sort {
+    fn name() -> String;
+
+    fn sort<T>(v: &mut [T])
+    where
+        T: Ord;
+
+    fn sort_by<T, F>(v: &mut [T], compare: F)
+    where
+        F: FnMut(&T, &T) -> std::cmp::Ordering;
+}
+
+mod ffi_types;
+mod known_good_stable_sort;
+mod patterns;
+mod tests;
+mod zipf;
@@ -0,0 +1,211 @@
+use std::env;
+use std::hash::Hash;
+use std::str::FromStr;
+use std::sync::OnceLock;
+
+use rand::prelude::*;
+use rand_xorshift::XorShiftRng;
+
+use crate::sort::zipf::ZipfDistribution;
+
+/// Provides a set of patterns useful for testing and benchmarking sorting algorithms.
+/// Currently limited to i32 values.
+
+// --- Public ---
+
+pub fn random(len: usize) -> Vec<i32> {
+    //     .
+    // : . : :
+    // :.:::.::
+
+    random_vec(len)
+}
+
+pub fn random_uniform<R>(len: usize, range: R) -> Vec<i32>
+where
+    R: Into<rand::distributions::Uniform<i32>> + Hash,
+{
+    // :.:.:.::
+
+    let mut rng: XorShiftRng = rand::SeedableRng::seed_from_u64(get_or_init_rand_seed());
+
+    // Abstracting over ranges in Rust :(
+    let dist: rand::distributions::Uniform<i32> = range.into();
+    (0..len).map(|_| dist.sample(&mut rng)).collect()
+}
+
+pub fn random_zipf(len: usize, exponent: f64) -> Vec<i32> {
+    // https://en.wikipedia.org/wiki/Zipf's_law
+
+    let mut rng: XorShiftRng = rand::SeedableRng::seed_from_u64(get_or_init_rand_seed());
+
+    // Abstracting over ranges in Rust :(
+    let dist = ZipfDistribution::new(len, exponent).unwrap();
+    (0..len).map(|_| dist.sample(&mut rng) as i32).collect()
+}
+
+pub fn random_sorted(len: usize, sorted_percent: f64) -> Vec<i32> {
+    //     .:
+    //   .:::. :
+    // .::::::.::
+    // [----][--]
+    //  ^      ^
+    //  |      |
+    // sorted  |
+    //     unsorted
+
+    // Simulate pre-existing sorted slice, where len - sorted_percent are the new unsorted values
+    // and part of the overall distribution.
+    let mut v = random_vec(len);
+    let sorted_len = ((len as f64) * (sorted_percent / 100.0)).round() as usize;
+
+    v[0..sorted_len].sort_unstable();
+
+    v
+}
+
+pub fn all_equal(len: usize) -> Vec<i32> {
+    // ......
+    // ::::::
+
+    (0..len).map(|_| 66).collect::<Vec<_>>()
+}
+
+pub fn ascending(len: usize) -> Vec<i32> {
+    //     .:
+    //   .:::
+    // .:::::
+
+    (0..len as i32).collect::<Vec<_>>()
+}
+
+pub fn descending(len: usize) -> Vec<i32> {
+    // :.
+    // :::.
+    // :::::.
+
+    (0..len as i32).rev().collect::<Vec<_>>()
+}
+
+pub fn saw_mixed(len: usize, saw_count: usize) -> Vec<i32> {
+    // :.  :.    .::.    .:
+    // :::.:::..::::::..:::
+
+    if len == 0 {
+        return Vec::new();
+    }
+
+    let mut vals = random_vec(len);
+    let chunks_size = len / saw_count.max(1);
+    let saw_directions = random_uniform((len / chunks_size) + 1, 0..=1);
+
+    for (i, chunk) in vals.chunks_mut(chunks_size).enumerate() {
+        if saw_directions[i] == 0 {
+            chunk.sort_unstable();
+        } else if saw_directions[i] == 1 {
+            chunk.sort_unstable_by_key(|&e| std::cmp::Reverse(e));
+        } else {
+            unreachable!();
+        }
+    }
+
+    vals
+}
+
+pub fn saw_mixed_range(len: usize, range: std::ops::Range<usize>) -> Vec<i32> {
+    //     :.
+    // :.  :::.    .::.      .:
+    // :::.:::::..::::::..:.:::
+
+    // ascending and descending randomly picked, with length in `range`.
+
+    if len == 0 {
+        return Vec::new();
+    }
+
+    let mut vals = random_vec(len);
+
+    let max_chunks = len / range.start;
+    let saw_directions = random_uniform(max_chunks + 1, 0..=1);
+    let chunk_sizes = random_uniform(max_chunks + 1, (range.start as i32)..(range.end as i32));
+
+    let mut i = 0;
+    let mut l = 0;
+    while l < len {
+        let chunk_size = chunk_sizes[i] as usize;
+        let chunk_end = std::cmp::min(l + chunk_size, len);
+        let chunk = &mut vals[l..chunk_end];
+
+        if saw_directions[i] == 0 {
+            chunk.sort_unstable();
+        } else if saw_directions[i] == 1 {
+            chunk.sort_unstable_by_key(|&e| std::cmp::Reverse(e));
+        } else {
+            unreachable!();
+        }
+
+        i += 1;
+        l += chunk_size;
+    }
+
+    vals
+}
+
+pub fn pipe_organ(len: usize) -> Vec<i32> {
+    //   .:.
+    // .:::::.
+
+    let mut vals = random_vec(len);
+
+    let first_half = &mut vals[0..(len / 2)];
+    first_half.sort_unstable();
+
+    let second_half = &mut vals[(len / 2)..len];
+    second_half.sort_unstable_by_key(|&e| std::cmp::Reverse(e));
+
+    vals
+}
+
+pub fn get_or_init_rand_seed() -> u64 {
+    *SEED_VALUE.get_or_init(|| {
+        env::var("OVERRIDE_SEED")
+            .ok()
+            .map(|seed| u64::from_str(&seed).unwrap())
+            .unwrap_or_else(rand_root_seed)
+    })
+}
+
+// --- Private ---
+
+static SEED_VALUE: OnceLock<u64> = OnceLock::new();
+
+#[cfg(not(miri))]
+fn rand_root_seed() -> u64 {
+    // Other test code hashes `panic::Location::caller()` and constructs a seed from that, in these
+    // tests we want to have a fuzzer like exploration of the test space, if we used the same caller
+    // based construction we would always test the same.
+    //
+    // Instead we use the seconds since UNIX epoch / 10, given CI log output this value should be
+    // reasonably easy to re-construct.
+
+    use std::time::{SystemTime, UNIX_EPOCH};
+
+    let epoch_seconds = SystemTime::now().duration_since(UNIX_EPOCH).unwrap().as_secs();
+
+    epoch_seconds / 10
+}
+
+#[cfg(miri)]
+fn rand_root_seed() -> u64 {
+    // Miri is usually run with isolation with gives us repeatability but also permutations based on
+    // other code that runs before.
+    use core::hash::{BuildHasher, Hash, Hasher};
+    let mut hasher = std::hash::RandomState::new().build_hasher();
+    core::panic::Location::caller().hash(&mut hasher);
+    hasher.finish()
+}
+
+fn random_vec(len: usize) -> Vec<i32> {
+    let mut rng: XorShiftRng = rand::SeedableRng::seed_from_u64(get_or_init_rand_seed());
+    (0..len).map(|_| rng.gen::<i32>()).collect()
+}
@@ -0,0 +1,208 @@
+// This module implements a Zipfian distribution generator.
+//
+// Based on https://github.com/jonhoo/rust-zipf.
+
+use rand::Rng;
+
+/// Random number generator that generates Zipf-distributed random numbers using rejection
+/// inversion.
+#[derive(Clone, Copy)]
+pub struct ZipfDistribution {
+    /// Number of elements
+    num_elements: f64,
+    /// Exponent parameter of the distribution
+    exponent: f64,
+    /// `hIntegral(1.5) - 1}`
+    h_integral_x1: f64,
+    /// `hIntegral(num_elements + 0.5)}`
+    h_integral_num_elements: f64,
+    /// `2 - hIntegralInverse(hIntegral(2.5) - h(2)}`
+    s: f64,
+}
+
+impl ZipfDistribution {
+    /// Creates a new [Zipf-distributed](https://en.wikipedia.org/wiki/Zipf's_law)
+    /// random number generator.
+    ///
+    /// Note that both the number of elements and the exponent must be greater than 0.
+    pub fn new(num_elements: usize, exponent: f64) -> Result<Self, ()> {
+        if num_elements == 0 {
+            return Err(());
+        }
+        if exponent <= 0f64 {
+            return Err(());
+        }
+
+        let z = ZipfDistribution {
+            num_elements: num_elements as f64,
+            exponent,
+            h_integral_x1: ZipfDistribution::h_integral(1.5, exponent) - 1f64,
+            h_integral_num_elements: ZipfDistribution::h_integral(
+                num_elements as f64 + 0.5,
+                exponent,
+            ),
+            s: 2f64
+                - ZipfDistribution::h_integral_inv(
+                    ZipfDistribution::h_integral(2.5, exponent)
+                        - ZipfDistribution::h(2f64, exponent),
+                    exponent,
+                ),
+        };
+
+        // populate cache
+
+        Ok(z)
+    }
+}
+
+impl ZipfDistribution {
+    fn next<R: Rng + ?Sized>(&self, rng: &mut R) -> usize {
+        // The paper describes an algorithm for exponents larger than 1 (Algorithm ZRI).
+        //
+        // The original method uses
+        //   H(x) = (v + x)^(1 - q) / (1 - q)
+        // as the integral of the hat function.
+        //
+        // This function is undefined for q = 1, which is the reason for the limitation of the
+        // exponent.
+        //
+        // If instead the integral function
+        //   H(x) = ((v + x)^(1 - q) - 1) / (1 - q)
+        // is used, for which a meaningful limit exists for q = 1, the method works for all
+        // positive exponents.
+        //
+        // The following implementation uses v = 0 and generates integral number in the range [1,
+        // num_elements]. This is different to the original method where v is defined to
+        // be positive and numbers are taken from [0, i_max]. This explains why the implementation
+        // looks slightly different.
+
+        let hnum = self.h_integral_num_elements;
+
+        loop {
+            use std::cmp;
+            let u: f64 = hnum + rng.gen::<f64>() * (self.h_integral_x1 - hnum);
+            // u is uniformly distributed in (h_integral_x1, h_integral_num_elements]
+
+            let x: f64 = ZipfDistribution::h_integral_inv(u, self.exponent);
+
+            // Limit k to the range [1, num_elements] if it would be outside
+            // due to numerical inaccuracies.
+            let k64 = x.max(1.0).min(self.num_elements);
+            // float -> integer rounds towards zero, so we add 0.5
+            // to prevent bias towards k == 1
+            let k = cmp::max(1, (k64 + 0.5) as usize);
+
+            // Here, the distribution of k is given by:
+            //
+            //   P(k = 1) = C * (hIntegral(1.5) - h_integral_x1) = C
+            //   P(k = m) = C * (hIntegral(m + 1/2) - hIntegral(m - 1/2)) for m >= 2
+            //
+            // where C = 1 / (h_integral_num_elements - h_integral_x1)
+            if k64 - x <= self.s
+                || u >= ZipfDistribution::h_integral(k64 + 0.5, self.exponent)
+                    - ZipfDistribution::h(k64, self.exponent)
+            {
+                // Case k = 1:
+                //
+                //   The right inequality is always true, because replacing k by 1 gives
+                //   u >= hIntegral(1.5) - h(1) = h_integral_x1 and u is taken from
+                //   (h_integral_x1, h_integral_num_elements].
+                //
+                //   Therefore, the acceptance rate for k = 1 is P(accepted | k = 1) = 1
+                //   and the probability that 1 is returned as random value is
+                //   P(k = 1 and accepted) = P(accepted | k = 1) * P(k = 1) = C = C / 1^exponent
+                //
+                // Case k >= 2:
+                //
+                //   The left inequality (k - x <= s) is just a short cut
+                //   to avoid the more expensive evaluation of the right inequality
+                //   (u >= hIntegral(k + 0.5) - h(k)) in many cases.
+                //
+                //   If the left inequality is true, the right inequality is also true:
+                //     Theorem 2 in the paper is valid for all positive exponents, because
+                //     the requirements h'(x) = -exponent/x^(exponent + 1) < 0 and
+                //     (-1/hInverse'(x))'' = (1+1/exponent) * x^(1/exponent-1) >= 0
+                //     are both fulfilled.
+                //     Therefore, f(x) = x - hIntegralInverse(hIntegral(x + 0.5) - h(x))
+                //     is a non-decreasing function. If k - x <= s holds,
+                //     k - x <= s + f(k) - f(2) is obviously also true which is equivalent to
+                //     -x <= -hIntegralInverse(hIntegral(k + 0.5) - h(k)),
+                //     -hIntegralInverse(u) <= -hIntegralInverse(hIntegral(k + 0.5) - h(k)),
+                //     and finally u >= hIntegral(k + 0.5) - h(k).
+                //
+                //   Hence, the right inequality determines the acceptance rate:
+                //   P(accepted | k = m) = h(m) / (hIntegrated(m+1/2) - hIntegrated(m-1/2))
+                //   The probability that m is returned is given by
+                //   P(k = m and accepted) = P(accepted | k = m) * P(k = m)
+                //                         = C * h(m) = C / m^exponent.
+                //
+                // In both cases the probabilities are proportional to the probability mass
+                // function of the Zipf distribution.
+
+                return k;
+            }
+        }
+    }
+}
+
+impl rand::distributions::Distribution<usize> for ZipfDistribution {
+    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> usize {
+        self.next(rng)
+    }
+}
+
+use std::fmt;
+impl fmt::Debug for ZipfDistribution {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> Result<(), fmt::Error> {
+        f.debug_struct("ZipfDistribution")
+            .field("e", &self.exponent)
+            .field("n", &self.num_elements)
+            .finish()
+    }
+}
+
+impl ZipfDistribution {
+    /// Computes `H(x)`, defined as
+    ///
+    ///  - `(x^(1 - exponent) - 1) / (1 - exponent)`, if `exponent != 1`
+    ///  - `log(x)`, if `exponent == 1`
+    ///
+    /// `H(x)` is an integral function of `h(x)`, the derivative of `H(x)` is `h(x)`.
+    fn h_integral(x: f64, exponent: f64) -> f64 {
+        let log_x = x.ln();
+        helper2((1f64 - exponent) * log_x) * log_x
+    }
+
+    /// Computes `h(x) = 1 / x^exponent`
+    fn h(x: f64, exponent: f64) -> f64 {
+        (-exponent * x.ln()).exp()
+    }
+
+    /// The inverse function of `H(x)`.
+    /// Returns the `y` for which `H(y) = x`.
+    fn h_integral_inv(x: f64, exponent: f64) -> f64 {
+        let mut t: f64 = x * (1f64 - exponent);
+        if t < -1f64 {
+            // Limit value to the range [-1, +inf).
+            // t could be smaller than -1 in some rare cases due to numerical errors.
+            t = -1f64;
+        }
+        (helper1(t) * x).exp()
+    }
+}
+
+/// Helper function that calculates `log(1 + x) / x`.
+/// A Taylor series expansion is used, if x is close to 0.
+fn helper1(x: f64) -> f64 {
+    if x.abs() > 1e-8 { x.ln_1p() / x } else { 1f64 - x * (0.5 - x * (1.0 / 3.0 - 0.25 * x)) }
+}
+
+/// Helper function to calculate `(exp(x) - 1) / x`.
+/// A Taylor series expansion is used, if x is close to 0.
+fn helper2(x: f64) -> f64 {
+    if x.abs() > 1e-8 {
+        x.exp_m1() / x
+    } else {
+        1f64 + x * 0.5 * (1f64 + x * 1.0 / 3.0 * (1f64 + 0.25 * x))
+    }
+}
@@ -1800,57 +1800,6 @@ fn brute_force_rotate_test_1() {
     }
 }
 
-#[test]
-#[cfg(not(target_arch = "wasm32"))]
-fn sort_unstable() {
-    use rand::Rng;
-
-    // Miri is too slow (but still need to `chain` to make the types match)
-    let lens = if cfg!(miri) { (2..20).chain(0..0) } else { (2..25).chain(500..510) };
-    let rounds = if cfg!(miri) { 1 } else { 100 };
-
-    let mut v = [0; 600];
-    let mut tmp = [0; 600];
-    let mut rng = crate::test_rng();
-
-    for len in lens {
-        let v = &mut v[0..len];
-        let tmp = &mut tmp[0..len];
-
-        for &modulus in &[5, 10, 100, 1000] {
-            for _ in 0..rounds {
-                for i in 0..len {
-                    v[i] = rng.gen::<i32>() % modulus;
-                }
-
-                // Sort in default order.
-                tmp.copy_from_slice(v);
-                tmp.sort_unstable();
-                assert!(tmp.windows(2).all(|w| w[0] <= w[1]));
-
-                // Sort in ascending order.
-                tmp.copy_from_slice(v);
-                tmp.sort_unstable_by(|a, b| a.cmp(b));
-                assert!(tmp.windows(2).all(|w| w[0] <= w[1]));
-
-                // Sort in descending order.
-                tmp.copy_from_slice(v);
-                tmp.sort_unstable_by(|a, b| b.cmp(a));
-                assert!(tmp.windows(2).all(|w| w[0] >= w[1]));
-            }
-        }
-    }
-
-    // Should not panic.
-    [0i32; 0].sort_unstable();
-    [(); 10].sort_unstable();
-    [(); 100].sort_unstable();
-
-    let mut v = [0xDEADBEEFu64];
-    v.sort_unstable();
-    assert!(v == [0xDEADBEEF]);
-}
-
 #[test]
 #[cfg(not(target_arch = "wasm32"))]
 #[cfg_attr(miri, ignore)] // Miri is too slow