サンプルソートを使用する

サンプルソート (sample sort) は、入力から少数の標本を取り出して整列し、その値を分割点(スプリッター)として全体をバケットに分け、各バケットを整列して連結する。

並列計算向けに設計された手法として知られ、プロセッサ数 p に対して p−1 個のスプリッターを選ぶ構成が典型である。

クイックソートが 1 つのピボットで左右に分けるのに対し、サンプルソートは複数のスプリッターで p 個(デモでは要素数に応じた複数個)のバケットに一度に仕分ける。分割の偏りを抑えやすく、各バケットを別プロセッサで同時に整列できる点が並列版の利点である。

  1. サンプリング: 配列から s 個の要素を(等間隔に)選び、サンプル集合を得る。文献では s ≈ √ns = p−1 などが用いられる。
  2. サンプルの整列: サンプルだけを整列し、昇順のスプリッター列 t₀ ≤ t₁ ≤ … を得る。
  3. 仕分け: 各要素 x を、t₀, t₁, … と比較して属するバケット番号を決める(例: x ≤ t₀ ならバケット 0、t₀ < x ≤ t₁ ならバケット 1、…)。
  4. バケット整列: 各バケットを独立に整列する(再帰的にサンプルソート、または十分小さければ挿入ソートなど)。
  5. 連結: バケット 0, 1, … の順に並べれば全体が昇順になる。
procedure sample_sort(A, s)
  if length(A) <= SMALL then
    insertion_sort(A)
    return
  S = choose_s_samples(A, s)
  sort(S)
  splitters = S
  buckets = empty list of s + 1 arrays
  for each x in A
    b = bucket_index(x, splitters)
    append x to buckets[b]
  for each bucket B in buckets
    sample_sort(B, s)
  A = concatenate(buckets)

プロセッサ p 台では O(n/p) 程度の並列時間を狙える。

比例拡張ソートのように「小さな整列済み標本から分割点を得る」発想はサンプルソートと共通する。バブルソートのように隣接交換だけで進む単純ソートや、クイックソートの単一ピボット再帰と比べると、標本に基づく多分割と並列化しやすい区切りが特徴的である。

以下のデモでは 15 要素から 4 個のサンプルを取り、スプリッター 4 個で 5 バケットに仕分けたあと、各バケットを挿入ソートで仕上げる。

外部ソートや分散整列の文脈では I/O 効率のよい分割としても使われる。サンプル数や再帰の打ち切り、等値要素の扱いが性能を左右する。

類似アルゴリズムとの相違点

ロムート分割型クイックソートが 1 つのピボットで 2 分割するのに対し、標本から複数スプリッターを求めて同時に多分割する。並列・外部整列向きである。バケットソートは値域を等分する点が異なる。

計算時間量および空間計算量を計測する

Size Average time Maximum time Average memory Maximum memory
256 0.000007 0.000047 3 6
512 0.000015 0.000069 7 10
1024 0.000029 0.000071 14 19
2048 0.000057 0.000129 27 32
4096 0.000108 0.000281 52 62
8192 0.000225 0.000655 101 120
16384 0.000446 0.003172 197 217
32768 0.000904 0.001899 388 429
65536 0.001820 0.003310 766 819
131072 0.003751 0.005948 1517 1612
262144 0.007610 0.014622 3012 3175
計測に使用したコードを表示する

set -euo pipefail

WORKDIR="$(mktemp -d)"
trap 'rm -rf "$WORKDIR"' EXIT

cat > "$WORKDIR/Dockerfile" <<'EOF'
FROM rust:1.95.0

WORKDIR /app

RUN mkdir -p src

RUN cat > Cargo.toml <<'CARGO'
[package]
name = "rust-benchmark"
version = "0.1.0"
edition = "2021"

[profile.release]
lto = true
codegen-units = 1
panic = "abort"
CARGO

RUN cat > src/main.rs <<'RUST'
use std::{
    alloc::{GlobalAlloc, Layout, System},
    env,
    process::Command,
    sync::atomic::{AtomicUsize, Ordering},
    time::{Duration, Instant},
};

/// Counts live heap bytes and the high-water mark so auxiliary sort buffers
/// (swap Vecs, etc.) are measured as explicit heap growth during the sort.
struct TrackingAllocator;

static LIVE_BYTES: AtomicUsize = AtomicUsize::new(0);
static PEAK_BYTES: AtomicUsize = AtomicUsize::new(0);

fn record_alloc(size: usize) {
    let live = LIVE_BYTES.fetch_add(size, Ordering::Relaxed) + size;
    PEAK_BYTES.fetch_max(live, Ordering::Relaxed);
}

unsafe impl GlobalAlloc for TrackingAllocator {
    unsafe fn alloc(&self, layout: Layout) -> *mut u8 {
        let ptr = System.alloc(layout);
        if !ptr.is_null() {
            record_alloc(layout.size());
        }
        ptr
    }

