O(n) 原地归并的不稳定版本

讲一个很有意思的学术算法，一个线性时间原地归并算法。

这一期是不稳定版本，下一期我们会讲稳定版本。以稳定版本为基础，归并排序可以同时满足稳定、原地、 $O(n\log n)$ 时间，可以说非常神奇了。

算法主要参考这两篇论文，细节上会有修改：

Practical in-place merging (opens new window) 这是不稳定版本。
Fast Stable Merging and Sorting in Constant Extra Space (opens new window) 这是稳定版本。

WikiSort (opens new window) 也值得一看，是稳定版本。

另外 An In-Place Sorting with O(n log n) Comparisons and O(n) Moves (opens new window) 这个排序比本文的结论更强，多了 $O(n)$ 次移动的限制，可能是理论上最完美的排序。但是过于复杂看不懂喵。

# 1. 原地归并是什么

原地合并两个有序数组，输入是一个长度为 n 的数组 arr 和一个位置 pos，arr[0] ... arr[pos - 1] 属于第一个有序数组 A，arr[pos] ... arr[n - 1] 属于有序数组 B，merge 的结果是合并两个有序数组为 arr[0] ... arr[n]。

“原地”其实有两种含义，一是算法的结果直接写回原数组，比如 C++ std::inplace_merge 就是这个含义；二是 $O(1)$ 的额外空间复杂度，有时会允许 $O(\log n)$ 的递归栈，这篇文章的“原地”都是这个含义。

# 2. 一些前置算法

双指针归并和区间旋转在我之前文章有讲，我就直接复制过来。

# 2.1. 双指针归并算法

两个指针指向两个数组的开始位置，不断把较小元素放到 buffer 数组，对应指针向后移动。最后把 buffer 里的元素移动回原数组。

这个过程如下图所示：

merge_to

这个算法在 C++ std::inplace_merge (opens new window) 里亦有记载（如果空间够的话）。

# 2.2. 缓冲区双指针归并算法

传统的双指针归并不是原地算法，需要进行改造。假设一个数组，按顺序是长度 r（或大于 r）的缓冲区，长度 l 的数组 A，长度 r 的数组 B，使用双指针归并将 A 和 B 归并在数组开头。

缓冲区的元素不能丢弃，因此使用交换来移动元素。

void buffered_merge(int* buf, int* left, int* mid, int* right) {
    int* a_ptr = left;
    int* b_ptr = mid;
    while (a_ptr < mid || b_ptr < right) {
        if (b_ptr == right || (a_ptr < mid && *a_ptr <= *b_ptr)) {
            std::swap(*buf, *a_ptr);
            buf++;
            a_ptr++;
        } else {
            std::swap(*buf, *b_ptr);
            buf++;
            b_ptr++;
        }
    }
}

这里表示归并结果的 buf 指针是不会追上表示 A 数组的 a_ptr 指针。因为如果 a_ptr 指针右移，buf 指针和 a_ptr 的距离不变；而 b_ptr 指针右移的次数只有 r 次（数组 B 的长度），刚好是缓冲区大小。

# 2.3. rotate 区间旋转算法

原地算法基本绕不开 rotate（区间旋转），把两个相邻区间 [A B] 原地变成 [B A]，保持区间内部顺序不变。这个方法不唯一，最经典的做法是三次翻转法（或手摇算法）：先分别翻转区间 A 和区间 B，再整体翻转整个区间 [A B]，就能得到 [B A]。

这个过程如下图所示：

rotate

这个算法在 C++ std::rotate (opens new window) 里亦有记载。当然标准库会根据数据量采用不同算法，这里就不深究了。

# 2.4. 旋转归并算法

旋转归并算法也是合并两个有序数组，基于 rotate，因此不需要额外空间。

一开始是两个有序数组 [A B]，不断进行下面步骤：

取 A 的第一个数 A[0]。
将 B 划分为小于 A[0] 的区间 B1、大于等于 A[0] 区间 B2。
通过 rotate 算法将 A 和 B1 交换。
A[0] 已经确定位置，剩下 A 的元素作为新的 A，B2 作为新的 B，回到步骤 1。

这个过程其实是用“批量”插入排序的方式完成 merge。

设 $l$ 是 A 的长度， $r$ 是 B 的长度。因为算法会循环 $l$ 次，所以 A 的每个数移动 $O(l)$ 次，B 的每个数移动 $O(1)$ 次，总时间复杂度 $O(l^2+r)$ 。

同理只要反一下就能做到 $O(l+r^2)$ 复杂度。

void rotate_merge(int* left, int* mid, int* right) {
    while (left < mid && mid < right) {
        int* split_b = mid;
        while (split_b < right && *split_b < *left) {
            split_b++;
        }
        std::rotate(left, mid, split_b);
        left += (split_b - mid) + 1;
        mid = split_b;
    }
}

void rotate_merge_reverse(int* left, int* mid, int* right) {
    std::rotate(left, mid, right);
    rotate_merge(left, left + (right - mid), right);
}

# 3. 不稳定原地归并 in O(n) time

这是一个分块算法，我们定义 $s=\lfloor \sqrt{n} \rfloor$ 。

# 3.1. 分块和对齐

把 A, B 数组划分为大小 s 的块，不整除的部分（非对齐部分）用 rotate 算法移动到 B 数组的后面。不整除的部分在最后一步会处理。

