シバーズソートを使用する

シバーズソート (shivers sort) は、入力にすでに存在する単調な連続区間(ラン)を検出し、ラン長の「レベル」に基づいてスタック最上段 2 本をいつ併合するかを決める自然マージソートである。

ティムソートより単純な、スタック上端 2 本だけを見て決めるマージ方針でありながら、入力長 n だけを見たときの最悪マージコストが n log₂(n) + O(n) に収まることが後に示された。

  1. ランの検出: 左から昇順または厳密な降順の連続区間を見つける。降順ランは反転して昇順にそろえる。
  2. ランの拡張: 長さが最小ラン長 min_run 未満なら挿入ソートで伸ばす(デモでは見やすさのため 4 に固定。計測実装では 32)。
  3. レベルの計算: ラン長 r に対し ℓ = ⌊log₂(r)⌋(パラメータ c = 1)をレベルとする。
  4. スタックに従ったマージ: ランを左から 1 本ずつ積み、スタック高さが 2 以上かつ ℓ_h ≥ ℓ_{h-1} なら最上段 2 本(R_{h-1}R_h)を併合する。
  5. 仕上げ: 入力のランをすべて積み終わったあと、スタックに 2 本以上残っていれば上から順に 2 本ずつ併合する。

条件 ℓ_h ≥ ℓ_{h-1}2^{⌊log₂|Y|⌋} ≤ |Z|Y = R_{h-1}Z = R_h)と同値である。長さそのものではなく 2 の冪へ丸めたレベルで比較するため、スタック上のレベル列をほぼ狭義減少に保ちやすい。

一方で、新しく積んだランがスタック上の既存ランよりはるかに長くても即座に併合しうる。そのためラン数 ρ やラン長分布への適応性は弱く、最悪マージコストは ω(n log₂(ρ)) になりうる(n に対しては最適に近い)。

procedure level(r)  // c = 1
  return floor(log2(r))

procedure shivers_sort(A)
  runs := run_decomposition(A)  // 反転・min_run 拡張済み
  S := empty stack
  while true
    h := height(S)
    if h >= 2 and level(len(R_h)) >= level(len(R_{h-1}))
      merge R_{h-1} and R_h on S
    else if runs is not empty
      push next run from runs onto S
    else
      break
  while height(S) >= 2
    merge R_{h-1} and R_h on S

安定ソートであり、最悪でもマージソートと同程度の O(n log n) を保つ。ただし整列済みに近い入力でも、適応型シバーズソートやパワーソートほどラン長エントロピーに追従する保証はない。

実装の核は「長さをレベルに丸めてから最上段 2 本を比べる」だけであり、ティムソートの 3 本比較やパワーソートの中点ベースのパワー計算より短い。その単純さの代償として、ラン構造への適応性は後続の改良版に譲る。

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

ティムソートはスタック上端のラン長そのものでマージを決める。パワーソートは隣接ランの中点からパワーを求め、ほぼ最適な二分マージ木に沿う。

シバーズソートはレベル比較を使う点で適応型シバーズソートと同じだが、マージの判定も併合もスタック上端 2 本(R_{h-1}R_h)だけに限る方針である。

適応型シバーズソートはティムソート寄りに R_{h-2}R_{h-1} をマージして最上段を残し、マージコストをラン長エントロピーに対して nH + O(n) へ近づけた。

ナチュラルマージソートは自然ランを使うが、マージ順は単純なペア併合にとどまる。

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

Size Average time Maximum time Average memory Maximum memory
256 0.000008 0.000039 2 2
512 0.000017 0.000060 4 4
1024 0.000043 0.000179 8 8
2048 0.000089 0.000228 17 17
4096 0.000173 0.000524 34 34
8192 0.000379 0.000745 68 68
16384 0.000730 0.001922 136 136
32768 0.001806 0.003574 272 272
65536 0.003374 0.008681 544 544
131072 0.006674 0.015976 1088 1088
262144 0.014530 0.034874 2176 2176
計測に使用したコードを表示する

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 as AtomicOrdering},
    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, AtomicOrdering::Relaxed) + size;
    PEAK_BYTES.fetch_max(live, AtomicOrdering::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(), AtomicOrdering::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(), AtomicOrdering::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;
        }
    }
}



#[derive(Clone, Copy)]
struct ShiversRun {
    lo: usize,
    hi: usize,
}

fn run_level(len: usize) -> u32 {
    debug_assert!(len > 0);
    (usize::BITS - 1 - len.leading_zeros()) as u32
}