    unsafe fn dealloc(&self, ptr: *mut u8, layout: Layout) {
        LIVE_BYTES.fetch_sub(layout.size(), Ordering::Relaxed);
        System.dealloc(ptr, layout);
    }

    unsafe fn alloc_zeroed(&self, layout: Layout) -> *mut u8 {
        let ptr = System.alloc_zeroed(layout);
        if !ptr.is_null() {
            record_alloc(layout.size());
        }
        ptr
    }

    unsafe fn realloc(&self, ptr: *mut u8, layout: Layout, new_size: usize) -> *mut u8 {
        let new_ptr = System.realloc(ptr, layout, new_size);
        if !new_ptr.is_null() {
            LIVE_BYTES.fetch_sub(layout.size(), Ordering::Relaxed);
            record_alloc(new_size);
        }
        new_ptr
    }
}

#[global_allocator]
static GLOBAL: TrackingAllocator = TrackingAllocator;
const MIN_POWER: u32 = 8;
const MAX_POWER: u32 = 18;
const RUNS: usize = 8192;
fn insertion_sort(a: &mut [usize]) {
    for i in 1..a.len() {
        let mut j = i;
        while j > 0 && a[j - 1] > a[j] {
            a.swap(j - 1, j);
            j -= 1;
        }
    }
}

fn partition_at(a: &mut [usize], lo: usize, hi: usize, pivot_idx: usize) -> usize {
    a.swap(pivot_idx, hi);
    let pivot = a[hi];
    let mut i = lo;
    for j in lo..hi {
        if a[j] < pivot {
            a.swap(i, j);
            i += 1;
        }
    }
    a.swap(i, hi);
    i
}

fn partition(a: &mut [usize], lo: usize, hi: usize) -> usize {
    partition_at(a, lo, hi, lo + (hi - lo) / 2)
}

fn quick_sort_range(a: &mut [usize], lo: usize, hi: usize) {
    if hi <= lo {
        return;
    }
    if hi - lo < 16 {
        insertion_sort(&mut a[lo..=hi]);
        return;
    }
    let p = partition(a, lo, hi);
    if p > 0 {
        quick_sort_range(a, lo, p - 1);
    }
    quick_sort_range(a, p + 1, hi);
}

fn quick_sort(a: &mut [usize]) {
    if let Some(hi) = a.len().checked_sub(1) {
        quick_sort_range(a, 0, hi);
    }
}



fn sample_sort(a: &mut [usize]) {
    if a.len() <= 32 {
        insertion_sort(a);
        return;
    }
    let sample_count = (a.len() as f64).sqrt() as usize;
    let step = (a.len() / sample_count.max(1)).max(1);
    let mut splitters: Vec<usize> = (step - 1..a.len())
        .step_by(step)
        .take(sample_count)
        .map(|i| a[i])
        .collect();
    quick_sort(&mut splitters);
    let mut buckets = vec![Vec::new(); splitters.len() + 1];
    for &value in a.iter() {
        let bucket = splitters.partition_point(|&splitter| value > splitter);
        buckets[bucket].push(value);
    }
    let mut pos = 0;
    for bucket in buckets.iter_mut() {
        sample_sort(bucket);
        for &value in bucket.iter() {
            a[pos] = value;
            pos += 1;
        }
    }
}


fn benchmark_sort(array: &mut [usize]) {

    sample_sort(array);

}

fn is_non_decreasing(a: &[usize]) -> bool {
    a.windows(2).all(|w| w[0] <= w[1])
}

fn same_multiset(a: &[usize], b: &[usize]) -> bool {
    if a.len() != b.len() {
        return false;
    }

    let mut left = a.to_vec();
    let mut right = b.to_vec();
    left.sort_unstable();
    right.sort_unstable();
    left == right
}

fn check_correctness_case(label: &str, mut input: Vec<usize>) {
    let original = input.clone();

    benchmark_sort(&mut input);

    if !is_non_decreasing(&input) {
        panic!("correctness case {}: output is not sorted", label);
    }

    if !same_multiset(&input, &original) {
        panic!("correctness case {}: elements were lost or added", label);
    }
}

fn few_unique_values(size: usize, unique: usize, seed: u64) -> Vec<usize> {
    let mut state = seed;

    (0..size)
        .map(|_| {
            state ^= state << 13;
            state ^= state >> 7;
            state ^= state << 17;
            (state as usize % unique) + 1
        })
        .collect()
}

