フランチェスキーニソートを使用する

フランチェスキーニソート (Franceschini sort) は比較回数・要素移動・補助記憶を同時に漸近最適へ近づけるインプレース整列の系統である。

古典的な未解決問題——最悪でも比較 O(n log n)・移動 O(n)・補助記憶 O(1) を両立できるか——に対し、肯定的な構成を与えたことで知られる。

安定版は同値の相対順序も保ちつつ同じ資源境界を狙う。実用上は定数倍と実装の重さからウィキソートやグレイルソートほどは使われないが、理論上の到達点として重要である。

下記のデモと計測コードは、フランチェスキーニソートが提案された論文の骨格を次のように簡略化した版である。

  1. バッファの切り出し: 順位おおよそ n/4 の要素をピボットにし、厳密に小さい要素を先頭へ集める。左側(アクティブ)は約 n/4、右側(バッファ)は約 3n/4 になる。
  2. バッファ付き部分整列: アクティブ区間をバッファ先頭と交換し、そこを高い分岐数の d 分木ヒープソートで整える。分岐数をおよそ n^{1/4} に取るとヒープの高さが定数に近く、要素あたりの移動が抑えられる。整列後、再び交換してアクティブ位置へ戻す。
  3. 残りへの再帰: ピボット以上の未整列側へ同じ手順を繰り返す。左はすでに整っており、かつ右のどの要素より小さいので、連結した配列全体が昇順になる。
  4. 小さな入力: 長さが小さいときは挿入ソート、または同じ d 分木ヒープへフォールバックする。

論文本体では、さらに標本とセグメント構造・ビット符号化(最小/最大要素ブロックの交換でポインタビットを作る)などで移動回数を O(n) に押し込む。計測コードはその外側の「四分割+バッファ+高分岐ヒープ」までを実装している。

procedure dary_heap_sort(A)
  d = roughly length(A)^(1/4)
  build_max_heap_with_branching_d(A)
  for end from length(A)-1 down to 1
    swap A[0] with A[end]
    sift_down(A, 0, end-1, d)

procedure sort_with_buffer(Active, Buffer)
  // |Buffer| >= |Active|
  swap Active with Buffer[0 .. |Active|)
  dary_heap_sort(Buffer[0 .. |Active|))
  swap back

procedure franceschini_sort(A)
  n = length(A)
  if n is small then
    insertion_or_dary_heap_sort(A)
    return
  pivot = select_kth(A, floor(n/4))
  split = stable_gather of elements strictly < pivot to the front
  if split = 0 or split > n - split then
    dary_heap_sort(A)
    return
  sort_with_buffer(A[0 .. split), A[split .. n))
  franceschini_sort(A[split .. n))

デモでは小さな配列向けに分岐数と閾値を下げている。本番の計測コードはより大きい入力で同じ骨格を動かす。

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

ウィキソートグレイルソートコタソートも原地安定な O(n log n) を狙うブロックマージ系だが、内部バッファやキータグで隣接ランを併合する。フランチェスキーニソートは順位分割でバッファ領域を切り出し、高分岐ヒープなどで移動回数そのものを漸近的に減らす点が異なる。

ヒープソートの二分ヒープは移動が Θ(n log n) になりやすい。こちらは分岐数を大きくして高さを抑え、論文の「移動 O(n)」側の直感に寄せている。

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

Size Average time Maximum time Average memory Maximum memory
256 0.000009 0.000053 70 76
512 0.000020 0.000071 58 64
1024 0.000040 0.000103 66 72
2048 0.000091 0.000353 74 80
4096 0.000185 0.000335 62 68
8192 0.000426 0.004516 74 80
16384 0.000959 0.005731 62 68
32768 0.002238 0.014780 58 64
65536 0.005032 0.013629 78 84
131072 0.011317 0.035816 74 80
262144 0.025697 0.048032 57 64
計測に使用したコードを表示する

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;


// Pedagogical Franceschini–Geffert skeleton:
// peel ~n/4 elements by rank, sort that prefix with the remaining suffix as a
// swap buffer via a high-arity heap (constant-ish height → few moves per element),
// then recurse on the unsorted suffix. The published algorithm adds sample /
// segment / bit-encoding machinery for full O(n)-move asymptotics; this port
// keeps the outer structure and the d-ary heap idea for measurement.

