マージ挿入ソートを使用する

マージ挿入ソート (merge-insertion sort) は、要素をペアに分け、各ペアの大きい方を再帰的に整列し、小さい方を二分探索で挿入する。

Ford-Johnson アルゴリズムとしても知られる。

要素をペアに分けて比較し、各ペアの大きい方を主系列 (main chain)、小さいを保留列 (pend) として分ける。主系列を再帰的に整列し、保留列の要素をヤーコプスタール数に基づく順序で二分探索挿入することで、比較回数の上界が抑えられる。

  1. ペアリング: 隣接要素を ⌊n/2⌋ 組のペアに分け、 (小, 大) の順に並べる。要素数が奇数なら最後の1要素は保留列に回す。
  2. 再帰: 各ペアの大きい方だけを取り出し、同じ手順で再帰的に整列する。
  3. 主系列の形成: 整列済みの大きい方の列の先頭に、先頭ペアの小さい方を置く。
  4. 保留列の挿入: 残りの小さい方を、ヤーコプスタール数で決めた順序で1つずつ主系列へ二分探索挿入する。保留列の各要素は、対応するペアの大きい方より左の区間だけを探索範囲とする。
procedure merge_insertion_sort(A)
  n = length(A)
  if n <= 1 then return
  if n = 2 then
    compare and swap A[0], A[1] if needed
    return
  form pairs (small, large) by comparing adjacent elements
  L = array of large values from each pair
  merge_insertion_sort(L)
  reorder pairs by sorted L
  chain = [small of first pair] followed by sorted larges
  pend = remaining smalls (and odd element if any)
  for each index i in jacobsthal_insertion_order(length(pend))
    limit = position of paired large for pend[i] in chain
    pos = binary_search(chain[0 .. limit), pend[i])
    insert pend[i] at pos in chain
  copy chain back into A

ヤーコプスタール数 JJ₀ = 0, J₁ = 1, Jₙ = Jₙ₋₁ + 2Jₙ₋₂ で定義され、J₂ 以降は 1, 3, 5, 11, 21, 43, … と続く。

保留列の1番目から数えて J₂, J₃, … の位置を先に挿入し、各ヤーコプスタール番号の間にある残りを降順で埋める。たとえば保留列が6要素なら挿入順は 1, 3, 2, 5, 4, 6 となる。

比較回数が理論的下界に近い O(n log n) となるが、一般に不安定である。

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

挿入ソートO(n²) 比較を、ペアリングと主系列への二分挿入で抑える。実装は複雑だが比較回数の最小化が目的である。

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

Size Average time Maximum time Average memory Maximum memory
256 0.000061 0.000748 78 84
512 0.000164 0.000915 98 104
1024 0.000473 0.001836 142 148
2048 0.001597 0.004256 170 176
4096 0.005933 0.015892 266 272
8192 0.022476 0.051114 406 412
16384 0.087843 0.133609 683 788
32768 0.386081 0.599429 1512 1556
計測に使用したコードを表示する

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 = 15;
const RUNS: usize = 8192;


fn ford_johnson_insert_order(m: usize) -> Vec<usize> {
    if m == 0 {
        return Vec::new();
    }
    let mut js = vec![0usize, 1];
    while *js.last().unwrap() < m {
        let l = js.len();
        js.push(js[l - 1] + 2 * js[l - 2]);
    }
    let mut order = Vec::new();
    let mut used = vec![false; m];
    let mut prev_j = 0usize;
    for &j in js.iter().skip(1) {
        if j > m {
            break;
        }
        if j > prev_j {
            for idx in (prev_j..=j - 1).rev() {
                if idx < m && !used[idx] {
                    order.push(idx);
                    used[idx] = true;
                }
            }
            prev_j = j;
        }
    }
    for idx in (0..m).rev() {
        if !used[idx] {
            order.push(idx);
        }
    }
    order
}

fn ford_johnson_reorder_pairs(
    pairs: &[(usize, usize)],
    sorted_larges: &[usize],
) -> Vec<(usize, usize)> {
    let mut out = Vec::with_capacity(sorted_larges.len());
    let mut taken = vec![false; pairs.len()];
    for &lg in sorted_larges {
        for (i, p) in pairs.iter().enumerate() {
            if !taken[i] && p.1 == lg {
                out.push(*p);
                taken[i] = true;
                break;
            }
        }
    }
    out
}

fn ford_johnson(a: &mut [usize]) {
    let n = a.len();
    if n <= 1 {
        return;
    }
    if n == 2 {
        if a[0] > a[1] {
            a.swap(0, 1);
        }
        return;
    }
    let pair_count = n / 2;
    let mut pairs: Vec<(usize, usize)> = Vec::with_capacity(pair_count);
    for i in 0..pair_count {
        let lo = 2 * i;
        let hi = lo + 1;
        if a[lo] > a[hi] {
            pairs.push((a[hi], a[lo]));
        } else {
            pairs.push((a[lo], a[hi]));
        }
    }
    let odd = if n % 2 == 1 { Some(a[n - 1]) } else { None };
    let mut larges: Vec<usize> = pairs.iter().map(|p| p.1).collect();
    ford_johnson(&mut larges);
    let sorted_pairs = ford_johnson_reorder_pairs(&pairs, &larges);
    let mut chain = Vec::with_capacity(n);
    chain.push(sorted_pairs[0].0);
    chain.extend(sorted_pairs.iter().map(|p| p.1));
    let mut pending_pairs: Vec<(usize, usize)> =
        sorted_pairs.iter().skip(1).copied().collect();
    if let Some(v) = odd {
        pending_pairs.push((v, usize::MAX));
    }
    let pending: Vec<usize> = pending_pairs.iter().map(|p| p.0).collect();
    for idx in ford_johnson_insert_order(pending.len()) {
        let val = pending[idx];
        let limit = if pending_pairs[idx].1 == usize::MAX {
            chain.len()
        } else {
            chain
                .iter()
                .position(|&x| x == pending_pairs[idx].1)
                .unwrap_or(chain.len())
        };
        let pos = chain[..limit].partition_point(|&x| x < val);
        chain.insert(pos, val);
    }
    a.copy_from_slice(&chain);
}


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

    ford_johnson(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