ペアリングヒープソートを使用する

ペアリングヒープソート (pairing heap sort) は、要素をペアリングヒープへ挿入したあと、最小値を繰り返し取り出して昇順にする整列である。

ペアリングヒープは、根が最小（または最大）となる多分岐木である。子は左子・右兄弟形式で並べ、ヒープ条件は「親のキーはどの子のキー以下」だけを課す。二項ヒープのように次数を揃えたり、フィボナッチヒープのようにランクを管理したりしない代わりに、合併 と 子のペアリング で構造を保つ。

  1. 合併: 2 本のヒープの根を比較し、キーの大きい方を小さい方の最左の子にする。比較は 1 回で、償却的に \(O(1)\)。
  2. 挿入: 単一節点のヒープを既存ヒープと合併する。
  3. 抽出: 根を外し、その子たちを左から 2 本ずつ合併（第 1 パス）したあと、できたヒープを右から順に合併（第 2 パス）して新しい根を得る。償却的に \(O(\log n)\)。
  4. 書き戻し: 取り出したキーを配列の先頭から順に書けば昇順になる。
procedure meld(H1, H2)
  if H1 is empty then
    return H2
  if H2 is empty then
    return H1
  if H1.key <= H2.key then
    make H2 the leftmost child of H1
    return H1
  else
    make H1 the leftmost child of H2
    return H2

procedure two_pass_meld(children)   // children は兄弟リスト
  // 第1パス: 左から隣り合う 2 本を合併
  pairs = empty list
  while children is not empty
    a = take first child
    if children is empty then
      append a to pairs
    else
      b = take first child
      append meld(a, b) to pairs
  // 第2パス: 右から順に合併
  H = empty
  for p in reverse(pairs)
    H = meld(p, H)
  return H

procedure pairing_heap_sort(A)
  H = empty pairing heap
  for x in A
    H = meld(H, singleton(x))
  for i from 0 to length(A) - 1
    A[i] = H.key
    H = two_pass_meld(children of H)

償却時間計算量は挿入・合併が \(O(1)\)、最小抽出が \(O(\log n)\) であり、全体では \(O(n \log n)\) になる。節点用に \(O(n)\) の追加記憶域が要る（インプレースではない）。等値キーの相対順序は合併時の規約に依存し、不安定である。実装が単純な一方で、ポインタ経由の多分岐木はキャッシュ効率では配列上のヒープソートに劣りやすい。

優先度付きキューとしてのペアリングヒープは、合併が 1 比較で書け、減少キーも「切り離して再合併」で扱える点が実務寄りの実装で選ばれやすい。整列用途ではその操作を「すべて挿入してからすべて取り出す」形に固定したものがペアリングヒープソートである。

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

ヒープソートは配列上の二分ヒープをインプレースで縮める。ペアリングヒープはポインタの多分岐木で、合併と子のペアリングが中心になる。

二項ヒープソートは次数の異なる二項木を二進加算のように結合する。ペアリングヒープは次数を持たず、抽出時に子を 2 本ずつペアにして畳み込む。

弱ヒープソートは配列上の不完全木と逆ビットで比較回数を抑える。ヒープ同士の合併を第一級には扱わない。

トーナメントソートも最小を繰り返し取り出すが、固定長のトーナメント木を更新する点が異なる。

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

Size Average time (s) Maximum time (s) Average memory (KiB) Maximum memory (KiB)
256 0.000028 0.000799 6 7
512 0.000059 0.000191 13 14
1024 0.000137 0.000387 27 28
2048 0.000265 0.000996 55 56
4096 0.000541 0.001161 110 112
8192 0.001122 0.002107 221 224
16384 0.002600 0.004976 442 448
32768 0.005757 0.014348 885 896
65536 0.013223 0.049456 1771 1792
131072 0.028593 0.115439 3540 3584
262144 0.075095 0.538889 7083 7168
計測に使用したコードを表示する 

#!/usr/bin/env swift
import Foundation

// This standalone Swift driver creates the same temporary Docker build
// context as the former shell wrapper.  The benchmark program itself remains
// embedded below so readers can copy one complete, reproducible file.
struct BenchmarkError: Error, CustomStringConvertible {
    let message: String

    var description: String { message }

    init(_ message: String) {
        self.message = message
    }
}

func runCommand(_ executable: String, _ arguments: [String]) throws {
    let process = Process()
    process.executableURL = URL(fileURLWithPath: "/usr/bin/env")
    process.arguments = [executable] + arguments
    process.standardInput = FileHandle.standardInput
    process.standardOutput = FileHandle.standardOutput
    process.standardError = FileHandle.standardError

    do {
        try process.run()
    } catch {
        throw BenchmarkError("Could not start \(executable): \(error)")
    }
    process.waitUntilExit()
    guard process.terminationStatus == 0 else {
        throw BenchmarkError(
            "Command failed (\(process.terminationStatus)): " +
            "\(executable) \(arguments.joined(separator: " "))"
        )
    }
}