fn merge_shivers_runs(a: &mut [usize], left: ShiversRun, right: ShiversRun) -> ShiversRun {
    let lo = left.lo;
    let hi = right.hi;
    let mid = left.hi + 1;
    let mut merged = Vec::with_capacity(hi - lo + 1);
    let (mut l, mut r) = (left.lo, mid);
    while l <= left.hi && r <= right.hi {
        if a[l] <= a[r] {
            merged.push(a[l]);
            l += 1;
        } else {
            merged.push(a[r]);
            r += 1;
        }
    }
    merged.extend_from_slice(&a[l..=left.hi]);
    merged.extend_from_slice(&a[r..=right.hi]);
    a[lo..=hi].copy_from_slice(&merged);
    ShiversRun { lo, hi }
}

fn prepare_shivers_run(a: &mut [usize], start: usize, min_run: usize) -> usize {
    let n = a.len();
    let mut i = start + 1;
    if i < n && a[i - 1] > a[i] {
        while i < n && a[i - 1] > a[i] {
            i += 1;
        }
        a[start..i].reverse();
    } else {
        while i < n && a[i - 1] <= a[i] {
            i += 1;
        }
    }
    let end = (start + min_run).min(n).max(i);
    insertion_sort(&mut a[start..end]);
    end
}

fn shivers_sort(a: &mut [usize]) {
    const MIN_RUN: usize = 32;
    let n = a.len();
    if n <= 1 {
        return;
    }

    let mut pending: Vec<ShiversRun> = Vec::new();
    let mut start = 0usize;
    while start < n {
        let end = prepare_shivers_run(a, start, MIN_RUN);
        pending.push(ShiversRun {
            lo: start,
            hi: end - 1,
        });
        start = end;
    }

    let mut stack: Vec<ShiversRun> = Vec::new();
    let mut next = 0usize;
    loop {
        let h = stack.len();
        if h >= 2 {
            let r_hm1 = stack[h - 2];
            let r_h = stack[h - 1];
            let ell_hm1 = run_level(r_hm1.hi - r_hm1.lo + 1);
            let ell_h = run_level(r_h.hi - r_h.lo + 1);
            if ell_h >= ell_hm1 {
                let right = stack.pop().unwrap();
                let left = stack.pop().unwrap();
                stack.push(merge_shivers_runs(a, left, right));
                continue;
            }
        }
        if next < pending.len() {
            stack.push(pending[next]);
            next += 1;
            continue;
        }
        break;
    }

    while stack.len() >= 2 {
        let right = stack.pop().unwrap();
        let left = stack.pop().unwrap();
        stack.push(merge_shivers_runs(a, left, right));
    }
}


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

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

// Skip cases larger than the algorithm's measured size cap (MAX_POWER). That
// cap exists because larger inputs are impractically slow; forcing them here
// would stall the published measurement script before any table rows print.
fn check_correctness_case_within_limit(label: &str, input: Vec<usize>) {
    if input.len() > (1usize << MAX_POWER) {
        return;
    }
    check_correctness_case(label, input);
}

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),
        );
    }
    // Small-input cutoffs (insertion sort below 32 elements, etc.) hide duplicate-key
    // bugs in the recursive path, so repeat the duplicate cases at the smallest
    // benchmark size, which every algorithm must handle within reasonable time.
    check_correctness_case("all_equal_len256", vec![7; 256]);
    for seed in 1..=4 {
        check_correctness_case(
            &format!("few_keys_len256_seed_{seed}"),
            few_unique_values(256, 4, seed),
        );
    }
    // Blit's equal-key second sweep used to copy the whole range into a fixed
    // 512-element swap; lengths above that must still sort without panicking.
    // Respect MAX_POWER so algorithms with a low measured-size cap (slow,
    // sleep) do not hang here for minutes or months.
    check_correctness_case_within_limit("all_equal_len600", vec![7; 600]);
    for seed in 1..=4 {
        check_correctness_case_within_limit(
            &format!("few_keys_len2048_seed_{seed}"),
            few_unique_values(2048, 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 bytes (explicit buffers such as swap).
/// Kept in bytes so the parent can average before rounding; converting to KiB here
/// would truncate sub-KiB buffers to 0 in every run and hide them from the average.
fn run_once(size: usize, seed: usize) -> (u128, usize) {
    let mut array = input_array(size, seed as u64);

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

    let start = Instant::now();

    benchmark_sort(&mut array);

    let elapsed = start.elapsed();
    let peak_bytes = PEAK_BYTES.load(AtomicOrdering::Relaxed);
    let aux_bytes = peak_bytes.saturating_sub(base_bytes);

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

    (micros(elapsed), aux_bytes)
}

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;
        // Memory is summed in bytes and converted to KiB once, after averaging.
        let avg_mem_kb = total_mem / RUNS / 1024;
        let max_mem_kb = max_mem / 1024;

        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_kb,
            max_mem_kb
        );
    }
}
RUST

RUN cargo build --release

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

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