只读排序的时空权衡

这期的算法问题是，输入只读可随机访问、输出只能顺序写的排序问题，允许使用额外空间。时空权衡 (Time-Space Trade-Offs) 问的是，在给定时间（或空间）复杂度下，能做到最小的空间（或时间）复杂度是多少？

（让文章严谨一点，本文一律采用比较模型 + Word RAM）

参考的论文是：

Upper Bounds for Time-Space Trade-offs in Sorting and Selection 这论文讲了三个算法，我们只讲其中的只读排序。
Optimal Time-Space Trade-Offs for Sorting 这篇就是本文的主线。

我们先初步给两个特解。

首先是遍历求最小值，然后反复遍历仅次于上一轮的最小值。空间只需要保存几个指针，即 $O(\log n)$ bit，时间是 $O(n^2)$ ，因此时间空间乘起来是 $O(n^2\log n)$ bit。

第二个方案是用 n 个指针指向输入元素，然后用任意一种 $O(n\log n)$ 的排序完成。这样空间是 $O(n\log n)$ bit，时间是 $O(n\log n)$ ，时间空间乘起来是 $O(n^2\log^2 n)$ bit。

本文的主线会是一个优雅的算法（在我看来），时间乘以空间是 $O(n^2)$ bit，已被证明是最优的。

1. 前置知识

你需要了解基础排序算法、二叉树、堆、位运算、复杂度的基础知识。

同时：

1.1. 位向量

位向量是每个位置存储 0 或 1 的数组，这里需要支持一次操作读写连续的多个比特。这里就只展示范围读的函数：

cpp

struct BitVector {
    int64_t size_ = 0;
    std::vector<uint64_t> data_;

    uint64_t load(int64_t pos, int64_t len) const {
        common::assert_or_throw(pos >= 0 && len >= 0 && pos + len <= size_);
        common::assert_or_throw(len <= n_word_bits);
        if (len == 0) {
            return 0;
        }
        int64_t off = pos % n_word_bits;
        int64_t n_low = n_word_bits - off;
        uint64_t low = data_[pos / n_word_bits] >> off;
        if (len <= n_low) {
            return low & low_mask(len);
        }
        return (low & low_mask(n_low)) | (data_[(pos + len) / n_word_bits] & low_mask(len - n_low))
                                             << n_low;
    }
};

2. 多路归并排序

第一篇论文讲的算法比较简单，分块后块内遍历、块间用堆做多路归并。（顺便讲一下，这和文章主线关系不大）

假设 $s$ 是额外空间 bit 数除以 $\log n$ 。

首先按 $\frac{n}{s}$ 大小分块，一共 $s$ 块。然后每块都遍历一遍求最小值，塞进堆里。

执行 $n$ 次这样的操作：

把堆里的最小值 min 拿出来，写到输出数组。
遍历 min 所在块，求大于 min 的最小值，塞进堆里。

我们来求一下复杂度。

堆单次操作的复杂度是 $O(\log s)$ ， $n$ 次就是 $O(n\log s)$ 。

遍历一块的复杂度是 $O(\frac{n}{s})$ ， $n$ 次遍历就是 $O(\frac{n^2}{s})$ 。

加起来 $O(n\log s + \frac{n^2}{s})$ 。考虑到 $s$ 超过 $\frac{n}{\log n}$ 就达到比较排序的下界 $O(n\log n)$ ，因此只考虑 $s \le \frac{n}{\log n}$ 。也就是 $n\log s<n\log n\le\frac{n^2}{s}$ ，最终化简为 $O(\frac{n^2}{s})$ 。

额外空间复杂度 $O(s\log n)$ bit。

最终时间乘以空间就是 $O(n^2\log n)$ bit。

贴个代码：

cpp

std::vector<RandomIt> heap;
for (int64_t i = 0; i < n_blocks; i++) {
    int64_t block_start = i * block_size;
    int64_t block_end = std::min(size, block_start + block_size);
    auto block_min = next_element(
        first + block_start, first + block_end, std::optional<RandomIt>{}, iter_proj);
    common::assert_or_throw(block_min.has_value());
    heap.push_back(*block_min);
}

std::ranges::make_heap(heap, std::greater{}, iter_proj);
for (int64_t i = 0; i < size; i++) {
    common::assert_or_throw(!heap.empty());
    *output = *heap.front();
    ++output;
    std::ranges::pop_heap(heap, std::greater{}, iter_proj);
    RandomIt min_it = heap.back();
    heap.pop_back();
    int64_t block_id = (min_it - first) / block_size;
    int64_t block_start = block_id * block_size;
    int64_t block_end = std::min(size, block_start + block_size);
    auto block_next = next_element(
        first + block_start, first + block_end, std::optional<RandomIt>{min_it}, iter_proj);
    if (block_next.has_value()) {
        heap.push_back(*block_next);
        std::ranges::push_heap(heap, std::greater{}, iter_proj);
    }
}

