ブロッククイックソートを使用する

ブロッククイックソート (BlockQuicksort) は分岐予測ミスを抑える分割を用いるクイックソートである。

ホーア分割型クイックソートと同じく左右から詰め合う分割だが、要素ごとの比較結果でその場で分岐して交換するのではなく、固定長のブロックをまとめて走査し、ずれている要素のオフセットをバッファに溜めてからまとめて交換する。

比較結果に依存する条件分岐をほぼ消し、平均の分岐予測ミスを \(\varepsilon n \log n + O(n)\)(ブロック長に依存する小さい \(\varepsilon\))程度に抑えられることが示されている。

  1. ピボットの選択: 部分配列の中央付近の要素をピボットとし、いったん末尾へ退避する。
  2. 走査(scanning): 左右に長さ \(B\)(論文では 128)のブロックを取り、ピボットと比較する。左ブロックでは \(A[i] \ge \mathrm{pivot}\) のオフセットを、右ブロックでは \(A[j] \le \mathrm{pivot}\) のオフセットをバッファへ書く。加算には比較結果を整数化した値を使い、比較ごとの分岐を避ける。
  3. 再配置(rearrangement): 両バッファから $$\min( L , R )$$ 組を取り出し、対応する要素を交換する。空になった側のブロックだけポインタを進める。
  4. 残り: 長さが \(2B\) 以下になったら、残りを同じ要領で片付け、ピボットを境界へ戻す。
  5. 再帰: ピボットより左・右を再帰する。十分短い区間は挿入ソートで仕上げる。
procedure block_quick_sort(A, lo, hi)
  if hi - lo < INSERTION_THRESHOLD then
    insertion_sort(A, lo, hi)
    return
  p = block_partition(A, lo, hi)
  block_quick_sort(A, lo, p - 1)
  block_quick_sort(A, p + 1, hi)

procedure block_partition(A, lo, hi)
  // hi is inclusive; pivot moves to A[hi]
  swap(A[lo + floor((hi - lo) / 2)], A[hi])
  pivot = A[hi]
  begin = lo
  last = hi - 1
  numL = numR = startL = startR = 0
  while last - begin + 1 > 2B
    if numL = 0 then
      startL = 0
      for i = 0 .. B - 1
        offsetsL[numL] = i
        numL = numL + (A[begin + i] >= pivot)   // branchless increment
    if numR = 0 then
      startR = 0
      for i = 0 .. B - 1
        offsetsR[numR] = i
        numR = numR + (A[last - i] <= pivot)
    num = min(numL, numR)
    for j = 0 .. num - 1
      swap(A[begin + offsetsL[startL + j]], A[last - offsetsR[startR + j]])
    numL = numL - num; numR = numR - num
    startL = startL + num; startR = startR + num
    if numL = 0 then begin = begin + B
    if numR = 0 then last = last - B
  // finish remaining ≤ 2B elements (same buffers), then place pivot
  ...
  return pivot_index

平均計算量は通常のクイックソートと同様 \(O(n \log n)\) で、ピボットが偏ると最悪 \(O(n^2)\) になりうる。追加メモリはオフセット用の \(O(B)\) と再帰スタック程度で、実質インプレースである。不安定である。

デモではブロック長を \(B = 4\) に下げ、走査と再配置が見えるようにしている(計測コードは \(B = 128\))。

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

ホーア分割は左右ポインタを 1 要素ずつ進め、比較のたびに「進む/止めて交換」の分岐が起きる。ブロック分割は走査と交換を分離し、比較結果をオフセットバッファへ畳み込む。

ロムート分割は片方向の境界更新が中心で、やはり比較ごとの分岐が多い。パターン撃退型クイックソートは偏った入力や等値だらけへの耐性を足すハイブリッドであり、分岐予測そのものへの対策ではない。

デュアルピボットサンプルソートはピボット数やバケツ分けでスキャン効率を上げる系統で、BlockQuicksort の「定数サイズバッファによる分岐削減」とは直交する改良である。

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

