スムースソートを使用する

スムースソート (smooth sort) は、レオナルド木を積み重ねた森を維持し、末尾を最大値として取り出していく。

ヒープソート (heap sort) と同様に最悪時間計算量は \(O(n \log n)\) で、追加の補助記憶域は \(O(1)\) に抑えられるインプレースなアルゴリズムだが、入力がすでに昇順に近いほど作業量が少なくなる適応的 (adaptive) 性質を持つ。

ヒープ上の交換で要素が動くため不安定である。

通常の二分ヒープでは、根が配列の一端にあり、最大値を取り出すたびに根と末尾を交換してヒープを底側から縮めていく。この配置は「ヒープとしての整理」と「昇順に書き出す向き」が噛み合わず、すでにソート済みの入力でも一度ヒープ状に混ぜ直すため、最良でも必ず \(O(n \log n)\) 程度の仕事が残る。

スムースソートでは、二分木ではなくレオナルド数に基づくレオナルド木を積み重ねた森を維持する。各木はヒープ条件を満たし、さらに各木の根の値は左から右へ弱い順(非減少)として並べる。そうすると常に右端が全体の最大になり、末尾から確定した最大値を「取り出す」処理がすでに整っている入力ではほとんど余計な比較を要しない。

レオナルド数 L(k) は

  • L(0) = L(1) = 1
  • L(k) = L(k-2) + L(k-1) + 1(\(k \ge 2\))

で定められ、1, 1, 3, 5, 9, 15, … と続く。サイズ L(k) のレオナルド木は、大きい方の子が左、小さい方が右になるよう二つのより小さいレオナルド木と根で構成される。

任意の長さ n は、高々 \(O(\log n)\) 個の互いに異なるレオナルド数の和として表せる。スムースソートはこの性質を利用して森に含まれる木の本数を常に対数オーダーに抑え、整列済みに近い入力では \(O(n)\) に近づき、最悪計算量でも \(O(n \log n)\) を保つ適応型ソートとなる。

  1. 森の表現: 左から右へ根が非減少になるよう、互いに異なるサイズのレオナルド木を複数本並べた森を維持する。どの位置に木の根があるかはビットマスクと末尾木のサイズ(オフセット)で表す。
  2. 第1段階(構築): インデックス i を 1 から n - 1 まで増やし、森の符号を更新したうえで新根位置 i に対し、単一木内の沈下 sift_in または隣木との整合 interheap_sift でヒープ条件を満たす。
  3. 第2段階(確定): 右端は常に全体の最大なので、i を n - 1 から 2 まで減らしながら末尾木を縮小または 2 子木へ分割し、生じた根に interheap_sift を適用して森の不変条件を保つ。
  4. 終了: 森に残る 2 要素以下は根が左から右へ非減少のため、すでに昇順であり追加の sift は不要。
procedure sift_in(A, rootIdx, size)
  // インデックス size のレオナルド木 1 本の中でヒープ条件を満たすよう根から下げる。
  // L[k]=L[k-2]+L[k-1]+1、子方向は二分木のサイズにより左右どちらかへ潜る。
  if size < 2 then return
  tmp = A[rootIdx]
  r = rootIdx
  sz = size
  loop
    right = r - 1
    left = right - L[sz - 2]
    if A[right] < A[left] then
      candidate = left
      nsz = sz - 1
    else
      candidate = right
      nsz = sz - 2
    if A[candidate] <= tmp then break
    A[r] = A[candidate]
    r = candidate
    sz = nsz
    if sz <= 1 then break
  A[r] = tmp

procedure interheap_sift(A, rootIdx, heap_state)
  // heap_state は「どの桁に木の根があるか」(マスク)と「現在見ている木の順 order」(オフセット)などを束ねたもの。
  tmp = A[rootIdx]
  r = rootIdx
  state = heap_state anchored at rootIdx
  loop while mask(state) <> 1
    maxValue = tmp
    if order(state) > 1 then
      right = r - 1
      left = right - L[order(state) - 2]
      maxValue = max(maxValue, A[left], A[right])
    next = r - L[order(state)]
    if A[next] <= maxValue then break
    A[r] = A[next]
    r = next
    // マスクを調整して左隣の木へカーソル移動
    march_one_tree_left(state)
  A[r] = tmp
  sift_in(A, r, order(state))

