位向量的只读选择算法

这一期来个相对轻松的算法（实现起来没花几天）。输入只读、可随机访问，我们要实现 $O(n)$ 时间复杂度、 $O(n)$ bit 额外空间的选择算法（求第 k 小）。因为上一期《位向量与 rank select》实现了位向量，剩下的部分就没那么难了。

参考论文：Selection from read-only memory with limited workspace

1. 前置知识

1.1. 位向量与 rank select

上一期《位向量与 rank select》讲的数据结构，用额外空间 $O(n)$ bit 实现位向量的常数时间 rank 和 select，其中：

位向量：每个位置存储 0 或 1 的数组。
rank：给定一个整数 $i$ ，求位向量前 $i$ 位有多少个 1。
select：给定一个整数 $k$ ，求位向量第 $k$ 个 1 的位置（假设“第 X 个”是从 0 开始计数）。

1.2. BFPRT 算法

BFPRT 又叫 Median of Medians，是一个比较有名的选择算法。在之前文章《O(n) 原地选择的不稳定版本》我们讨论过，这里就简要写一下。

第一步，中位数取中阶段，每 5 个数一组，取组内的中位数，然后把这些中位数用递归 BFPRT 的方法再取中位数作为 pivot。

第二步，淘汰阶段，根据 pivot 将数组进行三路划分，划分为小于 pivot、等于 pivot、大于 pivot 的三个区间。如果 k 落在了第 1 或 3 个区间就递归 BFPRT。

可以保证 $O(n)$ 的最坏时间复杂度，怎么证明可以看往期《O(n) 原地选择的不稳定版本》。

在代码里我们用 select_ref 表示这个算法：

cpp

template <typename RandomIt, typename Proj = std::identity>
void select_ref(RandomIt first, RandomIt nth, RandomIt last, Proj proj = {});

2. 只读选择算法

算法的框架依然是：选择一个近似中位数，比较它和 k（目标是求第 k 小）的关系并淘汰一边的元素。

这个算法独特的地方是，淘汰的元素用位向量标记，用位向量 select 可以快速访问候选（未淘汰）元素，非常好用。

但是支持 rank select 的位向量不支持修改（或者说新增淘汰标记），于是就有了——

2.1. Wavelet Stack

Wavelet Stack 和本文的算法非常契合。

一开始 Wavelet Stack 的所有元素都是候选。淘汰时，在栈上压入一个位向量，这个位向量的长度是淘汰前候选个数，对本轮淘汰的元素标记 0（淘汰标记）。也就是下面一层的位向量 1 的数量，是上一层位向量的长度。

由于每轮会淘汰大致一半的元素，可以预期 Wavelet Stack 的层数是 $O(\log n)$ ，空间是 $O(n)$ bit。现在缺少信息，真正的证明在后面进行。

唯一需要的查询操作是 select，从栈顶到栈底 select 一遍就可以了。时间复杂度和栈的大小呈正比。

cpp

struct WaveletStack {
    int64_t size_ = 0;
    int64_t n_word_bits_ = 0;
    std::vector<BitVectorType> stack_;

    template <typename Placement>
    void update(Placement placement) {
        auto bit_vector = BitVectorType::create(
            count(),
            [this, placement](int64_t k) {
                int64_t index = select(k);
                return placement(index);
            },
            n_word_bits_);
        stack_.push_back(bit_vector);
    }

    int64_t select(int64_t k) const {
        if (stack_.empty()) {
            assert_or_throw(k >= 0 && k < size_);
        }
        for (auto it = stack_.rbegin(); it != stack_.rend(); ++it) {
            k = it->select(k);
        }
        return k;
    }
};

2.2. 近似中位数

求 Median of Medians。

在 $O(n)$ bit 的限制下，我们最多存储 $O(\frac{n}{\log n})$ 个指针（因为一个地址是 $\log n$ bit）。

将候选者分为大小为 $\frac{n}{\log n}$ 的块，每块把候选者的指针保存在 buffer 里。用 BFPRT 取出 buffer 的中位数，保存到 medians 数组里。

块之间 buffer 是可以复用的，因此空间复杂度 $O(n)$ bit。medians 最多只有 $\log n$ 个指针，只要 $\log^2 n$ bit。

最后求 medians 数组的中位数作为近似中位数。

代码里的 BitVector 参数，上游会传 WaveletStack 进来。只不过对 median_of_medians 来说它是 WaveletStack 还是 BitVector 不重要。

cpp

template <typename BitVector, typename RandomIt, typename IterProj>
RandomIt median_of_medians(const BitVector& bit_vector, RandomIt first, IterProj iter_proj) {
    int64_t size = bit_vector.size();
    int64_t n_candidates = bit_vector.count();
    assert_or_throw(size > 0 && n_candidates > 0);
    int64_t block_size = std::max(size / std::max(ceil_log2(size), int64_t{1}), int64_t{1});
    int64_t n_blocks = (n_candidates + block_size - 1) / block_size;
    std::vector<RandomIt> medians;
    for (int64_t i = 0; i < n_blocks; i++) {
        int64_t block_start = i * block_size;
        int64_t block_end = std::min(block_start + block_size, n_candidates);
        std::vector<RandomIt> buffer;
        for (int64_t j = block_start; j < block_end; j++) {
            buffer.push_back(first + bit_vector.select(j));
        }
        int64_t mid_index = (static_cast<int64_t>(buffer.size()) - 1) / 2;
        select_ref(buffer.begin(), buffer.begin() + mid_index, buffer.end(), iter_proj);
        medians.push_back(buffer[mid_index]);
    }
    int64_t mid_index = (static_cast<int64_t>(medians.size()) - 1) / 2;
    select_ref(medians.begin(), medians.begin() + mid_index, medians.end(), iter_proj);
    return medians[mid_index];
}