3. M-Tree 排序

第二篇论文讲了 M-Tree 数据结构，这是一个类似锦标赛树的结构，可以快速求最小值、单点修改。

假设额外空间是 $O(s)$ bit。

3.1. state 索引

将数组按 $\frac{n}{s}$ 大小分块，一共 $s$ 块。建立一棵满二叉树 M-Tree，叶子节点对应一块。除了叶子每个节点存储一个比特标记 state，表示子树内的最小值位于左子树还是右子树。

查询：从根开始按 state 跳转到叶子，可以 $O(\log s)$ 的复杂度找到最小值所在块，然后 $O(\frac{n}{s})$ 的复杂度遍历这个块找最小值。这里只考虑 $s \le \Theta(\frac{n}{\log n})$ （为什么后面讲），所以最终复杂度 $O(\frac{n}{s})$ 。

更新：去掉最小值后，如何更新 state？只要更新一条树链，从叶子到根依次：求左右子节点各自的最小值（用上面查询的方法），比较两个最小值结果写入 state。需要更新 $O(\log s)$ 个节点，每个节点 $O(\frac{n}{s})$ 的查询最小值的开销，所以复杂度 $O(\frac{n\log s}{s})$ 。

完整算法 $n$ 次查询和更新，复杂度 $O(\frac{n^2\log s}{s})$ ，时间和空间乘起来就是 $O(n^2\log s)$ 。

3.2. μ 索引优化

这个复杂度还是大了，于是论文引入了 $\mu$ 索引的优化。

$\mu_i$ 即节点 $i$ 存储 $d-d_i$ bit（ $d$ 是根到叶子距离， $d_i$ 是根到节点 $i$ 的距离），用于索引最小值所在范围。越靠近根节点，范围越精确。

具体来说，节点先通过 state 标记确定子树局部最小值所在块（简称最小值块），然后将最小值块划分为 $\frac{s}{2^{d_i}}$ 个子块，那么最小值就位于第 $\mu_i$ 个（假设从 0 开始计数）。

$\mu$ 只在块内索引，所以和 state 是独立的功能。

另外有一点就是，每个节点的 $\mu$ 的信息是冗余的（除了二进制最低位，剩下位就是它的某个子节点的 $\mu$ ）。因此从信息量来说，每个节点存 1 bit 也是够的，只是复杂度会变大。本来 $O(1)$ 时间读一个 $\mu$ 的操作会变成 $O(d-d_i)$ 了。

这里有个不属于论文的小优化，我把 state 也改成相同的格式，节点 $i$ 存储 $d-d_i$ bit，这样 state 也能 $O(1)$ 直接查询了。

查询：只需要读根节点的 state 和 $\mu$ ，然后遍历子块找最小值，复杂度 $O(\frac{n}{s^2})$ 。

更新：从叶子到根的每个节点，对于它的左右子节点都进行查询，然后修改 state 和 $\mu$ 。复杂度 $O(\sum_{i=1}^{\log s}1+\frac{n}{2^is})=O(\log s+\frac{n}{s})$ 。同样，只考虑 $s \le \Theta(\frac{n}{\log n})$ ，所以最终复杂度只保留 $O(\frac{n}{s})$ 。

因此完整算法 $n$ 次查询和更新，复杂度 $O(\frac{n^2}{s})$ 。

state 或者 $\mu$ 的空间复杂度，即 $d+2(d-1)+2^2(d-2)+...+2^{d-1}=O(2^d)=O(s)$ （单位 bit）。总空间就是 $O(s)$ bit。

时间和空间乘起来就是 $O(n^2)$ bit。

考虑最小空间，要能存储数组的索引，因此 $s\ge \Theta(\log n)$ 。

最大空间的话， $s=\Theta(\frac{n}{\log n})$ 会达到比较排序的下界 $O(n\log n)$ 。在 $s>\Theta(\frac{n}{\log n})$ 的时候，更新操作的 $O(\log s) \ge O(\log n-\log \log n) = O(\log n)$ ，最后还是能得出不低于 $O(n\log n)$ 的结论的。因此 $s$ 太大了没意义。

3.3. 复杂度下界

论文将唯一元素问题 (unique elements problem) 规约到排序问题。唯一元素问题就是给定一个只读的数组，输出所有值唯一的元素。很显然，将所有元素排序后，按顺序比较就可以了。

因此排序问题和唯一元素问题一样难或更难。

