弱ヒープソートで配列を並び替える
弱ヒープソートを使用する
弱ヒープソート (weak-heap sort) は、配列を弱ヒープ(weak heap)に整えたあと、根(最大値)と末尾を入れ替えてヒープを縮めていく整列である。
手順はヒープソートと同じだが、二分ヒープより緩い順序条件と、左右の子の役割を示す逆ビット(reverse bit)を使う点が異なる。
弱ヒープを二分木として見ると次を満たし、この記事では根が最大となる最大ヒープ形となるように構築する。
- 根は左の子を持たない(右の子だけを持つ)。
- 各節点の値は、その右部分木に属するすべての値以上である(左の子側=兄弟列には直接の大小を課さない)。
- 葉は最下層かそのひとつ上にだけ現れる(完全二分木と同じ配置)。
配列表現では、節点 i の逆ビット r[i] が「左の子」と「右の子」のどちらを 2i / 2i+1 に割り当てるかを決める。r[i] = 0 なら左の子は 2i、右の子は 2i+1、r[i] = 1 なら入れ替わる。根(i = 0)は常に右の子 1 だけを見る。
多分岐ヒープとして見ると、右の子は「最初の子」、左の子は「次の兄弟」に対応し、二項ヒープの木を 1 本の不完全木にまとめた形になる。ある節点 j の多分岐上の親を区別祖先(distinguished ancestor)と呼び、ヒープ条件は「区別祖先の値が j の値以上」に落ちる。
- 構築: 逆ビットをすべて 0 にし、末尾から 1 まで各節点
jをその区別祖先と 結合 する。結合は 1 回の比較で、子の方が大きければ交換し、子側の逆ビットを反転する。全体でちょうどn - 1回の比較で弱ヒープになる。 - 抽出: 根とヒープ末尾を交換して最大値を確定する。
- 沈降: 新しい根について、右部分木の左背骨を葉まで下り、そこから親へ遡りながら根と結合を繰り返す。二分ヒープの沈降が各段で最大 2 比較なのに対し、弱ヒープでは高さぶんの比較で足りる。
- 反復: ヒープ長が 2 になるまで手順 2〜3 を繰り返し、最後に残った 2 要素を入れ替えて昇順を完成する。
procedure distinguished_ancestor(r, j)
while (j & 1) = r[j >> 1]
j = j >> 1
return j >> 1
procedure join(A, r, i, j) // i は区別祖先、最大ヒープ
if A[i] < A[j]
flip r[j]
swap A[i], A[j]
procedure weak_heap_sort(A)
n = length(A)
r[0..n) = 0
for j from n - 1 down to 1
join(A, r, distinguished_ancestor(r, j), j)
for end from n - 1 down to 2
swap A[0], A[end]
x = 1 // 根の右の子
while 2 * x + r[x] < end // 右部分木の左背骨を下る
x = 2 * x + r[x]
while x > 0
join(A, r, 0, x)
x = x >> 1
swap A[0], A[1]
最悪時間計算量は O(n log n) である。構築は n - 1 比較、抽出フェーズの比較回数はおよそ n ⌈log₂ n⌉ 前後に抑えられ、通常のヒープソート(沈降で最大約 2 n log n 比較)より比較が少なくなりやすい。
逆ビットに O(n) ビットの追加領域が要る(厳密なインプレースではない)。等値の扱いは結合時の規約依存で、一般に不安定である。
類似アルゴリズムとの相違点
ヒープソートは親子両方にヒープ条件を課し、沈降で最大 2 比較/段を使う。弱ヒープは右部分木だけに大小を課し、逆ビットで左右を入れ替えられるため、沈降の比較回数をおよそ半分に近づけられる。
二項ヒープソートは次数の異なる二項木の森として合併する。完全な弱ヒープ(要素数 2^k)は単一の二項木と同型だが、弱ヒープは不完全な 1 本の木のまま扱う。
トーナメントソートや敗者木ソートは比較結果を木に蓄えて最小を繰り返し取り出す方式で、配列上の弱ヒープ構築+末尾確定とは手順が異なる。
計算時間量および空間計算量を計測する
| Size | Average time | Maximum time | Average memory | Maximum memory |
|---|---|---|---|---|
| 256 | 0.000011 | 0.000127 | 0 | 0 |
| 512 | 0.000025 | 0.007365 | 0 | 0 |
| 1024 | 0.000051 | 0.000589 | 0 | 0 |
| 2048 | 0.000109 | 0.000357 | 0 | 0 |
| 4096 | 0.000233 | 0.010780 | 0 | 0 |
| 8192 | 0.000514 | 0.001743 | 1 | 1 |
| 16384 | 0.001121 | 0.008392 | 2 | 2 |
| 32768 | 0.002406 | 0.008379 | 4 | 4 |
| 65536 | 0.005158 | 0.007401 | 8 | 8 |
| 131072 | 0.011057 | 0.016855 | 16 | 16 |
| 262144 | 0.023785 | 0.097358 | 32 | 32 |
計測に使用したコードを表示する
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;
fn get_flag(r: &[u8], x: usize) -> usize {
((r[x >> 3] >> (x & 7)) & 1) as usize
}
fn toggle_flag(r: &mut [u8], x: usize) {
r[x >> 3] ^= 1 << (x & 7);
}
/// Join two equal-height weak heaps rooted at `i` (distinguished ancestor) and `j`.
/// Max-heap form: if `a[j]` is larger, promote it and flip the reverse bit at `j`.
fn join(a: &mut [usize], r: &mut [u8], i: usize, j: usize) {
if a[i] < a[j] {
toggle_flag(r, j);
a.swap(i, j);
}
}
fn distinguished_ancestor(r: &[u8], mut j: usize) -> usize {
while (j & 1) == get_flag(r, j >> 1) {
j >>= 1;
}
j >> 1
}
fn weak_heap_sort(a: &mut [usize]) {
let n = a.len();
if n <= 1 {
return;
}
let mut r = vec![0u8; (n + 7) / 8];
// Bottom-up construct: n - 1 joins with each node's distinguished ancestor.
for i in (1..n).rev() {
let g = distinguished_ancestor(&r, i);
join(a, &mut r, g, i);
}
// Extract maxima like heapsort; sift-down uses left-spine + upward joins.
for end in (2..n).rev() {
a.swap(0, end);
let mut x = 1usize;
while {
let y = 2 * x + get_flag(&r, x);
y < end
} {
x = 2 * x + get_flag(&r, x);
}
while x > 0 {
join(a, &mut r, 0, x);
x >>= 1;
}
}
a.swap(0, 1);
}
fn benchmark_sort(array: &mut [usize]) {
weak_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