只读算法系列第二篇。这次的神秘算法是选择算法（求第 k 小值），但是输入数组是只读的、可以随机访问。你只能申请 $O(1)$ 大小的工作空间。

在这样的条件下，可以做到复杂度 $O(n^{1+\epsilon})$，$\epsilon$ 是任意正数，这是目前最优结论。（看起来接近线性复杂度，其实是大于 $O(n\log n)$ 的，优化空间还有很多）

算法本身很简单，就是复杂度分析很难。完整复杂度分析可以看看原论文。

参考论文是 Selection from Read-only Memory and Sorting with Minimum Data Movement。

1. 算法框架

选择算法最基本的框架就是，确定两个过滤器，一个上界 L 一个下界 R。每一轮在位于 L 和 R 之间的元素（候选）里挑选一个 pivot，计算 pivot 的排名。根据排名和 k 的关系，更新 L / R 为 pivot 的后继 / 前驱。

因此挑选策略是算法的重点。

template <typename RandomIt, typename IterProj>
RandomIt select(RandomIt first, RandomIt last, RandomIt lower_it, RandomIt upper_it, int64_t k,
    int64_t n_layers, IterProj iter_proj) {
    while (true) {
        RandomIt pivot_it = /*...*/;  // 挑一个 [lower, upper] 范围里的元素
        int64_t rank = count_in_range(first, last, lower_it, pivot_it, iter_proj) - 1;
        if (rank == k) {
            return pivot_it;
        }
        if (rank < k) {
            lower_it = min_greater(first, last, pivot_it, iter_proj);
            k -= rank + 1;
        } else {
            upper_it = max_smaller(first, last, pivot_it, iter_proj);
        }
    }
}

2. 朴素算法 (A1)

朴素算法的挑选策略就很简单，挑第一个 $[L, R]$ 范围里的元素。

for (RandomIt it = first; it < last; it++) {
    if (in_range(it, lower_it, upper_it, iter_proj)) {
        pivot_it = it;
        break;
    }
}

这样最多进行 $n$ 轮，每轮都要 $O(n)$（遍历一遍），最终复杂度 $O(n^2)$。

3. 根号算法 (A2)

假设 $[L,R]$ 范围里有 r 个元素（候选），把数组分成大约 $\sqrt r$ 块，块大小是 $\frac{n}{\sqrt r}$。

块内有的元素是候选，有的不是。可以证明存在一块包含至少 $\sqrt r$ 个候选，因为假设所有块都小于 $\sqrt r$ 个候选，乘以块数就是小于 $r$ 个候选，矛盾。

在包含至少 $\sqrt r$ 个候选的块里，截取恰好 $\sqrt r$ 个候选的前缀区间，然后用朴素算法 (A1) 求中位数。这个中位数就是我们要的 pivot。

因为 pivot 一定大于 $\frac{\sqrt r}{2}$ 个元素、小于 $\frac{\sqrt r}{2}$ 个元素，每轮可以淘汰 $\frac{\sqrt r}{2}$ 个元素，即 $[L,R]$ 区间大小可以减少 $\frac{\sqrt r}{2}$。

---

复杂度分析环节。

选取合适的块需要 $O(n)$。朴素算法的复杂度是区间大小乘以候选数，即 $O(\frac{n}{\sqrt r}\cdot \sqrt r)=O(n)$。

设算法时间 $T(r)$，每轮淘汰 $\frac{\sqrt r}{2}$ 个元素，列公式 $T(r) = T(r-\frac{\sqrt r}{2})+\Theta(n)$。

• 注意到 $T(r-\frac{\sqrt r}{2})<T(r-\frac{\sqrt r}{2}+\frac{1}{16})=T((\sqrt r-\frac{1}{4})^2)$。
• 于是 $T(r)=T((\sqrt{r})^2)<T((\sqrt r-\frac{i}{4})^2)+\Theta(in)$。
• 令 $i=4\sqrt r$，得到 $T(r) = O(n\sqrt r)$，代入初始候选 $r = n$ 得到 $O(n\sqrt n)$。

所以复杂度就是 $O(n\sqrt n)$。

4. 通用递归算法 (As)

进一步，我们可以把数组分为大约 $r^{1/s}$ 块，块大小是 $\frac{n}{r^{1/s}}$。

找大于等于 $r^{1-1/s}$ 个候选的块，并截取 $r^{1-1/s}$ 个候选的区间。在区间里递归 $A_{s-1}$ 算法求中位数，作为 pivot。

---

复杂度是 $O(n^{1+1/s})$，至于为什么，这公式看得头疼，不推导了。

---

int64_t block_size =
    std::ceil(static_cast<double>(size) /
                std::pow(n_candidates, 1.0 / static_cast<double>(n_layers)));
int64_t n_blocks = (size + block_size - 1) / block_size;
int64_t recursive_max =
    std::ceil(std::pow(n_candidates, 1.0 - (1.0 / static_cast<double>(n_layers))));

RandomIt recursive_start = last;
RandomIt recursive_end = last;
int64_t recursive_count = 0;
for (int64_t i = 0; i < n_blocks; i++) {
    RandomIt block_start = first + (i * block_size);
    RandomIt block_end = first + std::min(size, (i + 1) * block_size);
    int64_t block_count =
        count_in_range(block_start, block_end, lower_it, upper_it, iter_proj);
    if (block_count > recursive_count) {
        recursive_start = block_start;
        recursive_end = block_end;
        recursive_count = block_count;
    }
}
while (recursive_count > recursive_max) {
    recursive_end--;
    if (in_range(recursive_end, lower_it, upper_it, iter_proj)) {
        recursive_count--;
    }
}
pivot_it = select_recursive(recursive_start, recursive_end, lower_it, upper_it,
    (recursive_count - 1) / 2, n_layers - 1, iter_proj);

5. 推论

论文的给出的最终复杂度其实是 $O(2^s s!n^{1+1/s})$。

假设 s 是常数，那么就是本文的主要结论 $O(n^{1+\epsilon})$，因为常数 s 可以任意大。

另一个推论是取 $s=O(\sqrt{\frac{\log n}{\log\log n}})$，可以达到 $O(n^{1+o(1)})$ 复杂度。

6. 最好的算法在哪里

这是个开放问题。

简单说明两个小知识：如果证明了解决某个问题不能小于 $O(f(n))$ 复杂度，那么 $O(f(n))$ 就是这个问题的复杂度下界；如果真的给出了 $O(f(n))$ 复杂度的算法，那么 $O(f(n))$ 就是这个问题的复杂度上界。

我们已经给出了 $O(n^{1+\epsilon})$ 的算法，这是已知最好的上界。

下界呢？很遗憾已知下界是平凡的，即输入复杂度 $O(n)$。

没人能回答 $O(n)$ 或者 $O(n\log n)$ 能不能解决这个问题。这也能说明只读算法的水是真的深。

7. 完整代码

完整实现和测试。

8. 结尾

这次算法的 $O(1)$ 空间其实是简化了问题，还有一些研究时间复杂度 - 空间复杂度曲线的论文，我打算写一篇《只读算法系列》介绍这个曲线是怎么回事。

另外，多趟选择算法我已经看完了，但是这实现好复杂啊。我研究算法一直是古法编程写着玩的，这下不得不启动 vibe coding 了。