而唯一元素问题的时空下界就是 $\Omega(n^2)$ bit（不知道怎么证，想了解的看论文 A General Sequential Time-Space Tradeoff for Finding Unique Elements）。排序问题的上界是就是这个下界，因此这个算法达到了最优。

3.4. 实现 M-Tree 排序

首先把 state 和 $\mu$ 用位向量紧凑保存，对于第 level 层第 offset 个用一个 data_offset 函数计算偏移（等差乘等比数列求和，直接问 AI 要一个公式）：

cpp

int64_t data_offset(int64_t level, int64_t offset) const {
    return ((height_ - level + 2) * (int64_t{1} << level)) - height_ - 2 +
            (offset * (height_ - level));
}

用 state 和 $\mu$ 定位要查找的子块，封装一个函数：

cpp

std::tuple<int64_t, int64_t> range(int64_t level, int64_t offset) const {
    common::assert_or_throw(level >= 0 && level <= height_);
    common::assert_or_throw(offset >= 0 && offset < (int64_t{1} << level));
    int64_t n_sub_blocks = int64_t{1} << (height_ - level);
    int64_t off = data_offset(level, offset);
    int64_t len = height_ - level;
    int64_t state_v = state_.load(off, len);
    int64_t mu_v = mu_.load(off, len);
    int64_t l = (((offset * n_sub_blocks) + state_v) * block_size_) +
                (mu_v * block_size_ / n_sub_blocks);
    int64_t r = l + (block_size_ / n_sub_blocks);
    return {l, r};
}

min 函数求某个节点下的最小值，用 range 定位子块后，遍历一遍即可，代码略。

update 函数用来更新某个节点，先用 min 函数找左子节点、右子节点的最小值，比较最小值并更新 state 和 $\mu$ ：

cpp

void update(int64_t level, int64_t offset) {
    common::assert_or_throw(level >= 0 && level < height_);
    common::assert_or_throw(offset >= 0 && offset < (int64_t{1} << level));
    auto lchild_min = min(level + 1, offset * 2);
    auto rchild_min = min(level + 1, (offset * 2) + 1);
    int64_t n_level_nodes = int64_t{1} << level;
    int64_t off = data_offset(level, offset);
    int64_t len = height_ - level;
    if (!lchild_min.has_value() && !rchild_min.has_value()) {
        state_.store(off, len, 0);
        mu_.store(off, len, 0);
    } else if (!rchild_min.has_value() ||
                (lchild_min.has_value() && iter_proj_(*lchild_min) < iter_proj_(*rchild_min))) {
        int64_t child_off = data_offset(level + 1, offset * 2);
        int64_t child_len = len - 1;
        state_.store(off, len, state_.load(child_off, child_len));
        mu_.store(off, len,
            (mu_.load(child_off, child_len) << 1) |
                ((lchild_min.value() - first_) / root_sub_block_size_ % n_blocks_ /
                    n_level_nodes % 2));
    } else {
        int64_t child_off = data_offset(level + 1, (offset * 2) + 1);
        int64_t child_len = len - 1;
        state_.store(off, len,
            state_.load(child_off, child_len) | (int64_t{1} << (height_ - level - 1)));
        mu_.store(off, len,
            (mu_.load(child_off, child_len) << 1) |
                ((rchild_min.value() - first_) / root_sub_block_size_ % n_blocks_ /
                    n_level_nodes % 2));
    }
}

初始化就是自底向上调一遍 update，代码略。

最后 pop 就是 min 求全局最小值，然后从最小值叶子到根 update 一遍，代码略。

最后实现排序，只要建树后不断 pop 就行了：

cpp

for (int64_t i = 0; i < size; i++) {
    std::optional<RandomIt> min_opt = tree.min();
    tree.pop();
    *output = *min_opt.value();
    ++output;
}

4. 完整代码

多路归并排序完整实现和测试。

M-Tree 排序完整实现和测试。

5. 结尾

我对论文有一个微小优化，就是上文提到的 state 读写优化。论文将 M-Tree 外面套了第一个算法（多路归并），也能得出相同的结论，只是复杂一点。

但不管怎么说， $O(n^2)$ bit 的结论非常优雅简洁。

时空权衡排序在理论计算机领域是很著名的问题，至少 AI (Gemini) 是这么说的。但是 AI 没说对复杂度是多少（好像被单位是 bit 还是字 (Word) 绕进去了），说明也没有那么著名。

只读排序的时空权衡

1. 前置知识 ​

1.1. 位向量 ​

2. 多路归并排序 ​

3. M-Tree 排序 ​

3.1. state 索引 ​

3.2. μ 索引优化 ​

3.3. 复杂度下界 ​

3.4. 实现 M-Tree 排序 ​

4. 完整代码 ​

5. 结尾 ​