スムースソートで配列を並び替える
スムースソートを使う
スムースソート (smooth sort) は、レオナルド木を積み重ねた森を維持し、末尾を最大値として取り出していく。
ヒープソート (heap sort) と同様に最悪時間計算量は O(n log n) で、追加の補助記憶域は O(1) に抑えられるインプレースなアルゴリズムだが、入力がすでに昇順に近いほど作業量が少なくなる適応的 (adaptive) 性質を持つ。
通常の二分ヒープでは、根が配列の一端にあり、最大値を取り出すたびに根と末尾を交換してヒープを底側から縮めていく。この配置は「ヒープとしての整理」と「昇順に書き出す向き」が噛み合わず、すでにソート済みの入力でも一度ヒープ状に混ぜ直すため、最良でも必ず O(n log n) 程度の仕事が残る。
スムースソートでは、二分木ではなくレオナルド数に基づくレオナルド木を積み重ねた森を維持する。各木はヒープ条件を満たし、さらに各木の根の値は左から右へ弱い順（非減少）として並べる。そうすると常に右端が全体の最大になり、末尾から確定した最大値を「取り出す」処理がすでに整っている入力ではほとんど余計な比較を要しない。
レオナルド数 L(k) は
L(0) = L(1) = 1L(k) = L(k-2) + L(k-1) + 1（k ≥ 2）
で定められ、1, 1, 3, 5, 9, 15, … と続く。サイズ L(k) のレオナルド木は、大きい方の子が左、小さい方が右になるよう二つのより小さいレオナルド木と根で構成される。
任意の長さ n は、高々 O(log n) 個の互いに異なるレオナルド数の和として表せる。スムースソートはこの性質を利用して森に含まれる木の本数を常に対数オーダーに抑え、整列済みに近い入力では O(n) に近づき、最悪計算量でも O(n log n) を保つ適応型ソートとなる。
- 森の表現: 左から右へ根が非減少になるよう、互いに異なるサイズのレオナルド木を複数本並べた森を維持する。どの位置に木の根があるかはビットマスクと末尾木のサイズ（オフセット）で表す。
- 第1段階（構築）: インデックス
iを 1 からn - 1まで増やし、森の符号を更新したうえで新根位置iに対し、単一木内の沈下sift_inまたは隣木との整合interheap_siftでヒープ条件を満たす。 - 第2段階（確定）: 右端は常に全体の最大なので、
iをn - 1から 2 まで減らしながら末尾木を縮小または 2 子木へ分割し、生じた根にinterheap_siftを適用して森の不変条件を保つ。 - 終了: 森に残る 2 要素以下は根が左から右へ非減少のため、すでに昇順であり追加の sift は不要。
procedure sift_in(A, rootIdx, size)
// インデックス size のレオナルド木 1 本の中でヒープ条件を満たすよう根から下げる。
// L[k]=L[k-2]+L[k-1]+1、子方向は二分木のサイズにより左右どちらかへ潜る。
if size < 2 then return
tmp = A[rootIdx]
r = rootIdx
sz = size
loop
right = r - 1
left = right - L[sz - 2]
if A[right] < A[left]
candidate = left
nsz = sz - 1
else
candidate = right
nsz = sz - 2
if A[candidate] <= tmp then break
A[r] = A[candidate]
r = candidate
sz = nsz
if sz <= 1 then break
A[r] = tmp
procedure interheap_sift(A, rootIdx, heap_state)
// heap_state は「どの桁に木の根があるか」（マスク）と「現在見ている木の順 order」（オフセット）などを束ねたもの。
tmp = A[rootIdx]
r = rootIdx
state = heap_state anchored at rootIdx
loop while mask(state) <> 1
maxValue = tmp
if order(state) > 1 then
right = r - 1
left = right - L[order(state) - 2]
maxValue = max(maxValue, A[left], A[right])
next = r - L[order(state)]
if A[next] <= maxValue then break
A[r] = A[next]
r = next
// マスクを調整して左隣の木へカーソル移動
march_one_tree_left(state)
A[r] = tmp
sift_in(A, r, order(state))
procedure smooth_sort(A)
n = length(A)
if n <= 1 then return
heap_state ← initial_singleton_forest_encoding()
// 第 1 段階（ヒープ化）
for i = 1 to n - 1
advance_leonardo_forest(heap_state, i)
if wide_bottom_condition(i, heap_state, n) then sift_in(A, i, order_at_tail(heap_state)) else interheap_sift(A, i, heap_state)
// 第 2 段階（右端から順に確定・森の更新）
for i = n - 1 down to 2
if order_at_tail(heap_state) < 2 then shrink_trailing_trees(heap_state)
else
split_rightmost_into_two_children(heap_state, i)
for child_root c in newborn_pair_roots()
interheap_sift(A, c, heap_state)
実装の複雑さと定数倍の大きさから、汎用ライブラリの sort として採用されることは稀である。
類似アルゴリズムとの相違点
ヒープソートと同様にインプレースで最悪計算量 O(n log n) だが、レオナルド木の森により整列済み入力では O(n) に近づく適応型ソートである。
計算時間量および空間計算量を計測する
| Size | Average time | Maximum time | Average memory | Maximum memory |
|---|---|---|---|---|
| 256 | 0.000012 | 0.000442 | 1662 | 1668 |
| 512 | 0.000028 | 0.000354 | 1666 | 1672 |
| 1024 | 0.000061 | 0.000355 | 1674 | 1680 |
| 2048 | 0.000130 | 0.000469 | 1689 | 1696 |
| 4096 | 0.000285 | 0.000971 | 1721 | 1728 |
| 8192 | 0.000597 | 0.000918 | 1786 | 1792 |
| 16384 | 0.001261 | 0.002096 | 1918 | 1924 |
| 32768 | 0.002709 | 0.006676 | 2178 | 2184 |
| 65536 | 0.005773 | 0.010908 | 2690 | 2696 |
| 131072 | 0.012651 | 0.017563 | 3714 | 3720 |
| 262144 | 0.027944 | 0.070136 | 5762 | 5768 |
計測に使用したコードを表示する
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 = 18;
const RUNS: usize = 8192;
const LEONARDO: [usize; 46] = [
1, 1, 3, 5, 9, 15, 25, 41, 67, 109, 177, 287, 465, 753, 1219, 1973, 3193,
5167, 8361, 13529, 21891, 35421, 57313, 92735, 150049, 242785, 392835,
635621, 1028457, 1664079, 2692537, 4356617, 7049155, 11405773, 18454929,
29860703, 48315633, 78176337, 126491971, 204668309, 331160281, 535828591,
866988873, 1402817465, 2269806339, 3672623805,
];
fn smooth_sift_in(a: &mut [usize], root_idx: usize, size: usize) {
if size < 2 {
return;
}
let tmp = a[root_idx];
let mut root = root_idx;
let mut sz = size;
loop {
let right = root - 1;
let left = right - LEONARDO[sz - 2];
let (next, next_size) = if a[right] < a[left] {
(left, sz - 1)
} else {
(right, sz - 2)
};
if a[next] <= tmp {
break;
}
a[root] = a[next];
root = next;
sz = next_size;
if sz <= 1 {
break;
}
}
a[root] = tmp;
}
fn smooth_interheap_sift(a: &mut [usize], root_idx: usize, mask: usize, offset: usize) {
let tmp = a[root_idx];
let mut root = root_idx;
let mut hmask = mask;
let mut hoffset = offset;
while hmask != 1 {
let mut max = tmp;
if hoffset > 1 {
let right = root - 1;
let left = right - LEONARDO[hoffset - 2];
max = max.max(a[left]).max(a[right]);
}
let next = root - LEONARDO[hoffset];
if a[next] <= max {
break;
}
a[root] = a[next];
root = next;
loop {
hmask >>= 1;
hoffset += 1;
if hmask & 1 != 0 {
break;
}
}
}
a[root] = tmp;
smooth_sift_in(a, root, hoffset);
}
fn smooth_sort(a: &mut [usize]) {
let n = a.len();
if n <= 1 {
return;
}
let mut mask = 1usize;
let mut offset = 1usize;
for i in 1..n {
if mask & 2 != 0 {
mask = (mask >> 2) | 1;
offset += 2;
} else if offset == 1 {
mask = (mask << 1) | 1;
offset = 0;
} else {
mask = (mask << (offset - 1)) | 1;
offset = 1;
}
let wide_bottom =
(mask & 2 != 0 && i + 1 < n)
|| (offset > 0 && 1 + i + LEONARDO[offset - 1] < n);
if wide_bottom {
smooth_sift_in(a, i, offset);
} else {
smooth_interheap_sift(a, i, mask, offset);
}
}
for i in (2..n).rev() {
if offset < 2 {
loop {
mask >>= 1;
offset += 1;
if mask & 1 != 0 {
break;
}
}
} else {
let ch1 = i - 1;
let ch0 = ch1 - LEONARDO[offset - 2];
mask &= !1;
for ch in [ch0, ch1] {
mask = (mask << 1) | 1;
offset -= 1;
smooth_interheap_sift(a, ch, mask, offset);
}
}
}
}
fn benchmark_sort(array: &mut [usize]) {
smooth_sort(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 {
// 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;
}
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