Size Average time (s) Maximum time (s) Average memory (KiB) Maximum memory (KiB)
256 0.000004 0.000081 0 0
512 0.000008 0.000041 0 0
1024 0.000017 0.000069 0 0
2048 0.000034 0.000078 0 0
4096 0.000070 0.000131 0 0
8192 0.000144 0.000461 0 0
16384 0.000295 0.000438 0 0
32768 0.000612 0.001124 0 0
65536 0.001264 0.001921 0 0
131072 0.002624 0.006469 0 0
262144 0.005401 0.011292 0 0
計測に使用したコードを表示する

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 insertion_sort(a: &mut [usize]) {
    for i in 1..a.len() {
        let mut j = i;
        while j > 0 && a[j - 1] > a[j] {
            a.swap(j - 1, j);
            j -= 1;
        }
    }
}



const BLOCK_SIZE: usize = 128;
const INSERTION_THRESHOLD: usize = 16;

/// Branch-light Hoare-style partition used by BlockQuicksort (Edelkamp & Weiß).
/// Returns the final index of the pivot.
fn quick_block_partition(a: &mut [usize], mut begin: usize, end: usize) -> usize {
    // end is exclusive; pivot starts at midpoint of [begin, end).
    let mid = begin + (end - begin) / 2;
    a.swap(mid, end - 1);
    let pivot = a[end - 1];
    let mut last = end - 2;

    let mut index_l = [0usize; BLOCK_SIZE];
    let mut index_r = [0usize; BLOCK_SIZE];
    let mut num_left = 0usize;
    let mut num_right = 0usize;
    let mut start_left = 0usize;
    let mut start_right = 0usize;

    while begin <= last && last - begin + 1 > 2 * BLOCK_SIZE {
        if num_left == 0 {
            start_left = 0;
            for j in 0..BLOCK_SIZE {
                index_l[num_left] = j;
                // left buffer: elements >= pivot (need to move right)
                num_left += usize::from(!(a[begin + j] < pivot));
            }
        }
        if num_right == 0 {
            start_right = 0;
            for j in 0..BLOCK_SIZE {
                index_r[num_right] = j;
                // right buffer: elements <= pivot (need to move left)
                num_right += usize::from(!(pivot < a[last - j]));
            }
        }

        let num = num_left.min(num_right);
        for j in 0..num {
            let li = begin + index_l[start_left + j];
            let ri = last - index_r[start_right + j];
            a.swap(li, ri);
        }
        num_left -= num;
        num_right -= num;
        start_left += num;
        start_right += num;
        if num_left == 0 {
            begin += BLOCK_SIZE;
        }
        if num_right == 0 {
            last -= BLOCK_SIZE;
        }
    }

    // Final (partial) scan of the remaining ≤ 2B elements.
    let (shift_l, shift_r) = if num_right == 0 && num_left == 0 {
        debug_assert!(begin <= last);
        let len = last - begin + 1;
        let shift_l = len / 2;
        let shift_r = len - shift_l;
        start_left = 0;
        start_right = 0;
        for j in 0..shift_l {
            index_l[num_left] = j;
            num_left += usize::from(!(a[begin + j] < pivot));
            index_r[num_right] = j;
            num_right += usize::from(!(pivot < a[last - j]));
        }
        if shift_l < shift_r {
            index_r[num_right] = shift_r - 1;
            num_right += usize::from(!(pivot < a[last - (shift_r - 1)]));
        }
        (shift_l, shift_r)
    } else if num_right != 0 {
        let shift_l = last - begin + 1 - BLOCK_SIZE;
        start_left = 0;
        for j in 0..shift_l {
            index_l[num_left] = j;
            num_left += usize::from(!(a[begin + j] < pivot));
        }
        (shift_l, BLOCK_SIZE)
    } else {
        let shift_r = last - begin + 1 - BLOCK_SIZE;
        start_right = 0;
        for j in 0..shift_r {
            index_r[num_right] = j;
            num_right += usize::from(!(pivot < a[last - j]));
        }
        (BLOCK_SIZE, shift_r)
    };

    let num = num_left.min(num_right);
    for j in 0..num {
        let li = begin + index_l[start_left + j];
        let ri = last - index_r[start_right + j];
        a.swap(li, ri);
    }
    num_left -= num;
    num_right -= num;
    start_left += num;
    start_right += num;
    if num_left == 0 {
        begin += shift_l;
    }
    if num_right == 0 {
        last = last.wrapping_sub(shift_r);
    }

    // Drain leftovers still recorded in one buffer.
    // `upper` is signed because the reference finish may leave it at -1.
    if num_left != 0 {
        let mut lower_i = (start_left + num_left - 1) as isize;
        let mut upper = last as isize - begin as isize;
        while lower_i >= start_left as isize && index_l[lower_i as usize] as isize == upper {
            upper -= 1;
            lower_i -= 1;
        }
        while lower_i >= start_left as isize {
            a.swap(
                (begin as isize + upper) as usize,
                begin + index_l[lower_i as usize],
            );
            upper -= 1;
            lower_i -= 1;
        }
        let pivot_pos = (begin as isize + upper + 1) as usize;
        a.swap(end - 1, pivot_pos);
        pivot_pos
    } else if num_right != 0 {
        let mut lower_i = (start_right + num_right - 1) as isize;
        let mut upper = last as isize - begin as isize;
        while lower_i >= start_right as isize && index_r[lower_i as usize] as isize == upper {
            upper -= 1;
            lower_i -= 1;
        }
        while lower_i >= start_right as isize {
            a.swap(
                (last as isize - upper) as usize,
                (last as isize - index_r[lower_i as usize] as isize) as usize,
            );
            upper -= 1;
            lower_i -= 1;
        }
        let pivot_pos = (last as isize - upper) as usize;
        a.swap(end - 1, pivot_pos);
        pivot_pos
    } else {
        a.swap(end - 1, begin);
        begin
    }
}