procedure smooth_sort(A)
  n = length(A)
  if n <= 1 then return
  heap_state ← initial_singleton_forest_encoding()
  // 第 1 段階(ヒープ化)
  for i = 1 to n - 1
    advance_leonardo_forest(heap_state, i)
    if wide_bottom_condition(i, heap_state, n) then sift_in(A, i, order_at_tail(heap_state)) else interheap_sift(A, i, heap_state)
  // 第 2 段階(右端から順に確定・森の更新)
  for i = n - 1 down to 2
    if order_at_tail(heap_state) < 2 then shrink_trailing_trees(heap_state)
    else
      split_rightmost_into_two_children(heap_state, i)
      for child_root c in newborn_pair_roots()
        interheap_sift(A, c, heap_state)

実装の複雑さと定数倍の大きさから、汎用ライブラリの sort として採用されることは稀である。

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

ヒープソートと同様にインプレースで最悪計算量 \(O(n \log n)\) だが、レオナルド木の森により整列済み入力では \(O(n)\) に近づく適応型ソートである。

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

Size Average time (s) Maximum time (s) Average memory (KiB) Maximum memory (KiB)
256 0.000022 0.000340 0 0
512 0.000043 0.000122 0 0
1024 0.000089 0.000169 0 0
2048 0.000190 0.000335 0 0
4096 0.000446 0.001147 0 0
8192 0.000860 0.001557 0 0
16384 0.001707 0.002865 0 0
32768 0.003663 0.014935 0 0
65536 0.008300 0.043144 0 0
131072 0.017308 0.052062 0 0
262144 0.036789 0.074240 0 0
計測に使用したコードを表示する

#!/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


let LEONARDO: [Int] = [
    1, 1, 3, 5, 9, 15, 25, 41, 67, 109, 177, 287, 465, 753, 1219, 1973, 3193,
    5167, 8361, 13529, 21891, 35421, 57313, 92735, 150049, 242785, 392835,
    635621, 1028457, 1664079, 2692537, 4356617, 7049155, 11405773, 18454929,
    29860703, 48315633, 78176337, 126491971, 204668309, 331160281, 535828591,
    866988873, 1402817465, 2269806339, 3672623805,
]

func smooth_sift_in(_ a: UnsafeMutableBufferPointer<Int>, _ root_idx: Int, _ size: Int) {
    if size < 2 {
        return
    }
    let tmp = a[root_idx]
    var root = root_idx
    var sz = size
    while true {
        let right = root - 1
        let left = right - LEONARDO[sz - 2]
        let next: Int
        let next_size: Int
        if a[right] < a[left] {
            next = left
            next_size = sz - 1
        } else {
            next = right
            next_size = sz - 2
        }
        if a[next] <= tmp {
            break
        }
        a[root] = a[next]
        root = next
        sz = next_size
        if sz <= 1 {
            break
        }
    }
    a[root] = tmp
}

func smooth_interheap_sift(
    _ a: UnsafeMutableBufferPointer<Int>,
    _ root_idx: Int,
    _ mask: Int,
    _ offset: Int
) {
    let tmp = a[root_idx]
    var root = root_idx
    var hmask = mask
    var hoffset = offset
    while hmask != 1 {
        var max = tmp
        if hoffset > 1 {
            let right = root - 1
            let left = right - LEONARDO[hoffset - 2]
            max = Swift.max(Swift.max(max, a[left]), a[right])
        }
        let next = root - LEONARDO[hoffset]
        if a[next] <= max {
            break
        }
        a[root] = a[next]
        root = next
        while true {
            hmask >>= 1
            hoffset += 1
            if hmask & 1 != 0 {
                break
            }
        }
    }
    a[root] = tmp
    smooth_sift_in(a, root, hoffset)
}

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

func smooth_sort(_ a: UnsafeMutableBufferPointer<Int>) {
    let n = a.count
    if n <= 1 {
        return
    }
    var mask = 1
    var offset = 1
    for i in 1..<n {
        if mask & 2 != 0 {
            mask = (mask >> 2) | 1
            offset += 2
        } else if offset == 1 {
            mask = (mask << 1) | 1
            offset = 0
        } else {
            mask = (mask << (offset - 1)) | 1
            offset = 1
        }
        let wide_bottom =
            (mask & 2 != 0 && i + 1 < n)
            || (offset > 0 && 1 + i + LEONARDO[offset - 1] < n)
        if wide_bottom {
            smooth_sift_in(a, i, offset)
        } else {
            smooth_interheap_sift(a, i, mask, offset)
        }
    }
    for i in (2..<n).reversed() {
        if offset < 2 {
            while true {
                mask >>= 1
                offset += 1
                if mask & 1 != 0 {
                    break
                }
            }
        } else {
            let ch1 = i - 1
            let ch0 = ch1 - LEONARDO[offset - 2]
            mask &= ~1
            for ch in [ch0, ch1] {
                mask = (mask << 1) | 1
                offset -= 1
                smooth_interheap_sift(a, ch, mask, offset)
            }
        }
    }
}


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

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