do {
    // The UUID avoids collisions when two benchmark copies are run at once.
    let workdir = FileManager.default.temporaryDirectory
        .appendingPathComponent("swift-sort-benchmark-\(UUID().uuidString)")
    try FileManager.default.createDirectory(at: workdir, withIntermediateDirectories: true)
    defer { try? FileManager.default.removeItem(at: workdir) }

    // A raw Swift string is used so the nested main.swift keeps its own
    // interpolation expressions such as \(seed) until Docker compiles it.
    let dockerfile = #"""
FROM swift:6.0

WORKDIR /app

RUN cat > alloc_track.c <<'ALLOC'
#define _GNU_SOURCE
#include <dlfcn.h>
#include <malloc.h>
#include <stdatomic.h>
#include <stddef.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>

static atomic_size_t live_bytes = 0;
static atomic_size_t peak_bytes = 0;

static void *(*real_malloc)(size_t) = NULL;
static void *(*real_calloc)(size_t, size_t) = NULL;
static void *(*real_realloc)(void *, size_t) = NULL;
static void (*real_free)(void *) = NULL;

static void init_reals(void) {
    if (real_malloc) {
        return;
    }
    real_malloc = (void *(*)(size_t))dlsym(RTLD_NEXT, "malloc");
    real_calloc = (void *(*)(size_t, size_t))dlsym(RTLD_NEXT, "calloc");
    real_realloc = (void *(*)(void *, size_t))dlsym(RTLD_NEXT, "realloc");
    real_free = (void (*)(void *))dlsym(RTLD_NEXT, "free");
}

static void record_alloc(size_t size) {
    size_t live = atomic_fetch_add(&live_bytes, size) + size;
    size_t peak = atomic_load(&peak_bytes);
    while (live > peak) {
        if (atomic_compare_exchange_weak(&peak_bytes, &peak, live)) {
            break;
        }
    }
}

void alloc_track_reset_peak(void) {
    atomic_store(&peak_bytes, atomic_load(&live_bytes));
}

size_t alloc_track_live(void) { return atomic_load(&live_bytes); }
size_t alloc_track_peak(void) { return atomic_load(&peak_bytes); }

void *malloc(size_t size) {
    init_reals();
    void *p = real_malloc(size);
    if (p) {
        record_alloc(malloc_usable_size(p));
    }
    return p;
}

void *calloc(size_t nmemb, size_t size) {
    init_reals();
    void *p = real_calloc(nmemb, size);
    if (p) {
        record_alloc(malloc_usable_size(p));
    }
    return p;
}

void *realloc(void *ptr, size_t size) {
    init_reals();
    size_t old_size = 0;
    if (ptr) {
        old_size = malloc_usable_size(ptr);
    }
    void *p = real_realloc(ptr, size);
    if (p) {
        atomic_fetch_sub(&live_bytes, old_size);
        record_alloc(malloc_usable_size(p));
    } else if (size == 0) {
        atomic_fetch_sub(&live_bytes, old_size);
    }
    return p;
}

void free(void *ptr) {
    init_reals();
    if (ptr) {
        atomic_fetch_sub(&live_bytes, malloc_usable_size(ptr));
        real_free(ptr);
    }
}

ALLOC

RUN cat > main.swift <<'SWIFT'
import Foundation
#if canImport(Glibc)
import Glibc
#elseif canImport(Darwin)
import Darwin
#endif

@_silgen_name("alloc_track_live") func alloc_track_live() -> Int
@_silgen_name("alloc_track_peak") func alloc_track_peak() -> Int
@_silgen_name("alloc_track_reset_peak") func alloc_track_reset_peak()

extension UnsafeMutableBufferPointer where Element == Int {
    func swapAt(_ i: Int, _ j: Int) {
        let t = self[i]; self[i] = self[j]; self[j] = t
    }
}

let MIN_POWER: Int = 8
let MAX_POWER: Int = 18
let RUNS: Int = 8192


fileprivate final class PairingNode {
    var key: Int
    var child: PairingNode?
    var sibling: PairingNode?

    init(key: Int, child: PairingNode? = nil, sibling: PairingNode? = nil) {
        self.key = key
        self.child = child
        self.sibling = sibling
    }
}

fileprivate func pairing_meld(_ a: PairingNode?, _ b: PairingNode?) -> PairingNode? {
    guard let a else { return b }
    guard let b else { return a }
    if a.key <= b.key {
        b.sibling = a.child
        a.child = b
        return a
    } else {
        a.sibling = b.child
        b.child = a
        return b
    }
}

fileprivate func pairing_two_pass_meld(_ first: PairingNode?) -> PairingNode? {
    var first = first
    var pairs = [PairingNode?]()
    while let a = first {
        first = a.sibling
        a.sibling = nil
        if let b = first {
            first = b.sibling
            b.sibling = nil
            pairs.append(pairing_meld(a, b))
        } else {
            pairs.append(a)
        }
    }

    var result: PairingNode? = nil
    for pair in pairs.reversed() {
        result = pairing_meld(pair, result)
    }
    return result
}

fileprivate func pairing_insert_key(_ heap: PairingNode?, _ key: Int) -> PairingNode? {
    let node = PairingNode(key: key)
    return pairing_meld(heap, node)
}

fileprivate func pairing_extract_min(_ heap: PairingNode?) -> (Int?, PairingNode?) {
    guard let root = heap else {
        return (nil, nil)
    }
    let key = root.key
    let children = root.child
    root.child = nil
    return (key, pairing_two_pass_meld(children))
}

func pairing_heap_sort(_ a: inout [Int]) {
    a.withUnsafeMutableBufferPointer { pairing_heap_sort($0) }
}

func pairing_heap_sort(_ a: UnsafeMutableBufferPointer<Int>) {
    var heap: PairingNode? = nil
    for i in 0..<a.count {
        heap = pairing_insert_key(heap, a[i])
    }
    for i in 0..<a.count {
        let (key, next) = pairing_extract_min(heap)
        heap = next
        guard let key else {
            fatalError("pairing heap exhausted early")
        }
        a[i] = key
    }
}


func benchmark_sort(_ array: inout [Int]) {

    pairing_heap_sort(&array)

}

func is_non_decreasing(_ a: [Int]) -> Bool {
    guard a.count >= 2 else { return true }
    for i in 1..<a.count {
        if a[i - 1] > a[i] { return false }
    }
    return true
}

func same_multiset(_ a: [Int], _ b: [Int]) -> Bool {
    if a.count != b.count {
        return false
    }

    var left = a
    var right = b
    left.sort()
    right.sort()
    return left == right
}

func check_correctness_case(_ label: String, _ input: [Int]) {
    var input = input
    let original = input

    benchmark_sort(&input)

    if !is_non_decreasing(input) {
        fatalError("correctness case \(label): output is not sorted")
    }

    if !same_multiset(input, original) {
        fatalError("correctness case \(label): elements were lost or added")
    }
}

// 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.
func check_correctness_case_within_limit(_ label: String, _ input: [Int]) {
    if input.count > (1 << MAX_POWER) {
        return
    }
    check_correctness_case(label, input)
}

func few_unique_values(_ size: Int, _ unique: Int, _ seed: UInt64) -> [Int] {
    var state = seed
    var result = [Int]()
    result.reserveCapacity(size)
    for _ in 0..<size {
        state ^= state << 13
        state ^= state >> 7
        state ^= state << 17
        result.append(Int(state % UInt64(unique)) + 1)
    }
    return result
}

func run_correctness_checks() {
    check_correctness_case("empty", [])
    check_correctness_case("single", [42])
    check_correctness_case("duplicates", [3, 1, 3, 2, 1, 2])
    check_correctness_case("sorted", [1, 2, 3, 4, 5])
    check_correctness_case("reverse", [5, 4, 3, 2, 1])
    check_correctness_case("all_equal", [7, 7, 7, 7])
    check_correctness_case("skewed_range", [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",
        [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(
            "few_keys_len32_seed_\(seed)",
            few_unique_values(32, 4, UInt64(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", [Int](repeating: 7, count: 256))
    for seed in 1...4 {
        check_correctness_case(
            "few_keys_len256_seed_\(seed)",
            few_unique_values(256, 4, UInt64(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", [Int](repeating: 7, count: 600))
    for seed in 1...4 {
        check_correctness_case_within_limit(
            "few_keys_len2048_seed_\(seed)",
            few_unique_values(2048, 4, UInt64(seed))
        )
    }
}


func shuffled(_ size: Int, seed: UInt64) -> [Int] {
    guard size > 0 else { return [] }

    var v = Array(1...size)
    var state = seed

    if size > 1 {
        for i in stride(from: size - 1, through: 1, by: -1) {
            state ^= state << 13
            state ^= state >> 7
            state ^= state << 17

            let j = Int(state % UInt64(i + 1))
            v.swapAt(i, j)
        }
    }

    return v
}

func micros(_ d: Duration) -> UInt64 {
    let c = d.components
    let fromSeconds = UInt64(c.seconds) * 1_000_000
    let fromAttos = UInt64(max(0, c.attoseconds / 1_000_000_000_000))
    return fromSeconds + fromAttos
}

func padLeft(_ value: String, _ width: Int) -> String {
    if value.count >= width {
        return value
    }
    return String(repeating: " ", count: width - value.count) + value
}

func formatSeconds(_ micros: UInt64) -> String {
    let whole = micros / 1_000_000
    let frac = micros % 1_000_000
    let fracStr = padLeft(String(frac), 6).replacingOccurrences(of: " ", with: "0")
    return "\(whole).\(fracStr)"
}

func input_array(_ size: Int, seed: UInt64) -> [Int] {
    shuffled(size, seed: 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.
func run_once(size: Int, seed: Int) -> (UInt64, Int) {
    var array = input_array(size, seed: UInt64(seed))

    let baseBytes = alloc_track_live()
    alloc_track_reset_peak()

    let start = ContinuousClock.now

    benchmark_sort(&array)

    let elapsed = ContinuousClock.now - start
    let peakBytes = alloc_track_peak()
    let auxBytes = max(0, peakBytes - baseBytes)

    let expected: [Int] = size > 0 ? Array(1...size) : []
    if array != expected {
        fatalError("sort failed with seed \(seed) for size \(size)")
    }

    return (micros(elapsed), auxBytes)
}

func run_child(_ args: [String]) {
    let size = Int(args[2])!
    let seed = Int(args[3])!
    let (elapsedUs, mem) = run_once(size: size, seed: seed)
    print("\(elapsedUs) \(mem)")
}

let args = CommandLine.arguments
if args.count > 1 && args[1] == "--run-once" {
    run_child(args)
} else {
    run_correctness_checks()

    let tableHeader =
        "| \(padLeft("Size", 10)) | " +
        "\(padLeft("Average time (s)", 16)) | " +
        "\(padLeft("Maximum time (s)", 16)) | " +
        "\(padLeft("Average memory (KiB)", 20)) | " +
        "\(padLeft("Maximum memory (KiB)", 20)) |"
    print(tableHeader)
    print("|-----------:|-----------------:|-----------------:|---------------------:|---------------------:|")

    for power in MIN_POWER...MAX_POWER {
        let size = 1 << power

        var totalTime: UInt64 = 0
        var maxTime: UInt64 = 0

        var totalMem = 0
        var maxMem = 0

        for seed in 1...RUNS {
            let process = Process()
            process.executableURL = URL(fileURLWithPath: args[0])
            process.arguments = ["--run-once", "\(size)", "\(seed)"]
            let stdout = Pipe()
            let stderr = Pipe()
            process.standardOutput = stdout
            process.standardError = stderr

            do {
                try process.run()
            } catch {
                fatalError("failed to run benchmark child process: \(error)")
            }
            process.waitUntilExit()

            if process.terminationStatus != 0 {
                let err = String(data: stderr.fileHandleForReading.readDataToEndOfFile(), encoding: .utf8) ?? ""
                fatalError("benchmark child process failed: \(err)")
            }

            let data = stdout.fileHandleForReading.readDataToEndOfFile()
            let stdoutText = String(data: data, encoding: .utf8) ?? ""
            let fields = stdoutText.split(whereSeparator: \.isWhitespace)
            guard fields.count >= 2,
                  let elapsedUs = UInt64(fields[0]),
                  let auxMem = Int(fields[1]) else {
                fatalError("invalid child process output: \(stdoutText)")
            }

            totalTime += elapsedUs
            if elapsedUs > maxTime {
                maxTime = elapsedUs
            }

            totalMem += auxMem
            if auxMem > maxMem {
                maxMem = auxMem
            }
        }

        let avgTime = totalTime / UInt64(RUNS)
        // Memory is summed in bytes and converted to KiB once, after averaging.
        let avgMemKb = totalMem / RUNS / 1024
        let maxMemKb = maxMem / 1024

        let tableRow =
            "| \(padLeft(String(size), 10)) | " +
            "\(padLeft(formatSeconds(avgTime), 16)) | " +
            "\(padLeft(formatSeconds(maxTime), 16)) | " +
            "\(padLeft(String(avgMemKb), 20)) | " +
            "\(padLeft(String(maxMemKb), 20)) |"
        print(tableRow)
    }
}
SWIFT

RUN clang -O2 -fPIC -shared alloc_track.c -o liballoc_track.so -ldl

RUN swiftc -Ounchecked -whole-module-optimization \
    main.swift \
    -o swift-benchmark \
    -L. -lalloc_track \
    -Xlinker -rpath -Xlinker /app

ENV LD_PRELOAD=/app/liballoc_track.so
CMD ["./swift-benchmark"]
"""#
    try dockerfile.write(
        to: workdir.appendingPathComponent("Dockerfile"),
        atomically: true,
        encoding: .utf8
    )

    // Keeping build and run as separate child processes preserves Docker's
    // normal output and the original image tag used by the benchmark skill.
    try runCommand("docker", ["build", "-t", "swift-benchmark", workdir.path])
    try runCommand("docker", ["run", "--rm", "--init", "swift-benchmark"])
} catch {
    fputs("\(error)\n", stderr)
    exit(1)
}