ファンエンデボアスソートを使用する

ファンエンデボアスソート (van emde boas sort) は、整数宇宙 [0, U) 上のファン・エムデ・ボアス木(vEB 木)へキーを挿入し、最小値から順に後続(successor)をたどって取り出す非比較ソートである。

vEB 木は最小値・最大値を定数時間で返し、挿入・後続探索を O(log log U) で行う。そのため n 個のキーを整列する全体の時間は O(n log log U) になる。宇宙サイズ U が入力長 n に近い整数データでは、比較ソートの Ω(n log n) より漸近的に有利になりうる。

構造の要点は、宇宙を高位桁と低位桁に再帰的に分割することである。U = 2^{2k}(またはそれに近い 2 の冪)のとき、各ノードは次を持つ。

  • min / max: その部分宇宙に含まれる最小・最大キー(木全体から切り離して保持する)。
  • cluster: 低位桁用の部分木を √U 本(実際には上側平方根本)。キー x の高位 high(x) がクラスタ番号、低位 low(x) がクラスタ内の位置になる。
  • summary: 「どのクラスタが空でないか」を表す、宇宙サイズ √U の vEB 木。

空でないクラスタだけを遅延確保すれば、疎なキー集合でも全宇宙分の配列を一気に確保しなくてよい。重複キーは出現回数を別配列で数え、vEB 木にはユニークなオフセットだけを入れる。

  1. 値域の正規化: 最小値 min を引き、オフセット 0 … max-min へ写す。宇宙サイズ U は値域幅以上の最小の 2 の冪(ただし 2 以上)とする。
  2. 集計と挿入: 各オフセットの出現回数を数え、回数が正のキーだけを vEB 木へ挿入する。
  3. 昇順取り出し: 木の最小値から始め、successor で次のキーへ進みながら、出現回数ぶん出力配列へ書き戻す。
procedure veb_insert(V, x)
  if V.min = NIL then
    V.min = V.max = x; return
  if x < V.min then swap x with V.min
  if V.u > 2 then
    h = high(x); l = low(x)
    if V.cluster[h] is empty then
      veb_insert(V.summary, h)
      V.cluster[h].min = V.cluster[h].max = l
    else
      veb_insert(V.cluster[h], l)
  if x > V.max then V.max = x

procedure veb_successor(V, x)
  if V.u = 2 then
    if x = 0 and V.max = 1 then return 1 else return NIL
  if V.min ≠ NIL and x < V.min then return V.min
  h = high(x); l = low(x)
  if low-part of cluster h has a key > l then
    return index(h, veb_successor(V.cluster[h], l))
  succ = veb_successor(V.summary, h)
  if succ = NIL then return NIL
  return index(succ, V.cluster[succ].min)

procedure van_emde_boas_sort(A)
  if length(A) ≤ 1 then return
  minVal = minimum(A); maxVal = maximum(A)
  span = maxVal - minVal + 1
  count[0..span-1] = 0
  for each x in A
    count[x - minVal] = count[x - minVal] + 1
  U = next_power_of_two(max(span, 2))
  V = empty vEB tree with universe U
  for v from 0 to span - 1
    if count[v] > 0 then veb_insert(V, v)
  idx = 0; cur = V.min
  while cur ≠ NIL
    repeat count[cur] times
      A[idx] = minVal + cur; idx = idx + 1
    cur = veb_successor(V, cur)

キー同士の大小比較は行わず、ビット分割と再帰的な summary 操作で順序を決める。同値は集計配列側でまとめて出力するため、入力を左から数えた実装では安定ソートになる。一方で補助構造は宇宙サイズに依存し、U が極端に大きいとメモリと定数倍が膨らむ。

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

カウンティングソート鳩の巣ソート は値域 k に対し O(n + k) で直接バケットを走査する。本アルゴリズムはバケット配列を線形走査する代わりに、vEB 木の O(log log U) 操作で「次に小さいキー」だけをたどる。ツリーソート が比較で二分探索木を育てるのに対し、 ここでは宇宙のビット分割が順序を決める。トライソート の桁トライとも「桁で空間を割る」点は近いが、 summary による空クラスタのスキップが vEB 木の特徴である。

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

Size Average time Maximum time Average memory Maximum memory
256 0.000025 0.000256 98 104
512 0.000046 0.000236 141 148
1024 0.000089 0.000305 198 204
2048 0.000179 0.000457 334 340
4096 0.000350 0.006436 550 556
8192 0.000816 0.002562 968 1004
16384 0.001450 0.003449 1923 1944
32768 0.003111 0.007243 3820 3860
65536 0.006850 0.019965 8168 8208
131072 0.014069 0.026511 16488 16528
262144 0.033484 0.134041 33003 33044
計測に使用したコードを表示する

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;


