置換選択ソートを使用する

置換選択ソート (replacement selection sort) は、限られた大きさの最小ヒープで入力を流し、平均でヒープ容量の約 2 倍の長さの整列済みランを生成し、それらを併合して全体を昇順にする整列である。

外部整列ではメモリに載らないファイルを扱うとき、単純にメモリ分だけ読んでクイックソートすると初期ラン長はメモリ容量 M に留まる。 置換選択では、出力した直前のキー以上の入力だけをヒープへ戻す(置換する)ことで、ランダム入力でも期待ラン長がおよそ 2M になる。 昇順に近い入力ではさらに長くなり、最悪(厳密な降順)では M まで縮む。

本記事のデモとベンチマークでは、主記憶上の配列を入力ストリームとみなし、容量 M の最小ヒープでランを作ったあと、ラン同士をマージして配列へ書き戻す。

  1. 充填: 入力から最大 M 個を読み、最小ヒープを構築する。
  2. 抽出と置換: ヒープの最小を現在ランへ出力する。次の入力が直前の出力以上なら根へ入れて沈降(現在ランに残す)。小さければ次ラン用の待避領域へ置き、ヒープは縮む。
  3. ラン区切り: ヒープが空になったら現在ランを確定し、待避していた要素でヒープを組み直して次ランを始める。入力が尽きるまで繰り返す。
  4. 併合: できたランをマージソート同様に併合し、1 本の昇順列にする。
procedure sift_down(H, i)
  // 最小ヒープ条件を満たすよう H[i] を沈降

procedure generate_runs(input, M)
  H = first min(M, length(input)) elements; heapify_min(H)
  frozen = empty; run = empty; i = |H|
  while true
    if H is empty
      if run nonempty then emit run; run = empty
      if frozen empty and i >= length(input) then break
      H = frozen; frozen = empty
      while i < length(input) and |H| < M
        append input[i] to H; i = i + 1
      heapify_min(H)
      continue
    out = extract_min(H)
    append out to run
    if i < length(input)
      next = input[i]; i = i + 1
      if next >= out then insert next into H
      else append next to frozen
  return all emitted runs

procedure replacement_selection_sort(A)
  runs = generate_runs(A, M)
  A = merge_all(runs)

ラン生成は各要素がヒープへ高々定数回出入りするため O(n log M)、併合はラン数を R とすると概ね O(n log R) で、合計は O(n log n) 程度になる。 ヒープと待避・ラン用に O(M + n) の追加領域が要り、一般に不安定である。デモでは M = 4、ベンチマークでは M = 32 とする。

外部整列の初期ラン生成として置換選択を使い、できたランをポリフェーズマージなどで併合するのが古典的な組み合わせである。メモリ全体をヒープに使えるならランは 1 本になり、振る舞いは最小ヒープからの連続抽出に近づく。

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

ヒープソートは配列全体をヒープ化しインプレースで縮める。置換選択は容量 M の窓だけをヒープに保ち、ストリームからランを伸ばす。

トーナメントソート敗者木ソートは「次の最小」を木で更新する構造が近く、外部マージの選択木としても使われる。置換選択はラン長を伸ばす生成法としての側面が強い。

ストランドソートも単調列を切り取って併合するが、ヒープによる置換は行わず、1 回の走査で拾える非減少部分列に限る。

ポリフェーズマージソートは固定長チャンクを初期ランとする実装が多い。置換選択で長い初期ランを渡せば、併合パス数をさらに抑えられる。

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

Size Average time Maximum time Average memory Maximum memory
256 0.000022 0.001439 4 5
512 0.000045 0.001127 10 11
1024 0.000087 0.001155 20 23
2048 0.000185 0.003634 40 44
4096 0.000350 0.003149 81 86
8192 0.000734 0.004720 163 171
16384 0.002018 0.056631 326 336
32768 0.004071 0.037409 652 668
65536 0.007650 0.050894 1305 1332
131072 0.015410 0.096450 2609 2654
262144 0.032595 0.099702 5216 5291
計測に使用したコードを表示する

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 merge_values(left: &[usize], right: &[usize]) -> Vec<usize> {
    let mut out = Vec::with_capacity(left.len() + right.len());
    let (mut l, mut r) = (0, 0);
    while l < left.len() && r < right.len() {
        if left[l] <= right[r] {
            out.push(left[l]);
            l += 1;
        } else {
            out.push(right[r]);
            r += 1;
        }
    }
    out.extend_from_slice(&left[l..]);
    out.extend_from_slice(&right[r..]);
    out
}