fn run_correctness_checks() {
    check_correctness_case("empty", vec![]);
    check_correctness_case("single", vec![42]);
    check_correctness_case("duplicates", vec![3, 1, 3, 2, 1, 2]);
    check_correctness_case("sorted", vec![1, 2, 3, 4, 5]);
    check_correctness_case("reverse", vec![5, 4, 3, 2, 1]);
    check_correctness_case("all_equal", vec![7, 7, 7, 7]);
    check_correctness_case("skewed_range", vec![1_000_000, 2, 1_000_001, 1, 999_999]);
    // Static-buffer Grail skips the in-buffer build when key collection is sparse
    // (ideal_buffer = false). Exercising that path catches regressions in buffer gating.
    check_correctness_case(
        "few_keys_len16",
        vec![2, 2, 2, 2, 2, 2, 2, 2, 4, 3, 1, 2, 3, 4, 1, 4],
    );
    // Seed 0 is a fixed point of the xorshift below, so it would degenerate into
    // yet another all-equal case instead of a 4-value mix. Start at 1.
    for seed in 1..=32 {
        check_correctness_case(
            &format!("few_keys_len32_seed_{seed}"),
            few_unique_values(32, 4, seed),
        );
    }
}


fn shuffled(size: usize, seed: u64) -> Vec<usize> {
    let mut v: Vec<usize> = (1..=size).collect();

    let mut state = seed;

    for i in (1..size).rev() {
        state ^= state << 13;
        state ^= state >> 7;
        state ^= state << 17;

        let j = (state as usize) % (i + 1);

        v.swap(i, j);
    }

    v
}

fn micros(d: Duration) -> u128 {
    d.as_micros()
}

fn input_array(size: usize, seed: u64) -> Vec<usize> {
    shuffled(size, seed)
}

/// Peak heap growth during `benchmark_sort`, in KiB (explicit buffers such as swap).
fn run_once(size: usize, seed: usize) -> (u128, usize) {
    let mut array = input_array(size, seed as u64);

    let base_bytes = LIVE_BYTES.load(Ordering::Relaxed);
    PEAK_BYTES.store(base_bytes, Ordering::Relaxed);

    let start = Instant::now();

    benchmark_sort(&mut array);

    let elapsed = start.elapsed();
    let peak_bytes = PEAK_BYTES.load(Ordering::Relaxed);
    let aux_kb = peak_bytes.saturating_sub(base_bytes) / 1024;

    let expected: Vec<usize> = (1..=size).collect();
    if array != expected {
        panic!(
            "sort failed with seed {} for size {}",
            seed,
            size
        );
    }

    (micros(elapsed), aux_kb)
}

fn run_child(args: &[String]) {
    let size = args[2].parse::<usize>().expect("invalid size");
    let seed = args[3].parse::<usize>().expect("invalid seed");
    let (elapsed_us, mem) = run_once(size, seed);
    println!("{} {}", elapsed_us, mem);
}

fn main() {
    let args: Vec<String> = env::args().collect();
    if args.get(1).is_some_and(|arg| arg == "--run-once") {
        run_child(&args);
        return;
    }

    run_correctness_checks();

    println!(
        "| {:>10} | {:>15} | {:>15} | {:>15} | {:>15} |",
        "Size",
        "Average time",
        "Maximum time",
        "Average memory",
        "Maximum memory"
    );

    println!(
        "|{:-<11}:|{:-<16}:|{:-<16}:|{:-<16}:|{:-<16}:|",
        "",
        "",
        "",
        "",
        ""
    );

    for power in MIN_POWER..=MAX_POWER {
        let size = 1usize << power;

        let mut total_time: u128 = 0;
        let mut max_time: u128 = 0;

        let mut total_mem: usize = 0;
        let mut max_mem: usize = 0;

        for seed in 1..=RUNS {
            let output = Command::new(env::current_exe().expect("failed to find current executable"))
                .arg("--run-once")
                .arg(size.to_string())
                .arg(seed.to_string())
                .output()
                .expect("failed to run benchmark child process");

            if !output.status.success() {
                panic!(
                    "benchmark child process failed: {}",
                    String::from_utf8_lossy(&output.stderr)
                );
            }

            let stdout = String::from_utf8(output.stdout)
                .expect("child process returned non-UTF-8 output");
            let mut fields = stdout.split_whitespace();
            let elapsed_us = fields
                .next()
                .expect("missing elapsed time")
                .parse::<u128>()
                .expect("invalid elapsed time");
            let aux_mem = fields
                .next()
                .expect("missing memory usage")
                .parse::<usize>()
                .expect("invalid memory usage");

            total_time += elapsed_us;

            if elapsed_us > max_time {
                max_time = elapsed_us;
            }

            total_mem += aux_mem;

            if aux_mem > max_mem {
                max_mem = aux_mem;
            }
        }

        let avg_time = total_time / RUNS as u128;
        let avg_mem = total_mem / RUNS;

        println!(
            "| {:>10} | {:>15} | {:>15} | {:>15} | {:>15} |",
            size,
            format!("{}.{:06}", avg_time / 1_000_000, avg_time % 1_000_000),
            format!("{}.{:06}", max_time / 1_000_000, max_time % 1_000_000),
            avg_mem,
            max_mem
        );
    }
}
RUST

RUN cargo build --release

CMD ["./target/release/rust-benchmark"]
EOF

docker build -t rust-benchmark "$WORKDIR"
docker run --rm --init rust-benchmark