ウィキソートで配列を並び替える

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;


const CACHE_SIZE: usize = 512;

#[derive(Clone, Copy)]
struct Range {
    start: usize,
    end: usize,
}

impl Range {
    fn new(start: usize, end: usize) -> Self {
        Self { start, end }
    }

    fn len(self) -> usize {
        self.end - self.start
    }
}

struct WikiIterator {
    size: usize,
    power_of_two: usize,
    numerator: usize,
    decimal: usize,
    denominator: usize,
    decimal_step: usize,
    numerator_step: usize,
}

impl WikiIterator {
    fn new(size: usize, min_level: usize) -> Self {
        let power_of_two = floor_power_of_two(size);
        let denominator = power_of_two / min_level;
        Self {
            size,
            power_of_two,
            numerator: 0,
            decimal: 0,
            denominator,
            decimal_step: size / denominator,
            numerator_step: size % denominator,
        }
    }

    fn begin(&mut self) {
        self.numerator = 0;
        self.decimal = 0;
    }

    fn next_range(&mut self) -> Range {
        let start = self.decimal;
        self.decimal += self.decimal_step;
        self.numerator += self.numerator_step;
        if self.numerator >= self.denominator {
            self.numerator -= self.denominator;
            self.decimal += 1;
        }
        Range::new(start, self.decimal)
    }

    fn finished(&self) -> bool {
        self.decimal >= self.size
    }

    fn next_level(&mut self) -> bool {
        self.decimal_step += self.decimal_step;
        self.numerator_step += self.numerator_step;
        if self.numerator_step >= self.denominator {
            self.numerator_step -= self.denominator;
            self.decimal_step += 1;
        }
        self.decimal_step < self.size
    }

    fn length(&self) -> usize {
        self.decimal_step
    }
}

fn floor_power_of_two(value: usize) -> usize {
    let mut x = value;
    x |= x >> 1;
    x |= x >> 2;
    x |= x >> 4;
    x |= x >> 8;
    x |= x >> 16;
    #[cfg(target_pointer_width = "64")]
    {
        x |= x >> 32;
    }
    x - (x >> 1)
}

fn wiki_insertion_sort(a: &mut [usize], range: Range) {
    for i in range.start + 1..range.end {
        let temp = a[i];
        let mut j = i;
        while j > range.start && temp < a[j - 1] {
            a[j] = a[j - 1];
            j -= 1;
        }
        a[j] = temp;
    }
}

fn reverse(a: &mut [usize], range: Range) {
    let len = range.len();
    for index in (0..len / 2).rev() {
        a.swap(range.start + index, range.end - index - 1);
    }
}

fn rotate(a: &mut [usize], amount: usize, range: Range, cache: &mut [usize], cache_size: usize) {
    if range.len() == 0 {
        return;
    }
    let split = range.start + amount;
    let range1 = Range::new(range.start, split);
    let range2 = Range::new(split, range.end);
    if range1.len() <= range2.len() {
        if range1.len() <= cache_size {
            cache[..range1.len()].copy_from_slice(&a[range1.start..range1.end]);
            a.copy_within(range2.start..range2.end, range1.start);
            a[range1.start + range2.len()..range1.start + range2.len() + range1.len()]
                .copy_from_slice(&cache[..range1.len()]);
            return;
        }
    } else if range2.len() <= cache_size {
        cache[..range2.len()].copy_from_slice(&a[range2.start..range2.end]);
        a.copy_within(range1.start..range1.end, range2.end - range1.len());
        a[range1.start..range1.start + range2.len()].copy_from_slice(&cache[..range2.len()]);
        return;
    }
    reverse(a, range1);
    reverse(a, range2);
    reverse(a, range);
}

fn merge_into(from: &[usize], a: Range, b: Range, into: &mut [usize]) {
    let mut a_index = a.start;
    let mut b_index = b.start;
    let mut insert = 0;
    loop {
        if from[b_index] >= from[a_index] {
            into[insert] = from[a_index];
            a_index += 1;
            insert += 1;
            if a_index == a.end {
                into[insert..insert + b.end - b_index].copy_from_slice(&from[b_index..b.end]);
                break;
            }
        } else {
            into[insert] = from[b_index];
            b_index += 1;
            insert += 1;
            if b_index == b.end {
                into[insert..insert + a.end - a_index].copy_from_slice(&from[a_index..a.end]);
                break;
            }
        }
    }
}

fn merge_external(a: &mut [usize], a_range: Range, b: Range, cache: &mut [usize]) {
    cache[..a_range.len()].copy_from_slice(&a[a_range.start..a_range.end]);
    let mut a_index = 0;
    let mut b_index = b.start;
    let mut insert = a_range.start;
    let a_last = a_range.len();
    let b_last = b.end;
    if b.len() > 0 && a_range.len() > 0 {
        loop {
            if a[b_index] >= cache[a_index] {
                a[insert] = cache[a_index];
                a_index += 1;
                insert += 1;
                if a_index == a_last {
                    break;
                }
            } else {
                a[insert] = a[b_index];
                b_index += 1;
                insert += 1;
                if b_index == b_last {
                    break;
                }
            }
        }
    }
    a[insert..insert + a_last - a_index].copy_from_slice(&cache[a_index..a_last]);
}

