ペアリングヒープソートを使用する

ペアリングヒープソート (pairing heap sort) は、要素をペアリングヒープへ挿入したあと、最小値を繰り返し取り出して昇順にする整列である。

ペアリングヒープは、根が最小(または最大)となる多分岐木である。子は左子・右兄弟形式で並べ、ヒープ条件は「親のキーはどの子のキー以下」だけを課す。二項ヒープのように次数を揃えたり、フィボナッチヒープのようにランクを管理したりしない代わりに、合併(meld)子のペアリング で構造を保つ。

  1. 合併: 2 本のヒープの根を比較し、キーの大きい方を小さい方の最左の子にする。比較は 1 回で、償却的に O(1)
  2. 挿入: 単一節点のヒープを既存ヒープと合併する。
  3. 抽出: 根を外し、その子たちを左から 2 本ずつ合併(第 1 パス)したあと、できたヒープを右から順に合併(第 2 パス)して新しい根を得る。償却的に O(log n)
  4. 書き戻し: 取り出したキーを配列の先頭から順に書けば昇順になる。
procedure meld(H1, H2)
  if H1 is empty
    return H2
  if H2 is empty
    return H1
  if H1.key <= H2.key
    make H2 the leftmost child of H1
    return H1
  else
    make H1 the leftmost child of H2
    return H2

procedure two_pass_meld(children)   // children は兄弟リスト
  // 第1パス: 左から隣り合う 2 本を meld
  pairs = empty list
  while children is not empty
    a = take first child
    if children is empty
      append a to pairs
    else
      b = take first child
      append meld(a, b) to pairs
  // 第2パス: 右から順に meld
  H = empty
  for p in reverse(pairs)
    H = meld(p, H)
  return H

procedure pairing_heap_sort(A)
  H = empty pairing heap
  for x in A
    H = meld(H, singleton(x))
  for i from 0 to length(A) - 1
    A[i] = H.key
    H = two_pass_meld(children of H)

償却時間計算量は挿入・合併が O(1)、最小抽出が O(log n) であり、全体では O(n log n) になる。節点用に O(n) の追加記憶域が要る(インプレースではない)。等値キーの相対順序は合併時の規約に依存し、一般に不安定である。実装が単純な一方で、ポインタ経由の多分岐木はキャッシュ効率では配列上のヒープソートに劣りやすい。

優先度付きキューとしてのペアリングヒープは、合併が 1 比較で書け、減少キーも「切り離して再合併」で扱える点が実務寄りの実装で選ばれやすい。整列用途ではその操作を「すべて挿入してからすべて取り出す」形に固定したものがペアリングヒープソートである。

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

ヒープソートは配列上の二分ヒープをインプレースで縮める。ペアリングヒープはポインタの多分岐木で、合併と子のペアリングが中心になる。

二項ヒープソートは次数の異なる二項木を二進加算のように結合する。ペアリングヒープは次数を持たず、抽出時に子を 2 本ずつペアにして畳み込む。

弱ヒープソートは配列上の不完全木と逆ビットで比較回数を抑える。ヒープ同士の合併を第一級には扱わない。

トーナメントソートも最小を繰り返し取り出すが、固定長のトーナメント木を更新する点が異なる。

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

Size Average time Maximum time Average memory Maximum memory
256 0.000018 0.000399 6 7
512 0.000039 0.000132 13 14
1024 0.000087 0.000370 27 28
2048 0.000177 0.000392 54 56
4096 0.000416 0.001518 110 112
8192 0.000897 0.002152 220 224
16384 0.002085 0.007597 442 448
32768 0.004604 0.012753 884 896
65536 0.012509 0.054332 1770 1792
131072 0.053413 0.492478 3540 3584
262144 0.065955 0.239336 7083 7168
計測に使用したコードを表示する

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


struct PairingNode {
    key: usize,
    child: Option<Box<PairingNode>>,
    sibling: Option<Box<PairingNode>>,
}

fn meld(
    a: Option<Box<PairingNode>>,
    b: Option<Box<PairingNode>>,
) -> Option<Box<PairingNode>> {
    match (a, b) {
        (None, x) | (x, None) => x,
        (Some(mut x), Some(mut y)) => {
            if x.key <= y.key {
                y.sibling = x.child.take();
                x.child = Some(y);
                Some(x)
            } else {
                x.sibling = y.child.take();
                y.child = Some(x);
                Some(y)
            }
        }
    }
}

fn two_pass_meld(mut first: Option<Box<PairingNode>>) -> Option<Box<PairingNode>> {
    let mut pairs = Vec::new();
    while let Some(mut a) = first.take() {
        first = a.sibling.take();
        if let Some(mut b) = first.take() {
            first = b.sibling.take();
            pairs.push(meld(Some(a), Some(b)));
        } else {
            pairs.push(Some(a));
        }
    }

    let mut result = None;
    for pair in pairs.into_iter().rev() {
        result = meld(pair, result);
    }
    result
}

fn insert_key(heap: Option<Box<PairingNode>>, key: usize) -> Option<Box<PairingNode>> {
    let node = Box::new(PairingNode {
        key,
        child: None,
        sibling: None,
    });
    meld(heap, Some(node))
}

fn extract_min(heap: Option<Box<PairingNode>>) -> (Option<usize>, Option<Box<PairingNode>>) {
    let Some(mut root) = heap else {
        return (None, None);
    };
    let key = root.key;
    let children = root.child.take();
    (Some(key), two_pass_meld(children))
}

fn pairing_heap_sort(a: &mut [usize]) {
    let mut heap = None;
    for &key in a.iter() {
        heap = insert_key(heap, key);
    }
    for slot in a.iter_mut() {
        let (key, next) = extract_min(heap);
        heap = next;
        *slot = key.expect("pairing heap exhausted early");
    }
}


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

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


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 KiB (explicit buffers such as swap).
fn run_once(size: usize, seed: usize) -> (u128, usize) {
    let mut array = input_array(size, seed as u64);

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

    let start = Instant::now();

    benchmark_sort(&mut array);

    let elapsed = start.elapsed();
    let peak_bytes = PEAK_BYTES.load(Ordering::Relaxed);
    let aux_kb = peak_bytes.saturating_sub(base_bytes) / 1024;

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

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

    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;
        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