対称分割ソートを使用する

対称分割ソート (symmetry partition sort, SPSort) は、比例拡張ソートと同様に接頭辞を伸ばしてから中央値で分割するが、整列済み標本を配列の両端(L | U | R)へ置いてから未整列部分 U を分割する。

クイックソートの最悪計算量 O(n²) を避けつつ分割処理を速くする比較ソートである。

PESort は整列済み接頭辞 S を左端に伸ばし、その中央値をピボットに未整列部分 U を分割する。

SPSort も同じ比例拡張の枠組みを使うが、分割に入る直前に S を中央値で左右に分け、左半分 L を左端、右半分 R を右端に置き、中央に U を残す L | U | R の形に整える。 LR はそれぞれピボット未満・以上の要素群として働き、分割ループの番兵(Sentinel)になるため、境界チェックを減らしつつ走査を速くできる。

  1. 接頭辞の初期化: 先頭 1 要素(または少数要素)を整列済み接頭辞 S とみなす。
  2. 比例拡張: |U| > p·|S| のあいだ、SU の先頭 p·|S| 要素をまとめて整列し、結果を新しい S とする。
  3. 対称配置: S の中央値位置で左右に分け、右半分 R を配列右端へ移動して L | U | R にする。
  4. 中央値ピボット: 中央値をピボットとして U を分割する(実装ではロムート型や三分割型などが用いられる)。
  5. 再帰: 左側(L とピボット未満の部分)・右側(ピボット以上と R)に同じ手順を適用する。十分短い区間は挿入ソートで仕上げる。
procedure symmetry_partition_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)
    symmetry_partition_sort(A, lo, chunk_end, p)
    s_end = chunk_end
  median = lo + floor((s_end - lo) / 2)
  r_len = s_end - median
  move A[median + 1 .. s_end] to A[hi - r_len + 1 .. hi]   // L | U | R
  pivot = partition(A, lo, hi, median)
  symmetry_partition_sort(A, lo, pivot - 1, p)
  symmetry_partition_sort(A, pivot + 1, hi, p)

パラメータ p は PESort と同じく「未整列部分 ÷ 整列済み接頭辞」の上限である。p ≈ 16 付近で平均性能がよく、小さくすると最悪ケースが改善しやすい。

Jing-Chao Chen の実装では 適応的な連続区間の検出等値キー領域を持つ分割 も組み合わされ、部分整列済み入力への適応性が高い。

PESort と同様に最悪計算量 O(n log n) を保ちつつ、対称配置で分割を速くするが、一般に不安定である。

以下のデモでは p = 2(小配列でも拡張・対称配置が見えるよう意図的に小さくしている)、小区間閾値 4、分割はロムート型である。

青枠は左端の整列済み標本 L、水色は右端の R、紫は中央値ピボット、オレンジは比較、緑は交換を表す。

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

比例拡張ソートと同系統だが、整列済み標本を L | U | R と両端へ配置してから未整列部分を分割していて、分割の局所性を改善しやすい。

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

Size Average time Maximum time Average memory Maximum memory
256 0.000006 0.000040 0 0
512 0.000018 0.000060 0 0
1024 0.000041 0.000126 0 0
2048 0.000092 0.000167 0 0
4096 0.000195 0.000369 0 0
8192 0.000567 0.000883 0 0
16384 0.001305 0.002071 0 0
32768 0.002827 0.004113 0 0
65536 0.005927 0.010777 0 0
131072 0.018084 0.033297 0 0
262144 0.040957 0.079740 0 0
計測に使用したコードを表示する

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 sp_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);
        sp_sort_range(a, lo, chunk_end);
        s_end = chunk_end;
    }
    let median = lo + (s_end - lo) / 2;
    let r_len = s_end - median;
    let r_start = median + 1;
    let r_dest = hi - r_len + 1;
    for i in 0..r_len {
        let from = r_start + i;
        let to = r_dest + i;
        if from != to {
            a.swap(from, to);
        }
    }
    let pivot = partition_at(a, lo, hi, median);
    if pivot > 0 {
        sp_sort_range(a, lo, pivot - 1);
    }
    sp_sort_range(a, pivot + 1, hi);
}

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


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

    symmetry_partition_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