fn merge_pair(
    a: &mut [usize],
    a_range: Range,
    b: Range,
    cache: &mut [usize],
    cache_size: usize,
) {
    if a[b.end - 1] < a[a_range.start] {
        rotate(
            a,
            a_range.len(),
            Range::new(a_range.start, b.end),
            cache,
            cache_size,
        );
    } else if a[b.start] < a[a_range.end - 1] {
        if a_range.len() + b.len() <= cache_size {
            cache[..a_range.len()].copy_from_slice(&a[a_range.start..a_range.end]);
            merge_external(a, a_range, b, cache);
        } else {
            let mut merged = Vec::with_capacity(a_range.len() + b.len());
            let (mut i, mut j) = (a_range.start, b.start);
            while i < a_range.end && j < b.end {
                if a[i] <= a[j] {
                    merged.push(a[i]);
                    i += 1;
                } else {
                    merged.push(a[j]);
                    j += 1;
                }
            }
            merged.extend_from_slice(&a[i..a_range.end]);
            merged.extend_from_slice(&a[j..b.end]);
            a[a_range.start..b.end].copy_from_slice(&merged);
        }
    }
}

fn wiki_sort(a: &mut [usize]) {
    let size = a.len();
    let mut cache = [0usize; CACHE_SIZE];
    let cache_size = CACHE_SIZE;

    if size < 4 {
        if size == 3 {
            if a[1] < a[0] {
                a.swap(0, 1);
            }
            if a[2] < a[1] {
                a.swap(1, 2);
                if a[1] < a[0] {
                    a.swap(0, 1);
                }
            }
        } else if size == 2 && a[1] < a[0] {
            a.swap(0, 1);
        }
        return;
    }

    let mut iterator = WikiIterator::new(size, 4);
    iterator.begin();
    while !iterator.finished() {
        let range = iterator.next_range();
        wiki_insertion_sort(a, range);
    }
    if size < 8 {
        return;
    }

    loop {
        if iterator.length() < cache_size {
            if (iterator.length() + 1) * 4 <= cache_size && iterator.length() * 4 <= size {
                iterator.begin();
                while !iterator.finished() {
                    let a1 = iterator.next_range();
                    let b1 = iterator.next_range();
                    let a2 = iterator.next_range();
                    let b2 = iterator.next_range();
                    let mut merged1_len = 0usize;
                    let mut merged2_len = 0usize;
                    if a[b1.end - 1] < a[a1.start] {
                        cache[b1.len()..b1.len() + a1.len()].copy_from_slice(&a[a1.start..a1.end]);
                        cache[..b1.len()].copy_from_slice(&a[b1.start..b1.end]);
                        merged1_len = a1.len() + b1.len();
                    } else if a[b1.start] < a[a1.end - 1] {
                        merge_into(a, a1, b1, &mut cache);
                        merged1_len = a1.len() + b1.len();
                    } else if !(a[b2.start] < a[a2.end - 1]) && !(a[a2.start] < a[b1.end - 1]) {
                        continue;
                    } else {
                        cache[..a1.len()].copy_from_slice(&a[a1.start..a1.end]);
                        cache[a1.len()..a1.len() + b1.len()].copy_from_slice(&a[b1.start..b1.end]);
                        merged1_len = a1.len() + b1.len();
                    }
                    let a1 = Range::new(a1.start, b1.end);
                    if a[b2.end - 1] < a[a2.start] {
                        cache[merged1_len + b2.len()..merged1_len + b2.len() + a2.len()]
                            .copy_from_slice(&a[a2.start..a2.end]);
                        cache[merged1_len..merged1_len + b2.len()].copy_from_slice(&a[b2.start..b2.end]);
                        merged2_len = a2.len() + b2.len();
                    } else if a[b2.start] < a[a2.end - 1] {
                        merge_into(a, a2, b2, &mut cache[merged1_len..]);
                        merged2_len = a2.len() + b2.len();
                    } else {
                        cache[merged1_len..merged1_len + a2.len()].copy_from_slice(&a[a2.start..a2.end]);
                        cache[merged1_len + a2.len()..merged1_len + a2.len() + b2.len()]
                            .copy_from_slice(&a[b2.start..b2.end]);
                        merged2_len = a2.len() + b2.len();
                    }
                    let a2 = Range::new(a2.start, b2.end);
                    let a3 = Range::new(0, merged1_len);
                    let b3 = Range::new(merged1_len, merged1_len + merged2_len);
                    if cache[b3.end - 1] < cache[a3.start] {
                        a[a1.start + merged2_len..a1.start + merged2_len + merged1_len]
                            .copy_from_slice(&cache[a3.start..a3.end]);
                        a[a1.start..a1.start + merged2_len].copy_from_slice(&cache[b3.start..b3.end]);
                    } else if cache[b3.start] < cache[a3.end - 1] {
                        merge_into(&cache, a3, b3, &mut a[a1.start..a1.start + merged1_len + merged2_len]);
                    } else {
                        a[a1.start..a1.start + merged1_len].copy_from_slice(&cache[a3.start..a3.end]);
                        a[a1.start + merged1_len..a1.start + merged1_len + merged2_len]
                            .copy_from_slice(&cache[b3.start..b3.end]);
                    }
                }
                iterator.next_level();
            } else {
                iterator.begin();
                while !iterator.finished() {
                    let a_range = iterator.next_range();
                    let b = iterator.next_range();
                    merge_pair(a, a_range, b, &mut cache, cache_size);
                }
            }
        } else {
            iterator.begin();
            while !iterator.finished() {
                let a_range = iterator.next_range();
                let b = iterator.next_range();
                merge_pair(a, a_range, b, &mut cache, cache_size);
            }
        }
        if !iterator.next_level() {
            break;
        }
    }
}


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

    wiki_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} | {:>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 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} | {:>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_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