ポリフェーズマージソートを使用する

ポリフェーズマージソート (poly-phase merge sort) は、複数のテープ(またはファイル)に分散した整列済みランを、フィボナッチ分布に沿って段階的にマージしていく。

テープ本数が少ない環境でも、各パスでほぼすべてのテープを稼働させ、バランスマージよりパス数を抑えられる場合がある。

本記事のデモとベンチマークでは、3 本の仮想テープを主記憶上のベクタで模擬する。実際の外部整列では置換選択などで初期ランを作り、ラン数がフィボナッチ数でないときはダミーランで分布を調整する。

  1. 初期ラン生成: 配列を固定長(例: 32 要素)の区間に区切り、各区間を整列してランとする。
  2. フィボナッチ分布: 3 本テープのうち 1 本を空け、残り 2 本へラン数比を連続するフィボナッチ数(例: {2, 3}, {3, 5})に近づけるよう分配する。
  3. ポリフェーズマージ: 2 本のソーステープから先頭ランを 1 組ずつ取り出し、空きテープへマージする。
  4. テープローテーション: 出力テープの役割を循環させ、再び 2 本から 1 本へのマージを繰り返す。
  5. 完了: 全要素が 1 本のテープ上の 1 ランにまとまったら配列へ書き戻す。
procedure create_runs(A, run_size)
  split A into chunks of run_size, sort each chunk into a run

procedure distribute_fibonacci(runs, tapes[3])
  target = smallest Fibonacci number >= length(runs)
  put runs on tape 1 and tape 2 in Fibonacci ratio; tape 0 empty

procedure polyphase_pass(tapes[3])
  while tape 1 and tape 2 both have runs
    merged = merge(tape1.pop_front(), tape2.pop_front())
    tape0.push_back(merged)
  rotate tape roles cyclically

procedure polyphase_merge_sort(A)
  runs = create_runs(A)
  tapes = distribute_fibonacci(runs)
  while total run count on all tapes > 1
    polyphase_pass(tapes)
  copy final run back into A

テープ本数が少ない外部整列向けで、フィボナッチ分布によりマージパス数を抑えられる。ラン生成とマージを合わせた時間はおおむね \(O(n \log n)\) で、安定なマージを使えば全体も安定ソートになる。作業領域はランとテープバッファに依存する。

テープドライブが高価だった時代のポリフェーズマージは、限られた I/O チャネルを稼働させ続ける典型例として学ぶ価値がある。

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

マージソートは 2 列のマージを繰り返す。ポリフェーズはテープが少ない外部整列向きに、フィボナッチ分布でマージ先を回転させる。

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

Size Average time (s) Maximum time (s) Average memory (KiB) Maximum memory (KiB)
256 0.000009 0.000107 4 4
512 0.000020 0.000116 12 12
1024 0.000040 0.000173 17 17
2048 0.000096 0.000295 51 51
4096 0.000269 0.000666 70 70
8192 0.000515 0.000984 204 204
16384 0.000927 0.003010 280 280
32768 0.002005 0.003154 816 816
65536 0.005373 0.018082 1120 1120
131072 0.010083 0.017612 2240 2240
262144 0.027018 0.042414 4480 4480
計測に使用したコードを表示する

#!/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
func merge_values(_ left: UnsafeBufferPointer<Int>, _ right: UnsafeBufferPointer<Int>) -> [Int] {
    var out = [Int]()
    out.reserveCapacity(left.count + right.count)
    var l = 0
    var r = 0
    while l < left.count && r < right.count {
        if left[l] <= right[r] {
            out.append(left[l])
            l += 1
        } else {
            out.append(right[r])
            r += 1
        }
    }
    while l < left.count {
        out.append(left[l])
        l += 1
    }
    while r < right.count {
        out.append(right[r])
        r += 1
    }
    return out
}

func merge_values(_ left: [Int], _ right: [Int]) -> [Int] {
    left.withUnsafeBufferPointer { leftBuf in
        right.withUnsafeBufferPointer { rightBuf in
            merge_values(leftBuf, rightBuf)
        }
    }
}



fileprivate let NUM_TAPES = 3
fileprivate let RUN_SIZE = 32

fileprivate func merge_runs(_ left: [Int], _ right: [Int]) -> [Int] {
    merge_values(left, right)
}

fileprivate func create_runs(_ a: UnsafeMutableBufferPointer<Int>, _ run_size: Int) -> [[Int]] {
    var runs: [[Int]] = []
    var i = 0
    while i < a.count {
        let end = min(i + run_size, a.count)
        var run = Array(a[i..<end])
        run.sort()
        runs.append(run)
        i = end
    }
    return runs
}

fileprivate func next_fibonacci_at_least(_ n: Int) -> (Int, Int) {
    var prev = 1
    var curr = 1
    while curr < n {
        let next = prev + curr
        prev = curr
        curr = next
    }
    return (prev, curr)
}

fileprivate func distribute_fibonacci(_ runs: [[Int]]) -> [[[Int]]] {
    var tapes: [[[Int]]] = [[], [], []]
    let n = runs.count
    if n == 0 {
        return tapes
    }
    if n == 1 {
        tapes[1].append(runs[0])
        return tapes
    }

    let (fib_prev, fib_target) = next_fibonacci_at_least(n)
    let dummies = fib_target - n
    let on_tape2 = max(0, fib_prev - dummies)
    let on_tape1 = n - on_tape2

    for (idx, run) in runs.enumerated() {
        if idx < on_tape1 {
            tapes[1].append(run)
        } else {
            tapes[2].append(run)
        }
    }
    return tapes
}

fileprivate func count_runs(_ tapes: [[[Int]]]) -> Int {
    tapes.reduce(0) { $0 + $1.count }
}

fileprivate func rotate_tapes(_ tapes: inout [[[Int]]]) {
    tapes.swapAt(0, 1)
    tapes.swapAt(1, 2)
}

fileprivate func polyphase_pass(_ tapes: inout [[[Int]]]) -> Bool {
    var merged = false
    while !tapes[1].isEmpty && !tapes[2].isEmpty {
        let left = tapes[1].removeFirst()
        let right = tapes[2].removeFirst()
        tapes[0].append(merge_runs(left, right))
        merged = true
    }
    return merged
}

fileprivate func merge_all_remaining(_ tapes: inout [[[Int]]]) -> [Int] {
    var all: [[Int]] = []
    for tape in tapes {
        all.append(contentsOf: tape)
    }
    while all.count > 1 {
        let a = all.removeFirst()
        let b = all.removeFirst()
        all.append(merge_runs(a, b))
    }
    return all.popLast() ?? []
}

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

func polyphase_merge_sort(_ a: UnsafeMutableBufferPointer<Int>) {
    if a.count <= 1 {
        return
    }
    let runs = create_runs(a, RUN_SIZE)
    if runs.count <= 1 {
        if let r = runs.first {
            for i in 0..<a.count {
                a[i] = r[i]
            }
        }
        return
    }

    var tapes = distribute_fibonacci(runs)
    var idle = 0
    while count_runs(tapes) > 1 {
        if polyphase_pass(&tapes) {
            rotate_tapes(&tapes)
            idle = 0
        } else {
            idle += 1
            if idle > NUM_TAPES * 4 {
                break
            }
            rotate_tapes(&tapes)
        }
    }

    let result: [Int]
    if count_runs(tapes) == 1 {
        result = tapes.compactMap { $0.first }.first ?? []
    } else {
        result = merge_all_remaining(&tapes)
    }
    for i in 0..<a.count {
        a[i] = result[i]
    }
}


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

    polyphase_merge_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)
}