ペイシェンスソートを使用する

ペイシェンスソート (patience sort) は、ソリティアでトランプを並べていく規則になぞらえたアルゴリズムで、左から並ぶ複数の山(パイル)に入力を順に載せ、その後山の一番上の値だけを読みながら昇順へ取り出していく。

  1. 積載: 入力配列を先頭から左へ見ていき、現在の値 x について、一番上の値が x より大きいなかで先頭にあるパイルを探し、その山の上へ x を載せる。

    該当する山がひとつもなければ、一番右に新しい山を増やしてそこに x だけ載せた状態から始める。

  2. 取出: すべての山が空になるまで、一番上の値のうち現在の最小値を持っている山ひとつを選び、その値を取り除いて結果列の末尾へ付ける。

    一番上の値だけを対象とするのでマージソートと同種のマージ処理の特殊形だと見なせる。

procedure patience_deal(elements)
  piles = sequence of piles, initially empty lists
  for each x in elements
    i = smallest index where piles[i] is non-empty AND top(piles[i]) > x
           (otherwise i = undefined)
    if i is undefined then
      piles.append([x])
    else
      push x onto piles[i]

procedure patience_collect(piles)
  result = []
  until all piles empty
    p = pile index minimizing top(piles[p]), breaking ties arbitrarily
    y = pop top from piles[p]
    append y to result
  return result

取出では山の総数によらず、一番上の値の参照にヒープなどを使える。積載のルールだけを切り離すと、山の個数は最長増加部分列の長さに対応することが知られている。

振り分けと取出を合わせた時間はおおむね \(O(n \log n)\) で、山数を \(k\) とすれば取出側はヒープで \(O(n \log k)\) にできる。安定ソートではなく、山と出力用に \(O(n)\) の追加領域が要る。

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

アンシャッフルソートもパイルへ載せてから整列する点は似ているが、ペイシェンスソートは山の一番上の値が新要素より厳密に大きい山にだけ載せ、積載が終わったあと各山の一番上に見えている値のうち最小を繰り返し取り出して昇順列を組み立てる。山の個数は最長増加部分列の長さに対応することが知られている。

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

Size Average time (s) Maximum time (s) Average memory (KiB) Maximum memory (KiB)
256 0.000027 0.000146 3 4
512 0.000061 0.000230 7 7
1024 0.000139 0.000654 12 15
2048 0.000324 0.000545 25 28
4096 0.000704 0.001135 48 54
8192 0.001741 0.004222 98 107
16384 0.004425 0.007648 190 206
32768 0.011365 0.017142 381 407
65536 0.029831 0.048388 756 789
131072 0.085482 0.136742 1506 1584
262144 0.257253 0.517995 3016 3094
計測に使用したコードを表示する

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 patience_sort(a: &mut [usize]) {
    let mut piles: Vec<Vec<usize>> = Vec::new();
    for &value in a.iter() {
        let mut lo = 0;
        let mut hi = piles.len();
        while lo < hi {
            let mid = (lo + hi) / 2;
            if *piles[mid].last().unwrap() >= value {
                hi = mid;
            } else {
                lo = mid + 1;
            }
        }
        if lo == piles.len() {
            piles.push(vec![value]);
        } else {
            piles[lo].push(value);
        }
    }
    for slot in a.iter_mut() {
        let mut min = 0;
        for i in 1..piles.len() {
            if piles[i].last().is_some()
                && (piles[min].last().is_none()
                    || piles[i].last().unwrap() < piles[min].last().unwrap())
            {
                min = i;
            }
        }
        *slot = piles[min].pop().unwrap();
    }
}


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

    patience_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