置換選択ソートで配列を並び替える
置換選択ソートを使用する
置換選択ソート (replacement selection sort) は、限られた大きさの最小ヒープで入力を流し、平均でヒープ容量の約 2 倍の長さの整列済みランを生成し、それらを併合して全体を昇順にする整列である。
外部整列ではメモリに載らないファイルを扱うとき、単純にメモリ分だけ読んでクイックソートすると初期ラン長はメモリ容量 M に留まる。
置換選択では、出力した直前のキー以上の入力だけをヒープへ戻す(置換する)ことで、ランダム入力でも期待ラン長がおよそ 2M になる。
昇順に近い入力ではさらに長くなり、最悪(厳密な降順)では M まで縮む。
本記事のデモとベンチマークでは、主記憶上の配列を入力ストリームとみなし、容量 M の最小ヒープでランを作ったあと、ラン同士をマージして配列へ書き戻す。
- 充填: 入力から最大
M個を読み、最小ヒープを構築する。 - 抽出と置換: ヒープの最小を現在ランへ出力する。次の入力が直前の出力以上なら根へ入れて沈降(現在ランに残す)。小さければ次ラン用の待避領域へ置き、ヒープは縮む。
- ラン区切り: ヒープが空になったら現在ランを確定し、待避していた要素でヒープを組み直して次ランを始める。入力が尽きるまで繰り返す。
- 併合: できたランをマージソート同様に併合し、1 本の昇順列にする。
procedure sift_down(H, i)
// 最小ヒープ条件を満たすよう H[i] を沈降
procedure generate_runs(input, M)
H = first min(M, length(input)) elements; heapify_min(H)
frozen = empty; run = empty; i = |H|
while true
if H is empty
if run nonempty then emit run; run = empty
if frozen empty and i >= length(input) then break
H = frozen; frozen = empty
while i < length(input) and |H| < M
append input[i] to H; i = i + 1
heapify_min(H)
continue
out = extract_min(H)
append out to run
if i < length(input)
next = input[i]; i = i + 1
if next >= out then insert next into H
else append next to frozen
return all emitted runs
procedure replacement_selection_sort(A)
runs = generate_runs(A, M)
A = merge_all(runs)
ラン生成は各要素がヒープへ高々定数回出入りするため O(n log M)、併合はラン数を R とすると概ね O(n log R) で、合計は O(n log n) 程度になる。
ヒープと待避・ラン用に O(M + n) の追加領域が要り、一般に不安定である。デモでは M = 4、ベンチマークでは M = 32 とする。
外部整列の初期ラン生成として置換選択を使い、できたランをポリフェーズマージなどで併合するのが古典的な組み合わせである。メモリ全体をヒープに使えるならランは 1 本になり、振る舞いは最小ヒープからの連続抽出に近づく。
類似アルゴリズムとの相違点
ヒープソートは配列全体をヒープ化しインプレースで縮める。置換選択は容量 M の窓だけをヒープに保ち、ストリームからランを伸ばす。
トーナメントソートや敗者木ソートは「次の最小」を木で更新する構造が近く、外部マージの選択木としても使われる。置換選択はラン長を伸ばす生成法としての側面が強い。
ストランドソートも単調列を切り取って併合するが、ヒープによる置換は行わず、1 回の走査で拾える非減少部分列に限る。
ポリフェーズマージソートは固定長チャンクを初期ランとする実装が多い。置換選択で長い初期ランを渡せば、併合パス数をさらに抑えられる。
計算時間量および空間計算量を計測する
| Size | Average time | Maximum time | Average memory | Maximum memory |
|---|---|---|---|---|
| 256 | 0.000022 | 0.001439 | 4 | 5 |
| 512 | 0.000045 | 0.001127 | 10 | 11 |
| 1024 | 0.000087 | 0.001155 | 20 | 23 |
| 2048 | 0.000185 | 0.003634 | 40 | 44 |
| 4096 | 0.000350 | 0.003149 | 81 | 86 |
| 8192 | 0.000734 | 0.004720 | 163 | 171 |
| 16384 | 0.002018 | 0.056631 | 326 | 336 |
| 32768 | 0.004071 | 0.037409 | 652 | 668 |
| 65536 | 0.007650 | 0.050894 | 1305 | 1332 |
| 131072 | 0.015410 | 0.096450 | 2609 | 2654 |
| 262144 | 0.032595 | 0.099702 | 5216 | 5291 |
計測に使用したコードを表示する
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 as AtomicOrdering},
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, AtomicOrdering::Relaxed) + size;
PEAK_BYTES.fetch_max(live, AtomicOrdering::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(), AtomicOrdering::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(), AtomicOrdering::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 merge_values(left: &[usize], right: &[usize]) -> Vec<usize> {
let mut out = Vec::with_capacity(left.len() + right.len());
let (mut l, mut r) = (0, 0);
while l < left.len() && r < right.len() {
if left[l] <= right[r] {
out.push(left[l]);
l += 1;
} else {
out.push(right[r]);
r += 1;
}
}
out.extend_from_slice(&left[l..]);
out.extend_from_slice(&right[r..]);
out
}
const HEAP_SIZE: usize = 32;
fn sift_down(heap: &mut [usize], mut i: usize) {
let n = heap.len();
loop {
let left = 2 * i + 1;
let right = left + 1;
let mut smallest = i;
if left < n && heap[left] < heap[smallest] {
smallest = left;
}
if right < n && heap[right] < heap[smallest] {
smallest = right;
}
if smallest == i {
break;
}
heap.swap(i, smallest);
i = smallest;
}
}
fn sift_up(heap: &mut [usize], mut i: usize) {
while i > 0 {
let parent = (i - 1) / 2;
if heap[i] >= heap[parent] {
break;
}
heap.swap(i, parent);
i = parent;
}
}
fn heapify_min(heap: &mut [usize]) {
if heap.len() <= 1 {
return;
}
for i in (0..heap.len() / 2).rev() {
sift_down(heap, i);
}
}
fn heap_push(heap: &mut Vec<usize>, value: usize) {
heap.push(value);
let i = heap.len() - 1;
sift_up(heap, i);
}
fn heap_pop_min(heap: &mut Vec<usize>) -> usize {
let n = heap.len();
debug_assert!(n > 0);
let min = heap[0];
let last = heap.pop().unwrap();
if !heap.is_empty() {
heap[0] = last;
sift_down(heap, 0);
}
min
}
fn generate_runs(input: &[usize], mem: usize) -> Vec<Vec<usize>> {
let n = input.len();
let mut runs = Vec::new();
if n == 0 {
return runs;
}
let m = mem.min(n).max(1);
let mut i = 0usize;
let mut heap = Vec::with_capacity(m);
while i < n && heap.len() < m {
heap.push(input[i]);
i += 1;
}
heapify_min(&mut heap);
let mut frozen = Vec::with_capacity(m);
let mut run = Vec::new();
loop {
if heap.is_empty() {
if !run.is_empty() {
runs.push(std::mem::take(&mut run));
}
if frozen.is_empty() && i >= n {
break;
}
heap = std::mem::take(&mut frozen);
while i < n && heap.len() < m {
heap.push(input[i]);
i += 1;
}
if heap.is_empty() {
break;
}
heapify_min(&mut heap);
continue;
}
let out = heap_pop_min(&mut heap);
run.push(out);
if i < n {
let next = input[i];
i += 1;
if next >= out {
heap_push(&mut heap, next);
} else {
frozen.push(next);
}
}
}
runs
}
fn merge_all_runs(runs: Vec<Vec<usize>>) -> Vec<usize> {
if runs.is_empty() {
return Vec::new();
}
let mut queue = runs;
while queue.len() > 1 {
let mut next = Vec::with_capacity((queue.len() + 1) / 2);
let mut idx = 0;
while idx + 1 < queue.len() {
next.push(merge_values(&queue[idx], &queue[idx + 1]));
idx += 2;
}
if idx < queue.len() {
next.push(queue[idx].clone());
}
queue = next;
}
queue.pop().unwrap_or_default()
}
fn replacement_selection_sort(a: &mut [usize]) {
let n = a.len();
if n <= 1 {
return;
}
let runs = generate_runs(a, HEAP_SIZE);
let sorted = merge_all_runs(runs);
a.copy_from_slice(&sorted);
}
fn benchmark_sort(array: &mut [usize]) {
replacement_selection_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),
);
}
// Small-input cutoffs (insertion sort below 32 elements, etc.) hide duplicate-key
// bugs in the recursive path, so repeat the duplicate cases at the smallest
// benchmark size, which every algorithm must handle within reasonable time.
check_correctness_case("all_equal_len256", vec![7; 256]);
for seed in 1..=4 {
check_correctness_case(
&format!("few_keys_len256_seed_{seed}"),
few_unique_values(256, 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 bytes (explicit buffers such as swap).
/// Kept in bytes so the parent can average before rounding; converting to KiB here
/// would truncate sub-KiB buffers to 0 in every run and hide them from the average.
fn run_once(size: usize, seed: usize) -> (u128, usize) {
let mut array = input_array(size, seed as u64);
let base_bytes = LIVE_BYTES.load(AtomicOrdering::Relaxed);
PEAK_BYTES.store(base_bytes, AtomicOrdering::Relaxed);
let start = Instant::now();
benchmark_sort(&mut array);
let elapsed = start.elapsed();
let peak_bytes = PEAK_BYTES.load(AtomicOrdering::Relaxed);
let aux_bytes = peak_bytes.saturating_sub(base_bytes);
let expected: Vec<usize> = (1..=size).collect();
if array != expected {
panic!(
"sort failed with seed {} for size {}",
seed,
size
);
}
(micros(elapsed), aux_bytes)
}
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;
// Memory is summed in bytes and converted to KiB once, after averaging.
let avg_mem_kb = total_mem / RUNS / 1024;
let max_mem_kb = max_mem / 1024;
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_kb,
max_mem_kb
);
}
}
RUST
RUN cargo build --release
CMD ["./target/release/rust-benchmark"]
EOF
docker build -t rust-benchmark "$WORKDIR"
docker run --rm --init rust-benchmark