fn quick_block_sort_range(a: &mut [usize], lo: usize, hi: usize) {
    if hi <= lo {
        return;
    }
    if hi - lo < INSERTION_THRESHOLD {
        insertion_sort(&mut a[lo..=hi]);
        return;
    }
    let pivot_pos = quick_block_partition(a, lo, hi + 1);
    if pivot_pos > lo {
        quick_block_sort_range(a, lo, pivot_pos - 1);
    }
    if pivot_pos < hi {
        quick_block_sort_range(a, pivot_pos + 1, hi);
    }
}

fn quick_block_sort(a: &mut [usize]) {
    if let Some(hi) = a.len().checked_sub(1) {
        quick_block_sort_range(a, 0, hi);
    }
}


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

    quick_block_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);
    }
}

// Skip cases larger than the algorithm's measured size cap (MAX_POWER). That
// cap exists because larger inputs are impractically slow; forcing them here
// would stall the published measurement script before any table rows print.
fn check_correctness_case_within_limit(label: &str, input: Vec<usize>) {
    if input.len() > (1usize << MAX_POWER) {
        return;
    }
    check_correctness_case(label, input);
}

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),
        );
    }
    // Blit's equal-key second sweep used to copy the whole range into a fixed
    // 512-element swap; lengths above that must still sort without panicking.
    // Respect MAX_POWER so algorithms with a low measured-size cap (slow,
    // sleep) do not hang here for minutes or months.
    check_correctness_case_within_limit("all_equal_len600", vec![7; 600]);
    for seed in 1..=4 {
        check_correctness_case_within_limit(
            &format!("few_keys_len2048_seed_{seed}"),
            few_unique_values(2048, 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} | {:>16} | {:>16} | {:>20} | {:>20} |",
        "Size",
        "Average time (s)",
        "Maximum time (s)",
        "Average memory (KiB)",
        "Maximum memory (KiB)"
    );

    println!(
        "|{:-<11}:|{:-<17}:|{:-<17}:|{:-<21}:|{:-<21}:|",
        "",
        "",
        "",
        "",
        ""
    );

    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} | {:>16} | {:>16} | {:>20} | {:>20} |",
            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