struct VebTree {
    universe: usize,
    min: Option<usize>,
    max: Option<usize>,
    summary: Option<Box<VebTree>>,
    cluster: Vec<Option<Box<VebTree>>>,
}

fn lower_sqrt(universe: usize) -> usize {
    1usize << (universe.trailing_zeros() / 2)
}

impl VebTree {
    fn new(universe: usize) -> Self {
        debug_assert!(universe.is_power_of_two());
        debug_assert!(universe >= 2);

        if universe == 2 {
            return Self {
                universe,
                min: None,
                max: None,
                summary: None,
                cluster: Vec::new(),
            };
        }

        let lower = lower_sqrt(universe);
        let upper = universe / lower;

        Self {
            universe,
            min: None,
            max: None,
            summary: Some(Box::new(Self::new(upper))),
            cluster: (0..upper).map(|_| None).collect(),
        }
    }

    fn high(&self, x: usize) -> usize {
        x / lower_sqrt(self.universe)
    }

    fn low(&self, x: usize) -> usize {
        x % lower_sqrt(self.universe)
    }

    fn index(&self, high: usize, low: usize) -> usize {
        high * lower_sqrt(self.universe) + low
    }

    fn minimum(&self) -> Option<usize> {
        self.min
    }

    fn maximum(&self) -> Option<usize> {
        self.max
    }

    fn cluster_mut(&mut self, i: usize) -> &mut VebTree {
        let lower = lower_sqrt(self.universe);
        self.cluster[i].get_or_insert_with(|| Box::new(Self::new(lower)))
    }

    fn empty_insert(&mut self, x: usize) {
        self.min = Some(x);
        self.max = Some(x);
    }

    fn insert(&mut self, mut x: usize) {
        if self.min.is_none() {
            self.empty_insert(x);
            return;
        }

        if x < self.min.unwrap() {
            let old_min = self.min.unwrap();
            self.min = Some(x);
            x = old_min;
        }

        if self.universe > 2 {
            let h = self.high(x);
            let l = self.low(x);
            if self.cluster[h]
                .as_ref()
                .and_then(|c| c.minimum())
                .is_none()
            {
                self.summary.as_mut().unwrap().insert(h);
                self.cluster_mut(h).empty_insert(l);
            } else {
                self.cluster_mut(h).insert(l);
            }
        }

        if x > self.max.unwrap() {
            self.max = Some(x);
        }
    }

    fn successor(&self, x: usize) -> Option<usize> {
        if self.universe == 2 {
            return if x == 0 && self.max == Some(1) {
                Some(1)
            } else {
                None
            };
        }

        if let Some(min) = self.min {
            if x < min {
                return Some(min);
            }
        } else {
            return None;
        }

        let h = self.high(x);
        let l = self.low(x);
        let max_low = self.cluster[h].as_ref().and_then(|c| c.maximum());
        if max_low.is_some_and(|m| l < m) {
            let offset = self.cluster[h].as_ref().unwrap().successor(l).unwrap();
            return Some(self.index(h, offset));
        }

        let succ_cluster = self.summary.as_ref().unwrap().successor(h)?;
        let offset = self.cluster[succ_cluster]
            .as_ref()
            .and_then(|c| c.minimum())
            .unwrap();
        Some(self.index(succ_cluster, offset))
    }
}

fn van_emde_boas_sort(a: &mut [usize]) {
    if a.len() <= 1 {
        return;
    }

    let min = *a.iter().min().unwrap();
    let max = *a.iter().max().unwrap();
    let span = max - min + 1;
    let mut count = vec![0usize; span];

    for &x in a.iter() {
        count[x - min] += 1;
    }

    let universe = span.next_power_of_two().max(2);
    let mut tree = VebTree::new(universe);

    for (offset, &c) in count.iter().enumerate() {
        if c > 0 {
            tree.insert(offset);
        }
    }

    let mut idx = 0;
    let mut cur = tree.minimum();
    while let Some(v) = cur {
        let value = min + v;
        for _ in 0..count[v] {
            a[idx] = value;
            idx += 1;
        }
        cur = tree.successor(v);
    }
}


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

    van_emde_boas_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],
    );
    for seed in 0..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 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;
    }

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