マージ挿入ソートで配列を並び替える

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

マージ挿入ソート (merge-insertion sort) は、Lester Randolph Ford Jr. と Selmer Martin Johnson が 1959 年に提案した比較ソートで、Ford-Johnson アルゴリズム としても知られる。

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

ペアリング: 隣接要素を ⌊n/2⌋ 組のペアに分け、 (小, 大) の順に並べる。要素数が奇数なら最後の1要素は保留列に回す。
再帰: 各ペアの大きい方だけを取り出し、同じ手順で再帰的に整列する。
主系列の形成: 整列済みの大きい方の列の先頭に、先頭ペアの小さい方を置く。
保留列の挿入: 残りの小さい方を、ヤーコプスタール数で決めた順序で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

ヤーコプスタール数 J は J₀ = 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) で、小さな n において比較ソートの理論的下界に近い上界を達成する。実装では主系列や保留列用に O(n) の追加配列を使うことが多く、安定ソートではない。

挿入ソート単体の O(n²) と比べ、ペアリングと主系列による再帰で比較を節約する。マージソートやクイックソートのように実装が単純ではない反面、比較回数の最小化 を目的としたアルゴリズムとして使用される。

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

Size	Average time	Maximum time	Average memory	Maximum memory
256	0.000048	0.000962	1682	1688
512	0.000128	0.001177	1702	1708
1024	0.000406	0.001633	1738	1744
2048	0.001584	0.022723	1802	1808
4096	0.005493	0.037684	1930	1936
8192	0.019735	0.054744	2122	2128
16384	0.066834	0.556116	2559	2560
32768	0.282587	1.051091	3587	3592

計測に使用したコードを表示する

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 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 {
    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 run_once(size: usize, seed: usize) -> (u128, usize) {
    let expected: Vec<usize> = (1..=size).collect();
    let mut array = shuffled(size, seed as u64);

    let start = Instant::now();

    benchmark_sort(&mut array);

    let elapsed = start.elapsed();

    if array != expected {
        panic!(
            "sort failed with seed {} for size {}",
            seed,
            size
        );
    }

    (micros(elapsed), memory_usage_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;
    }

    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 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 += mem;

            if mem > max_mem {
                max_mem = 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