const HEAP_SIZE: usize = 32;

fn sift_down(heap: &mut [usize], mut i: usize) {
    let n = heap.len();
    loop {
        let left = 2 * i + 1;
        let right = left + 1;
        let mut smallest = i;
        if left < n && heap[left] < heap[smallest] {
            smallest = left;
        }
        if right < n && heap[right] < heap[smallest] {
            smallest = right;
        }
        if smallest == i {
            break;
        }
        heap.swap(i, smallest);
        i = smallest;
    }
}

fn sift_up(heap: &mut [usize], mut i: usize) {
    while i > 0 {
        let parent = (i - 1) / 2;
        if heap[i] >= heap[parent] {
            break;
        }
        heap.swap(i, parent);
        i = parent;
    }
}

fn heapify_min(heap: &mut [usize]) {
    if heap.len() <= 1 {
        return;
    }
    for i in (0..heap.len() / 2).rev() {
        sift_down(heap, i);
    }
}

fn heap_push(heap: &mut Vec<usize>, value: usize) {
    heap.push(value);
    let i = heap.len() - 1;
    sift_up(heap, i);
}

fn heap_pop_min(heap: &mut Vec<usize>) -> usize {
    let n = heap.len();
    debug_assert!(n > 0);
    let min = heap[0];
    let last = heap.pop().unwrap();
    if !heap.is_empty() {
        heap[0] = last;
        sift_down(heap, 0);
    }
    min
}

fn generate_runs(input: &[usize], mem: usize) -> Vec<Vec<usize>> {
    let n = input.len();
    let mut runs = Vec::new();
    if n == 0 {
        return runs;
    }

    let m = mem.min(n).max(1);
    let mut i = 0usize;
    let mut heap = Vec::with_capacity(m);
    while i < n && heap.len() < m {
        heap.push(input[i]);
        i += 1;
    }
    heapify_min(&mut heap);

    let mut frozen = Vec::with_capacity(m);
    let mut run = Vec::new();

    loop {
        if heap.is_empty() {
            if !run.is_empty() {
                runs.push(std::mem::take(&mut run));
            }
            if frozen.is_empty() && i >= n {
                break;
            }
            heap = std::mem::take(&mut frozen);
            while i < n && heap.len() < m {
                heap.push(input[i]);
                i += 1;
            }
            if heap.is_empty() {
                break;
            }
            heapify_min(&mut heap);
            continue;
        }

        let out = heap_pop_min(&mut heap);
        run.push(out);

        if i < n {
            let next = input[i];
            i += 1;
            if next >= out {
                heap_push(&mut heap, next);
            } else {
                frozen.push(next);
            }
        }
    }

    runs
}

fn merge_all_runs(runs: Vec<Vec<usize>>) -> Vec<usize> {
    if runs.is_empty() {
        return Vec::new();
    }
    let mut queue = runs;
    while queue.len() > 1 {
        let mut next = Vec::with_capacity((queue.len() + 1) / 2);
        let mut idx = 0;
        while idx + 1 < queue.len() {
            next.push(merge_values(&queue[idx], &queue[idx + 1]));
            idx += 2;
        }
        if idx < queue.len() {
            next.push(queue[idx].clone());
        }
        queue = next;
    }
    queue.pop().unwrap_or_default()
}

fn replacement_selection_sort(a: &mut [usize]) {
    let n = a.len();
    if n <= 1 {
        return;
    }
    let runs = generate_runs(a, HEAP_SIZE);
    let sorted = merge_all_runs(runs);
    a.copy_from_slice(&sorted);
}


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

    replacement_selection_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),
        );
    }
}


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