fn franceschini_branch_factor(len: usize) -> usize {
    if len <= 2 {
        return 2;
    }
    // Aim for heap height around 4: d ≈ n^(1/4).
    let mut d = 2usize;
    while d * d * d * d < len {
        d += 1;
        if d > 64 {
            break;
        }
    }
    d.max(2)
}

fn franceschini_child(parent: usize, which: usize, d: usize) -> usize {
    parent * d + 1 + which
}

fn franceschini_sift_down(a: &mut [usize], mut root: usize, end: usize, d: usize) {
    loop {
        let first = franceschini_child(root, 0, d);
        if first > end {
            break;
        }
        let mut best = first;
        let last = (first + d - 1).min(end);
        for child in first + 1..=last {
            if a[child] > a[best] {
                best = child;
            }
        }
        if a[root] >= a[best] {
            break;
        }
        a.swap(root, best);
        root = best;
    }
}

fn franceschini_dary_heap_sort(a: &mut [usize]) {
    let n = a.len();
    if n <= 1 {
        return;
    }
    let d = franceschini_branch_factor(n);
    let last_parent = (n - 2) / d;
    for start in (0..=last_parent).rev() {
        franceschini_sift_down(a, start, n - 1, d);
    }
    for end in (1..n).rev() {
        a.swap(0, end);
        if end > 1 {
            franceschini_sift_down(a, 0, end - 1, d);
        }
    }
}

fn franceschini_insertion_sort(a: &mut [usize]) {
    for i in 1..a.len() {
        let key = a[i];
        let mut j = i;
        while j > 0 && a[j - 1] > key {
            a[j] = a[j - 1];
            j -= 1;
        }
        a[j] = key;
    }
}

fn franceschini_partition_at(a: &mut [usize], left: usize, right: usize, pivot_index: usize) -> usize {
    a.swap(pivot_index, right);
    let pivot = a[right];
    let mut store = left;
    for i in left..right {
        if a[i] < pivot {
            a.swap(store, i);
            store += 1;
        }
    }
    a.swap(store, right);
    store
}

fn franceschini_quickselect(a: &mut [usize], mut left: usize, mut right: usize, k: usize) {
    while left < right {
        let mid = left + (right - left) / 2;
        // Median-of-three pivot index.
        if a[right] < a[left] {
            a.swap(left, right);
        }
        if a[mid] < a[left] {
            a.swap(left, mid);
        }
        if a[right] < a[mid] {
            a.swap(mid, right);
        }
        let pivot_index = franceschini_partition_at(a, left, right, mid);
        if k == pivot_index {
            return;
        } else if k < pivot_index {
            if pivot_index == 0 {
                return;
            }
            right = pivot_index - 1;
        } else {
            left = pivot_index + 1;
        }
    }
}

fn franceschini_sort_with_buffer(active: &mut [usize], buffer: &mut [usize]) {
    let m = active.len();
    if m == 0 {
        return;
    }
    debug_assert!(buffer.len() >= m);
    for i in 0..m {
        std::mem::swap(&mut active[i], &mut buffer[i]);
    }
    if m <= 32 {
        franceschini_insertion_sort(&mut buffer[..m]);
    } else {
        franceschini_dary_heap_sort(&mut buffer[..m]);
    }
    for i in 0..m {
        std::mem::swap(&mut active[i], &mut buffer[i]);
    }
}

fn franceschini_rec(a: &mut [usize]) {
    let n = a.len();
    if n <= 1 {
        return;
    }
    if n <= 64 {
        if n <= 32 {
            franceschini_insertion_sort(a);
        } else {
            franceschini_dary_heap_sort(a);
        }
        return;
    }

    let rank = n / 4;
    franceschini_quickselect(a, 0, n - 1, rank);
    let pivot = a[rank];

    let mut split = 0usize;
    for i in 0..n {
        if a[i] < pivot {
            a.swap(split, i);
            split += 1;
        }
    }

    // Need a non-empty active prefix and a buffer at least as large.
    if split == 0 || split > n - split {
        franceschini_dary_heap_sort(a);
        return;
    }

    {
        let (left, right) = a.split_at_mut(split);
        franceschini_sort_with_buffer(left, right);
    }

    franceschini_rec(&mut a[split..]);
}

fn franceschini_sort(a: &mut [usize]) {
    franceschini_rec(a);
}


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

    franceschini_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