適応型シバーズソートを使用する

適応型シバーズソート (adaptive shivers sort) は、入力にすでに存在する単調な連続区間(ラン)を検出し、ラン長の「レベル」に基づいてマージ順を決める自然マージソートである。

ティムソートと同じく、スタック上端 3 本のラン長だけを見てマージを決める方針を保ちつつ、シバーズソート由来のレベル比較で最悪時のマージコストを抑えやすくした折衷案と捉えられる。マージコストがラン長エントロピーに対して nH + O(n) に収まることが示されている。

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

ティムソートはラン長そのものの大小関係でマージを決めるのに対し、適応型シバーズソートは長さを 2 の冪の段階(レベル)に丸めてから比較する。その結果、同程度の長さのランを早めにまとめる方針が明示的になり、最悪ケースの定数が改善しやすい。

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

procedure adaptive_shivers_sort(A)
  runs := run_decomposition(A)  // 反転・min_run 拡張済み
  S := empty stack
  while true
    h := height(S)
    if h >= 3 and level(len(R_{h-2})) <= max(level(len(R_{h-1})), level(len(R_h)))
      merge R_{h-2} and R_{h-1} 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) に近づく。最悪でもマージソートと同程度の O(n log n) を保つ。

ティムソートからの差し替えはマージ抑制条件をレベル比較に変える程度で済む一方、パワーソートのような中点ベースのパワー計算は不要で、実装と証明の両方が短い点が利点である。

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

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

適応型シバーズソートはティムソート同様、スタック上端 3 本だけを見てマージを決めるが、長さを ⌊log₂(r)⌋ に丸めてから比較する。

古典的なシバーズソートは条件成立時に最上段 2 本をマージするのに対し、適応型はティムソート寄りに R_{h-2}R_{h-1} をマージして最上段を残す。その差が最悪マージコストの改善につながる。

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

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

Size Average time Maximum time Average memory Maximum memory
256 0.000008 0.000065 2 2
512 0.000018 0.000072 4 4
1024 0.000039 0.000394 8 8
2048 0.000084 0.000160 17 17
4096 0.000177 0.000272 34 34
8192 0.000379 0.000864 68 68
16384 0.000857 0.001699 136 136
32768 0.001624 0.013000 272 272
65536 0.003461 0.006211 544 544
131072 0.007483 0.020498 1088 1088
262144 0.015821 0.032127 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 adaptive_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 >= 3 {
            let r_hm2 = stack[h - 3];
            let r_hm1 = stack[h - 2];
            let r_h = stack[h - 1];
            let ell_hm2 = run_level(r_hm2.hi - r_hm2.lo + 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_hm2 <= ell_hm1.max(ell_h) {
                let top = stack.pop().unwrap();
                let mid = stack.pop().unwrap();
                let left = stack.pop().unwrap();
                let merged = merge_shivers_runs(a, left, mid);
                stack.push(merged);
                stack.push(top);
                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]) {

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

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.
    check_correctness_case("all_equal_len600", vec![7; 600]);
    for seed in 1..=4 {
        check_correctness_case(
            &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