スムースソートを使う

スムースソート (smooth sort) は、1981年に Edsger W. Dijkstra が提案した 比較ソート である。

ヒープソート (heap sort) と同様に 最悪時間計算量は O(n log n) で、追加の補助記憶域は O(1) に抑えられる インプレース なアルゴリズムだが、入力がすでに昇順に近いほど作業量が少なくなる適応的 (adaptive) 性質 を持つ。

通常の二分ヒープでは、根が配列の一端にあり、最大値を取り出すたびに 根と末尾を交換して ヒープを底側から縮めていく。この配置は「ヒープとしての整理」と「昇順に書き出す向き」が噛み合わず、すでにソート済みの入力でも一度ヒープ状に混ぜ直すため、最良でも必ず Ω(n log n) 程度の仕事が残る。

スムースソートでは、二分木ではなく レオナルド数 に基づく レオナルド木 を積み重ねた 森 を維持する。レオナルド数は Edsger W. Dijkstra による造語でフィボナッチ数列の亜種として Leonardo Fibonacci にちなんで名付けられた。

各木はヒープ条件を満たし、さらに 各木の根の値は左から右へ弱い順（非減少） として並べる。そうすると 常に右端が全体の最大 になり、末尾から確定した最大値を「取り出す」処理が すでに整っている入力ではほとんど余計な比較を要しない。

レオナルド数 L(k) は

  • L(0) = L(1) = 1
  • L(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) に近づくことから付けられている。

スムースソートの処理の概略は、単一レオナルド木での根の沈下 sift_in と、複数の木があるときに根が隣へ繰られないよう整える interheap_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 として採用されることは稀である。

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

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 {
    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 run_once(size: usize, seed: usize) -> (u128, usize) {
    let expected: Vec<usize> = (1..=size).collect();
    let mut array = shuffled(size, seed as u64);

    let start = Instant::now();

    benchmark_sort(&mut array);

    let elapsed = start.elapsed();

    if array != expected {
        panic!(
            "sort failed with seed {} for size {}",
            seed,
            size
        );
    }

    (micros(elapsed), memory_usage_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;
    }

    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 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 += mem;

            if mem > max_mem {
                max_mem = 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