比例拡張ソートを使用する

比例拡張ソート (proportion extend sort, PESort) は、整列済み接頭辞 S と未整列部分 U に分け、|U| > p·|S| の間は SU の先頭をまとめて整列して S を伸ばし、十分短い U だけが残ったら S の中央値で分割して再帰する。

クイックソートの分割を土台に、最悪計算量 O(n²) を避けようとする比較ソートである。

クイックソートは分割が 1 回の走査で済み、現代 CPU のメモリ階層と相性がよい。一方、ピボット選びが偏ると片側だけ巨大な部分問題が残り、最悪時間計算量が O(n²) になりうる。

比例拡張ソートは整列済み接頭辞 S と 未整列部分 U に区間を分けて考え、|U|p·|S| を超える間は SU の先頭 p·|S| 要素をまとめて整列して S を伸ばす。十分短い U だけが残ったら S の中央値をピボットに分割し、左右へ再帰する。サンプルソートに近い「小さな整列済み標本からピボットを得る」発想である。

  1. 接頭辞の初期化: 先頭 1 要素(または少数要素)を整列済み接頭辞 S とみなす。
  2. 比例拡張: |U| > p·|S| のあいだ、SU の先頭 p·|S| 要素をまとめて整列し、結果を新しい S とする。
  3. 中央値ピボット: S の中央値をピボットとして部分配列を分割する(実装では Lomuto 型などクイックソートと同様の分割が用いられる)。
  4. 再帰: ピボットより小さい側・大きい側に同じ手順を適用する。十分短い区間は挿入ソートで仕上げる。
procedure proportion_extend_sort(A, lo, hi, p)
  if hi - lo <= INSERTION_THRESHOLD then
    insertion_sort(A, lo, hi)
    return
  s_end = lo
  while s_end < hi and (hi - s_end) > p * (s_end - lo + 1)
    chunk_end = min(s_end + p * (s_end - lo + 1), hi)
    proportion_extend_sort(A, lo, chunk_end, p)
    s_end = chunk_end
  median = lo + floor((s_end - lo) / 2)
  pivot = partition(A, lo, hi, median)
  proportion_extend_sort(A, lo, pivot - 1, p)
  proportion_extend_sort(A, pivot + 1, hi, p)

パラメータ p は「未整列部分 ÷ 整列済み接頭辞」の上限である(文献によって定義が異なるので実装時は論文・ソースの取り方に注意する)。p ≈ 16 付近で平均性能がよく、p を小さくすると最悪ケースが改善しやすい。p が無限大に近いとクイックソート、未整列部分が 1 要素に近いと二分探索挿入ソートに近い振る舞いになる。

クイックソートの最悪計算量 O(n²) を避けつつ、最悪計算量 O(n log n) を保つが、一般に不安定である。

以下のデモでは小配列でも拡張フェーズが見えるように p = 2、小区間閾値 4、分割は Lomuto 型で実行する。青枠は整列済み接頭辞、紫は中央値ピボット、オレンジは比較、緑は交換を表す。

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

クイックソートと同様に分割するが、整列済み接頭辞を伸ばしてから中央値で切るため最悪計算量 O(n log n) を狙える。対称分割ソートは標本を両端へ置いてから分割する。

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

Size Average time Maximum time Average memory Maximum memory
256 0.000006 0.000390 1661 1668
512 0.000019 0.000474 1666 1672
1024 0.000045 0.000513 1673 1680
2048 0.000102 0.000424 1690 1696
4096 0.000218 0.000450 1722 1728
8192 0.000609 0.001594 1786 1792
16384 0.001453 0.002719 1918 1924
32768 0.003233 0.005752 2178 2184
65536 0.006769 0.013338 2690 2696
131072 0.019993 0.134023 3713 3720
262144 0.046617 0.065233 5762 5768
計測に使用したコードを表示する

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::{
    env,
    process::Command,
    time::{Duration, Instant},
};
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 pe_sort_range(a: &mut [usize], lo: usize, hi: usize) {
    const P: usize = 16;
    if hi <= lo {
        return;
    }
    if hi - lo < 16 {
        insertion_sort(&mut a[lo..=hi]);
        return;
    }
    let mut s_end = lo;
    while s_end < hi && hi - s_end > P * (s_end - lo + 1) {
        let chunk_end = (s_end + P * (s_end - lo + 1)).min(hi);
        pe_sort_range(a, lo, chunk_end);
        s_end = chunk_end;
    }
    let median = lo + (s_end - lo) / 2;
    let pivot = partition_at(a, lo, hi, median);
    if pivot > 0 {
        pe_sort_range(a, lo, pivot - 1);
    }
    pe_sort_range(a, pivot + 1, hi);
}

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


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

    proportion_extend_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],
    );
    for seed in 0..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 memory_usage_kb() -> usize {
    // VmHWM (peak RSS, KiB). Reported memory subtracts a per-size baseline that only
    // holds the input array, so the table reflects auxiliary space during sorting.
    let contents = std::fs::read_to_string("/proc/self/status")
        .unwrap_or_default();

    for line in contents.lines() {
        if let Some(rest) = line.strip_prefix("VmHWM:") {
            let kb = rest
                .split_whitespace()
                .next()
                .unwrap_or("0")
                .parse::<usize>()
                .unwrap_or(0);

            return kb;
        }
    }

    0
}

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

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

fn run_baseline(size: usize) -> usize {
    let _hold = input_array(size, 1);
    memory_usage_kb()
}

fn run_once(size: usize, seed: usize) -> (u128, usize) {
    let mut array = input_array(size, seed as u64);

    let start = Instant::now();

    benchmark_sort(&mut array);

    let elapsed = start.elapsed();
    let mem = memory_usage_kb();

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

    (micros(elapsed), mem)
}

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

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 == "--baseline-once") {
        run_baseline_child(&args);
        return;
    }
    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 baseline_output = Command::new(env::current_exe().expect("failed to find current executable"))
            .arg("--baseline-once")
            .arg(size.to_string())
            .output()
            .expect("failed to run benchmark baseline process");

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

        let baseline_stdout = String::from_utf8(baseline_output.stdout)
            .expect("baseline process returned non-UTF-8 output");
        let baseline_mem = baseline_stdout
            .split_whitespace()
            .next()
            .expect("missing baseline memory usage")
            .parse::<usize>()
            .expect("invalid baseline memory usage");

        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 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;
            }

            let aux_mem = mem.saturating_sub(baseline_mem);

            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