フィボナッチヒープソートで配列を並び替える

フィボナッチヒープソートを使用する

フィボナッチヒープソート (fibonacci heap sort) は、要素をフィボナッチヒープへ挿入したあと、最小値を繰り返し取り出して昇順にする整列である。

フィボナッチヒープは、根が最小となる木の森である。各木はヒープ条件（親のキーが子以下）を満たし、根リスト上では同じ次数（子の本数）の木が複数あってもよい。二項ヒープと違い、挿入や合併では次数をすぐには揃えない。代わりに最小抽出のときに 統合（consolidate） で同次数の木を結合し、結果として根の次数は高々 1 本になる。木の形の上限がフィボナッチ数と結びつくためこの名がある。

挿入: 各要素を次数 0 の単独の根として根リストへ加える。統合は行わない。償却的に O(1)。
抽出: 根リストから最小キーの根を外し、その子たちを根リストへ移す。そのあと次数ごとの表で同次数の根を結合（キーの小さい方を親にする）して森を整える。償却的に O(log n)。
書き戻し: 取り出したキーを配列の先頭から順に書けば昇順になる。

procedure link(y, x)   // y.key >= x.key
  make y a child of x
  x.degree = x.degree + 1

procedure consolidate(H)
  // 次数 d ごとに高々 1 本の根が残るよう link で畳み込む
  for each root x in H
    while there is another root y with y.degree = x.degree
      if x.key > y.key
        swap x, y
      link(y, x)
  return the new root list

procedure fibonacci_heap_sort(A)
  H = empty fibonacci heap   // 根リスト
  for x in A
    add singleton(x) to root list of H
  for i from 0 to length(A) - 1
    m = minimum root in H
    remove m from root list
    add children of m to root list
    H = consolidate(H)
    A[i] = m.key

償却時間計算量は全体で O(n log n) であり、節点用に O(n) の追加記憶域が要る（インプレースではない）。等値キーの相対順序は結合時の規約に依存し、一般に不安定である。減少キーや削除を多用する優先度付きキューでは償却性能が強みになるが、整列だけなら統合のコストが毎回の抽出に乗り、ポインタ経由の森はキャッシュ効率では配列上のヒープソートに劣りやすい。

優先度付きキューとしてのフィボナッチヒープは、挿入・合併・減少キーが償却的にほぼ定数時間で書ける点が理論上の強みである。整列用途では減少キーを使わず、「すべて挿入してからすべて取り出す」形に固定したものがフィボナッチヒープソートである。

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

ヒープソートは配列上の二分ヒープをインプレースで縮める。フィボナッチヒープはポインタの森で、抽出時の統合が中心になる。

二項ヒープソートは挿入のたびに同次数を結合する。フィボナッチヒープは挿入を遅延し、抽出時にまとめて統合する。

ペアリングヒープソートは多分岐木の合併と子の二パス・ペアリングを使う。フィボナッチヒープは次数表による統合で森を整える。

左傾ヒープソートは単一の二分木とヌルパス長で左傾性を保つ。フィボナッチヒープは複数の根を持つ森である。

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

Size	Average time	Maximum time	Average memory	Maximum memory
256	0.000067	0.000307	10	10
512	0.000152	0.000372	20	20
1024	0.000337	0.000992	40	40
2048	0.000754	0.002167	80	80
4096	0.001687	0.003353	160	160
8192	0.003736	0.008631	320	320
16384	0.007950	0.017043	640	640
32768	0.017218	0.056944	1280	1280
65536	0.035678	0.123685	2560	2560
131072	0.085601	0.420435	5120	5120
262144	0.193887	0.415481	10240	10240

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

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 FibNode {
    key: usize,
    degree: u32,
    child: Option<Box<FibNode>>,
    sibling: Option<Box<FibNode>>,
}

fn link(mut child: Box<FibNode>, mut parent: Box<FibNode>) -> Box<FibNode> {
    child.sibling = parent.child.take();
    parent.child = Some(child);
    parent.degree += 1;
    parent
}

fn roots_to_vec(mut head: Option<Box<FibNode>>) -> Vec<Box<FibNode>> {
    let mut roots = Vec::new();
    while let Some(mut node) = head.take() {
        head = node.sibling.take();
        roots.push(node);
    }
    roots
}

fn vec_to_roots(roots: Vec<Box<FibNode>>) -> Option<Box<FibNode>> {
    let mut head: Option<Box<FibNode>> = None;
    let mut tail = &mut head;
    for node in roots {
        *tail = Some(node);
        tail = &mut tail.as_mut().unwrap().sibling;
    }
    head
}

fn consolidate(head: Option<Box<FibNode>>) -> Option<Box<FibNode>> {
    let roots = roots_to_vec(head);
    if roots.is_empty() {
        return None;
    }

    let mut degree_table: Vec<Option<Box<FibNode>>> = Vec::new();

    for mut x in roots {
        loop {
            let d = x.degree as usize;
            if d >= degree_table.len() {
                degree_table.resize_with(d + 1, || None);
            }
            if degree_table[d].is_none() {
                degree_table[d] = Some(x);
                break;
            }
            let y = degree_table[d].take().unwrap();
            x = if x.key <= y.key {
                link(y, x)
            } else {
                link(x, y)
            };
        }
    }

    let mut new_roots = Vec::new();
    for slot in degree_table {
        if let Some(node) = slot {
            new_roots.push(node);
        }
    }
    vec_to_roots(new_roots)
}

fn insert_key(heap: Option<Box<FibNode>>, key: usize) -> Option<Box<FibNode>> {
    let mut node = Box::new(FibNode {
        key,
        degree: 0,
        child: None,
        sibling: None,
    });
    node.sibling = heap;
    Some(node)
}

fn extract_min(heap: Option<Box<FibNode>>) -> (Option<usize>, Option<Box<FibNode>>) {
    let Some(head) = heap else {
        return (None, None);
    };

    let mut roots = roots_to_vec(Some(head));
    let mut min_i = 0;
    for i in 1..roots.len() {
        if roots[i].key < roots[min_i].key {
            min_i = i;
        }
    }

    let mut min_node = roots.remove(min_i);
    let key = min_node.key;
    let children = roots_to_vec(min_node.child.take());
    roots.extend(children);
    (Some(key), consolidate(vec_to_roots(roots)))
}

fn fibonacci_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("fibonacci heap exhausted early");
    }
}


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

    fibonacci_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