void inplace_merge(int* left, int* mid, int* right) {
    size_t len1 = mid - left;
    size_t len2 = right - mid;
    size_t s = std::sqrt(len1 + len2);
    size_t len1_aligned = len1 / s * s;
    size_t len2_aligned = len2 / s * s;
    size_t aligned = len1_aligned + len2_aligned;

    std::rotate(left + len1_aligned, mid, mid + len2_aligned);

    /* 后续算法 */
}

# 3.2. 块间排序

所有块按首元素升序排序，块内的顺序不变。

由于交换两个块的代价很大，我们需要每个块只要交换 $O(1)$ 次的排序，用什么排序呢？没错，就是选择排序。

比较次数 $(n/s)^2 \approx (\sqrt{n})^2=n$ ，交换次数 $s(n/s)=n$ ，因此总复杂度 $O(n)$ 。

void blocked_select_sort(int* left, int* right, size_t s) {
    for (int* i = left; i < right; i += s) {
        int* min_block = i;
        for (int* j = i + s; j < right; j += s) {
            if (*min_block > *j) {
                min_block = j;
            }
        }
        if (min_block != i) {
            std::swap_ranges(i, i + s, min_block);
        }
    }
}

# 3.3. 块间合并

这是最核心的一个步骤。

我们需要大小为 s 的缓冲区来完成这一步骤，考虑把第一块的数划为缓冲区。

利用第 i 块为缓冲区对第 i + 1 块、第 i + 2 块进行双指针归并。注意不能丢失缓冲区的数，所以要用缓冲区双指针归并算法。这一步后缓冲区的数会转移到第 i + 2 块的位置。
交换第 i + 1 块和第 i + 2 块，如果是 i + 2 是最后一块就不用。
i 增加 1，回到步骤 1。

块间合并完后除了缓冲区，所有块都整体有序。

为什么算法可行呢？因为，我们每次归并的是首元素最小的两个块（不考虑缓冲区）。

如果首元素最小的两块一开始属于同一数组，假设是数组 A。显然数组 B 的最小值大于等于第二块首元素，也就大于等于前一块的所有数。因此第一块就是最小的 s 个数。

如果首元素最小的两块一开始属于不同的数组（一块属于 A，一块属于 B），这两块已经包含了数组 A 和数组 B 的最小的 s 个数，因此最小的 s 个数只会出现在这两块里。

单次双指针归并 $O(s)$ 次交换，一共 $O(n/s)$ 次归并，所以总复杂度 $O(n)$ 。

void blocked_merge(int* left, int* right, size_t s) {
    size_t n = right - left;
    for (size_t i = 0; i + s * 2 < n; i += s) {
        buffered_merge(left + i, left + i + s, left + i + s * 2, left + i + s * 3);
        if (i + s * 3 < n) {
            std::swap_ranges(left + i + s, left + i + s * 2, left + i + s * 2);
        }
    }
}

# 3.4. 处理尾部元素

我们把缓冲区移动到后面，和非对齐部分一起统称尾部元素。非对齐部分长度不大于 2s，缓冲区长度 s，因此尾部元素长度不超过 3s。

首先对尾部元素进行排序，用任意一个平方复杂度的排序即可，例如冒泡。这里的复杂度是 $O((3s)^2)=O(n)$ 。

用旋转归并算法把前面部分（长度不超过 n）和尾部元素（长度不超过 3s）合并为一个有序数组。这里的复杂度是 $O(l+r^2)=O(n+(3s)^2)=O(n)$ 。

void inplace_merge(int* left, int* mid, int* right) {
    size_t len1 = mid - left;
    size_t len2 = right - mid;
    size_t s = std::floor(std::sqrt(len1 + len2));
    size_t len1_aligned = len1 / s * s;
    size_t len2_aligned = len2 / s * s;
    size_t aligned = len1_aligned + len2_aligned;

    /* 前面几个步骤 */

    std::sort(left + aligned - s, right);
    rotate_merge_reverse(left, left + aligned - s, right);
}

# 4. 补个代码

#include <algorithm>
#include <cmath>
#include <print>

/* 粘贴上文的 buffered_merge */
/* 粘贴上文的 rotate_merge 和 rotate_merge_reverse */
/* 粘贴上文的 blocked_select_sort */
/* 粘贴上文的 blocked_merge */

void inplace_merge(int* left, int* mid, int* right) {
    size_t len1 = mid - left;
    size_t len2 = right - mid;
    size_t s = std::floor(std::sqrt(len1 + len2));
    size_t len1_aligned = len1 / s * s;
    size_t len2_aligned = len2 / s * s;
    size_t aligned = len1_aligned + len2_aligned;

    // 分块和对齐
    std::rotate(left + len1_aligned, mid, mid + len2_aligned);
    // 块间排序
    blocked_select_sort(left, left + aligned, s);
    // 块间合并
    blocked_merge(left, left + aligned, s);
    // 处理尾部元素
    std::sort(left + aligned - s, right);  // 这里应该用冒泡排序，偷懒了
    rotate_merge_reverse(left, left + aligned - s, right);
}

int main() {
    int arr[10] = {2, 5, 7, 9, 1, 3, 4, 6, 8, 10};
    inplace_merge(arr, arr + 4, arr + 10);
    for (int i : arr) {
        std::print("{} ", i);
    }
    std::println("");
}

# 5. 稳定原地归并预告

不妨分析一下上面算法里什么操作破坏了稳定性。第一，“块间排序”里的选择排序是不稳定排序。第二，“块间合并”会把缓冲区打乱，同样导致不稳定。

我们会利用缓冲区作为 label 把算法稳定下来，思路类似带 index 的快速排序。细节有很多，就放到下一期讲了。

← 一个 PRAM 并行算法模拟器