2.3. 算法主体

最后来实现前文讲到的算法框架。

cpp

while (true) {
    RandomIt pivot_it = median_of_medians(wavelet_stack, first, iter_proj);
    int64_t pivot_rank = 0;
    for (int64_t i = 0; i < wavelet_stack.count(); i++) {
        if (iter_proj(first + wavelet_stack.select(i)) < iter_proj(pivot_it)) {
            pivot_rank++;
        }
    }
    if (pivot_rank == k) {
        return pivot_it;
    }
    if (pivot_rank < k) {
        wavelet_stack.update(
            [&](int64_t i) { return iter_proj(first + i) > iter_proj(pivot_it); });
        k -= pivot_rank + 1;
    } else {
        wavelet_stack.update(
            [&](int64_t i) { return iter_proj(first + i) < iter_proj(pivot_it); });
    }
}

2.4. 重复元素的处理

因为输入数组是可随机访问的，因此可以看到我们工作空间不存元素而是指针。除了省空间外，还有一个好处就是，比较元素发现相同，可以继续比较指针。这样的比较方式直接可以视为所有元素互不相同了。

代码只要一行：

cpp

auto iter_proj = [proj](RandomIt it) { return std::pair{proj(*it), it}; };

2.5. 复杂度

淘汰近似中位数左边时，因为近似中位数是块中位数的中位数，因此较小的块中位数（占比一半）所在块的一半被淘汰。右边也同理。因此近似中位数落在 $0.25n$ 到 $0.75n$ 的范围里。

易得 Wavelet Stack 需要的空间不超过 $0.75n+0.75^2n+0.75^3n+...=O(n)$ （单位 bit），并且 select 总时间是 $O(n+2\cdot0.75n+3\cdot 0.75^2n+4\cdot 0.75^3n)=O(n)$ 。

3. 完整代码

完整实现和测试。

4. 扩展算法

在工作空间不足 $O(n)$ bit 时，论文对这个算法进行了调整。用到了 Frederickson 算法，是我之后打算讲的算法，这里就当作黑盒了。

假设工作空间为 $O(s)$ bit，用多趟只读选择算法（Frederickson 算法）把候选数降到 s（假设这些是中间候选），然后就能执行位向量的只读选择算法了。但问题是，候选元素分散在各个位置，需要一个结构可以快速索引。

答案依旧是分桶（其实就是分块），把原数组分成 $\frac{n}{s}$ 大小的桶，最多 $s$ 个桶。构造 C 为计数向量（位向量），每个桶记录若干个 1（数量为桶内的中间候选数），桶之间用 1 个 0 分隔。这样计数向量最多只要 $2s-1$ bit。

构造 H 为 Wavelet Stack，用 $O(s)$ bit 维护 s 个中间候选元素的淘汰情况。

然后就是遍历操作，在位向量的只读选择算法里，select 都只用于遍历。

假设上一次执行了 select(i - 1)，这次的 select(i)：先用 H.select(i) 拿到中间候选的索引 j，然后用 C.select(j) - j 拿到桶号。如果和上一次是同一个桶，就在原数组里从 select(i - 1) 开始顺序扫描，直到找到候选；否则在新桶的起始位置开始顺序扫描。

最终复杂度是 $O(n\log^*\frac{n}{s}+\frac{n\log n}{\log s})$ 。这里要说明一下，左边看起来始终比右边小（比如 $s$ 是常数时），实则不然：令 $s=\sqrt n$ ，左边是 $O(n\log^*\sqrt n)$ ，右边直接 $O(n)$ ，右边略胜一筹。所以公式没法简化了。

这个算法我就简单介绍一下，不打算细讲了。

5. 最好的结果在哪里

那么， $O(n)$ 时间复杂度下，最小额外空间是多少？这是一个开放问题。

位向量的只读选择算法给出了 $O(n)$ bit 的算法，这是已知最好的上界。

对于下界，Comparison-Based Time–Space Lower Bounds for Selection 这篇论文证明了不能小于 $\Omega(n^{\epsilon})$ bit， $\epsilon$ 可以任意小。好像叫亚多项式。比如 $\Omega(\log n)$ 这种对数复杂度是小于亚多项式的。

往期《只读 O(1) 空间选择算法》也是开放问题，我拿 deepseek-v4.1f 尝试突破已知的上下界，可想而知失败了。至少说明这类问题不简单。

6. 结尾

只读选择问题的最小空间、最小时间都讲了，介于两者的中间那块还没讲。目前还不知道中间的最好算法是什么，之后再看看吧。

位向量的只读选择算法

1. 前置知识 ​

1.1. 位向量与 rank select ​

1.2. BFPRT 算法 ​

2. 只读选择算法 ​

2.1. Wavelet Stack ​

2.2. 近似中位数 ​

2.3. 算法主体 ​

2.4. 重复元素的处理 ​

2.5. 复杂度 ​

3. 完整代码 ​

4. 扩展算法 ​

5. 最好的结果在哪里 ​

6. 结尾 ​