浅谈高级筛法

CheeseFunction · 2025-10-12 13:56:43 · 算法·理论

修复：

\text{[update 2025/10/20 17:18] 修改一处文本，但不影响理解}
\text{[update 2025/10/20 18:17] 修改三处标题，但不影响理解}

Parts1: Min_25:

简单来说，Min_25 是一种快速求一个积性函数 f(x) 前缀和的筛法，对于这个 f(x)，我们有以下要求和限制：

首先对于 f(p) 其中 p 为质数，可以被表示为低次多项式或者可以快速求出。
其次对于 f(p^k) 其中 p 为质数，也能够快速计算。

那么到底怎么做呢，我们现在来说。首先，我们考虑将质数部分先处理掉，定义一个辅助函数 g(n,j), 其定义为在埃拉托斯特尼筛法，也就是埃氏筛中未被前 j 个质数筛掉（即最小质因子 > p_j）的数的 f'(x) 之和，这里为了后面的讲解，使用一个完全积性函数 f'(x) 近似原函数 f(x) 在质数位置的值。现在我们容易得出，对于 g(n,0) 有：

g(n,0) = \sum_{i=2}^n f'(i)

考虑如果 f'(i) 是一个简单多项式的化，我们可以直接用公式求和。接下来对于所有 g(n,j) 我们考虑从 g(n,0) 推广得出以下递推式：

g(n, j) = \begin{cases} g(n, j-1) - f'(p_j) \cdot \left[ g\left( \left\lfloor \frac{n}{p_j} \right\rfloor, j-1 \right) - \sum_{i=1}^{j-1} f'(p_i) \right], & \text{If } p_j^2 \le n \\ g(n, j-1), & \text{If } p_j^2 > n \end{cases}

对于 g(n,j) 的上限，我们认为 j 应该尽可能取到一个与 \sqrt{n} 最接近且严格大于其的质数，然后我们就会发现这个极限的 g(n,j) 就是对于所有 p \le n 且 p 为质数情况下，f'(p) 之和，搞定！接下来考虑计算合数部分，我们定义一个 s(n,i) 表示所有最小质因子大于 p_i 的数的 f(i) 之和。那么我们的目标就是计算 s(n,0)，最后加上 f(1) (如果需要的话), 我们考虑 s(n,i) 递归式为

S(n,i) = (g(n,j_{\max}) - \sum_{j = 1}^j f(p_j)) + \sum_{k > i} \sum_{e \ge 1} [f(p_k^e) \cdot (S(\lfloor \frac{n}{p^e_k}\rfloor , k) + [e > 1])]

其中我们计算的是所有大于 p_i 的质数的 f(p) 的和，后面的那个一长串二重求和是枚举合数的最小质因子 p_k，以及其指数。然后 S(\lfloor \frac{n}{p^e_k}\rfloor , k) 表示递归计算除去最小质因子 p_k^e 后剩余部分（其最小质因子大于 p_k）的 f 值。然后后面的 [e > 1] 是一个指示函数，主要是考虑 p_k^e 本身用的。

递归边界为 n < p_i 时满足 s(n,i) = 0，最终就是 \sum_{i = 1}^n f(i) = s(n,0) + f(1)。

现在我们考虑证明 Min_25 时间复杂度，先证明对于所有 S(m,i)，其中 m 只有最多 \sqrt{n} 种取值，原因非常简单，在计算 S(n,i) 的时候，参数 m 总是形如 \lfloor \frac{n}{k} \rfloor，就是数论分块的证明。当然这个是小菜，我们接下来质数分布的估计，需要前置知识切比雪夫定理，也就是我们规定对于第 k 个质数 p_k \sim k \ln k 且 \pi(n) \sim \frac{n}{\ln n}。这个就是定理，非常简单，而且根据上面的递推式，也很好证明 S(m,i) 所有有意义的 i 满足 p_i \le \sqrt{m}。

然后现在我们说， Min_25 时间复杂度为 O(\frac{n^{\frac{3}{4}}}{\log n}) 级别的，我们分析递归计算 S(m,i) 的总代价，我们定义 V = \{\left\lfloor \frac{n}{k} \right\rfloor : 1 \le k \le n\}, 然后根据 p_i 与 m 的大小关系，分为：

Case 1: p_i \le m^{\frac{1}{4}}
Case 2: m^{\frac{1}{4}} < p_i \le \sqrt{n}

现在我们逐步分析：

对于 Case 1 的，我们有：首先，根据递归式，我们知道 m 的取值只有 \sqrt{n} 中（定理 1），然后，对于所有 m，满足 p_i \le m^{\frac{1}{4}} 的 i 的个数为 \pi(m^{\frac{1}{4}}) = O(\frac{m^{\frac{1}{4}}}{\log n}) 大概是这个量级的，这个十个粗略估计的结果，但是也足以判断 Case 1 是一个小部分，我们忽略不计。
接下来对于 Case 2 这个关键的地方，我们做以下事情：首先由递归式可以快速推导出我们在考虑上面的东西时，只需要考虑 e = 1 的情况，因此对于每个 S(m,i) (Case 2)，产生的调用总次数为 \pi(\sqrt{m}) - i 量级。现在将计算总量。我们令 T(m,i) 表示以 S(m,i) 为根的递归子树中的状态数，我们有：

T(m, i) = 1 + \sum_{\substack{k > i \\ p_k \le \sqrt{m}}} T(\lfloor \frac{m}{p_k} \rfloor, k)

但是我们关心的时 T(n,0)，可以尝试使用积分近似，有：

\sum_{\substack{m \in V \\ m \ge n^{\frac{1}{2}}}} [\pi(\sqrt{m}) - \pi(m^{\frac{1}{4}})] \sim \sum_{\substack{m \in V \\ m \ge n^{\frac{1}{2}}}} \left[ \frac{2\sqrt{m}}{\log m} - \frac{4m^{\frac{1}{4}}}{\log m} \right]

考虑到 m 较为密集，用积分近似求和：

\int_{n^{\frac{1}{2}}}^n \left( \frac{2\sqrt{x}}{\log x} - \frac{4x^{\frac{1}{4}}}{\log x} \right) \frac{dx}{x}

这里 \frac{dx}{x} 是因为 m 的取值密度约为 \frac{dm}{m}，然后考虑计算：

\int_{n^{\frac{1}{2}}}^n \frac{2\sqrt{x}}{x\log x} dx = \int_{n^{\frac{1}{2}}}^n \frac{2}{x^{\frac{1}{2}}\log x} dx

然后令 t = \log x 且 dx = e^tdt 则：

= \int_{\frac{1}{2}\log n}^{\log n} \frac{2e^{\frac{t}{2}}}{t} dt \sim \frac{4n^{\frac{1}{2}}}{\log n}

然后后面这一项太小了，直接不看，因此我们估计时间复杂度在 O(\frac{n^{\frac{3}{4}}}{\log n}), 当然这个不是准确值，也有别的说法，但是大概在这个量级是没错的。

Part 2: Meissel-Lehmer

严格意义上来说，Meissel-Lehmer 没有 Min_25 那么实用，但是还是很棒的，我们来详细说一下，Meissel-Lehmer 是求：

\pi(x) = \#\{p \le x | p \text{ is prime}\}

你可能说，哎呀这个筛法不是轻轻松松吗，实则要是我 n 大一点点直接原地升天，这个时候，我们就要升级算法了，我们来说一下算法流程，首先，定义：

\phi(x, a) = \#\{n \le x : p \mid n \implies p > p_a\}

其中 p_a 表示第 a 个素数，其中 p_1 = 2。换句话说，\phi(x,a) 是 1 \dots x 中所有不被前 a 个素数整除的正整数的个数，即：用埃氏筛筛掉前 a 个素数后剩下的数的个数，和上面 Min_25 是差不多的（但是就是要换个说法，主要是写作时间线和参考的文献不太一样，就有些问题多多谅解）显然，这些数包括：

所有大于 p_a 的素数（\le x）。
以及 1。

所以容易推导得到：

\pi(x) = \phi(x, a) + a - 1, \quad \text{If } p_a \le \sqrt{x} < p_{a+1}.

接下来我们就要来推导 \phi 的公式了，你只需要记住有：

\phi(x, a) = \phi(x, a-1) - \phi\left( \frac{x}{p_a}, a-1 \right)

从 a-1 到 a 时，我们筛掉那些不被前 a-1 个素数整除但能被 p_a 整除的数。\ 这些数可以写成 p_a \cdot m，其中 m 不被前 a-1 个素数整除，且 p_a m \le x，即 m \le \lfloor \frac{x}{p_a} \rfloor。\ 满足 m 不被前 a-1 个素数整除的个数正是 \phi\left( \frac{x}{p_a}, a-1 \right)。

但是先需要定义 \phi(x,0) = 0, 也就是没筛任何数的时候，所有正整数留下。

但是直接递归计算 \phi 还是无法令人接受，但是我们有一个非常关键的观察：

定义 P_k(x,a) 为：

P_k(x, a) = \#\{ n \le x : n = q_1 q_2 \dots q_k,\ q_1 > q_2 > \dots > q_k > p_a \}

取 a = \pi(x^{\frac{1}{3}})，b = \pi(\sqrt{x})。

因为每个不被前 a 个素数整除的数 \le 1 且所有素数大于 p_a。

此时我们容易写出：

\phi(x, a) = 1 + \sum_{k=1}^\infty P_k(x, a)

因为每个不被前 k 个素数整除的 >1 的数，其所有素因子 > p_a，且素因子个数 k \geq 1。

显然我们可以打个表：

P_1(x, a) = \pi(x) - a
P_2(x, a) = \#\{ pq \le x : p > q > p_a \}
P_3(x, a) = \#\{ pqr \le x : p > q > r > p_a \}

然后瞎猜有：

\phi(x, a) = 1 + [\pi(x) - a] + P_2(x, a) + P_3(x, a) + \dots

接下来我们分开处理，取 a = \pi(x^{\frac{1}{3}})，则：

p_{a+1} > x^{\frac{1}{3}}

然后对于 k \geq 3 最小 P_3 对应 p_{a+1} p_{a+2} p_{a+3} > (x^{\frac{1}{3}})^3 = x 因此：

P_k(x, a) = 0 \quad \text{If } k \ge 3

所以：

\phi(x, a) = 1 + \pi(x) - a + P_2(x, a)

整理以后：

\pi(x) = \phi(x, a) + a - 1 - P_2(x, a) \tag{1}

然后考虑固定第 i 个素数，其中 a < i \le b = \pi(\sqrt{x})，有 p_j > p_i 使得 p_i p_j \le x 也就是 p_j \le x/p_i, 容易写出这样的 p_j 一共有：

\pi\left( \frac{x}{p_i} \right) - \pi(p_i)

个，但是 \pi(p_i) = i，由此：

P_2(x, a) = \sum_{i=a+1}^{b} \left[ \pi\left( \frac{x}{p_i} \right) - i \right] \tag{2}

最后，将 (1) 和 (2) 写出来，就变成了：

\pi(x) = \phi(x, \pi(\frac{1}{3})) + \pi(\frac{1}{3}) - 1 - \sum_{i=\pi(x^{\frac{1}{3}})+1}^{\pi(\sqrt{x})} \left[ \pi\left( \frac{x}{p_i} \right) - i \right]

计算步骤很简单，预处理素数到 \sqrt{x} 然后用 DP 或别的递归方法计算 \phi(x,a)，然后计算 P_2 带入求和得到 \pi(x)。

现在是时间复杂度证明时间！

定义计算 \phi(x,a) 时出现的不同第一个参数的集合为：

S = \left\{ \left\lfloor \frac{x}{n} \right\rfloor : n = p_{i_1} p_{i_2} \cdots p_{i_k},\ i_1 > i_2 > \cdots > i_k,\ n \le x,\ k \ge 0 \right\}.

每个 n 是前 a 个素数中若干个的乘积（无平方因子）。

有引理如 |S| = O\left( \frac{x^{\frac{2}{3}}}{\log^2 x} \right)，证明很简单，如下：

将 S 按 n \le x^{\frac{1}{3}} 与 n > x^{\frac{1}{3}} 分成两部分。

若 n \le x^{\frac{1}{3}}，则 \lfloor \frac{x}{n} \rfloor \ge x^{\frac{2}{3}}，且不同的 n 的数量不超过 x^{\frac{1}{3}}，因此对应的 m 的数量 \le x^{\frac{1}{3}}。
若 n > x^{\frac{1}{3}}，则 m = \lfloor \frac{x}{n} \rfloor < x^{\frac{2}{3}}。\ 对于固定的 m，n = \lfloor \frac{x}{m} \rfloor 唯一确定，因此这样的 m 的数量 \le x^{\frac{2}{3}}。\ 但需注意 n 必须是无平方因子数且素因子来自前 a 个素数。\ 记 m = \lfloor \frac{x}{n} \rfloor，则 n \approx \frac{x}{m}。\ 当 m \le x^{\frac{2}{3}} 时，n \ge x^{\frac{1}{3}}。\ 对于 m \le x^{\frac{1}{2}}，n \ge x^{\frac{1}{2}}，这样的 n 的数量由大素数积控制，可用组合数上界：\ 不超过 \binom{a}{r} 对某个 r，且 n \ge x^{\frac{1}{3}} 时 r 较大，总数经细致平衡可得 O\left( \frac{x^{\frac{2}{3}}}{\log^2 x} \right)。

综合两部分，|S| = O\left( \frac{x^{\frac{2}{3}}}{\log^2 x} \right)。

然后又有引理有：计算 \phi(m, a) 对所有 m \in S 的时间为 O\left( \frac{x^{\frac{2}{3}}}{\log^2 x} \cdot a \right)。

证明依旧简单，但是因为篇幅问题，不证明了。接下来证明 P_2 的时间复杂度：

P_2 = \sum_{i=a+1}^{b} \left[ \pi\left( \frac{x}{p_i} \right) - i \right].

其中项数 b-a = \pi(\sqrt{x}) - \pi(x^{1/3}) = \Theta\left( \frac{\sqrt{x}}{\log x} \right)。

需要计算每个 \pi(\frac{x}{p_i})，其中 p_i \ge x^{\frac{1}{3}}, 然后考虑预计算 \pi(m) 对于所有 m \le x^{\frac{2}{3}} 可以使用埃氏筛完成，因此 P_2 的求和耗时为：

O\left( \frac{\sqrt{x}}{\log x} \right).

接下来汇总：

注意到瓶颈过高，要优化（Lagarias–Miller–Odlyzko 1985）：通过只保留 m = \lfloor \frac{x}{t} \rfloor 对于 t \le x^{\frac{1}{3}}，则计算 \phi 的时间为：

O(a \cdot |S|) = O\left( \frac{x^{\frac{1}{3}}}{\log x} \cdot x^{\frac{1}{3}} \right) = O\left( \frac{x^{\frac{2}{3}}}{\log x} \right)

预计算一些小 m 的 \phi 可以在 O(x^{\frac{2}{3}} \log \log x) 内完成。

最后汇总得到：

O\left( \frac{x^{\frac{2}{3}}}{\log x} + x^{\frac{2}{3}} \log \log x + \frac{\sqrt{x}}{\log x} \right) = O\left( \frac{x^{\frac{2}{3}}}{\log x} \right).

如果你还是优化大佬的话，兴许可以来到：

O\left( \frac{x^{\frac{2}{3}}}{\log^2 x} \right)

还是比 Min_25 快多了。

Part3: Powerful Number 筛

Powerful Number 筛，我们也称作 PN 筛，是最近几年才出现的一种非常厉害的筛法。

定义：正整数 n 为 Powerful Number 当且仅当 n 的所有质因子 p 都有 p^2 \mid n。

性质有下：

PN 的数量，我们记 \text{PN}(N) 为不超过 n 的 PN 的数量，有 \text{PN}(N) \sim \frac{\zeta(\frac{3}{2})}{\zeta(3)} \sqrt{N} \approx 2.173 \sqrt{N} 如果你想写好一点，当然也可以说 \text{PN}(N) \approx \frac{\sqrt{N}}{\log N} 量级。

那么这玩意对我筛 f(n) 前缀和有什么用啊，其核心思路在于我们需要先找一个辅助的函数 g(n) 满足：

看到这个辅助函数，也许你就会想到卷积，确实，定义 h 满足：

f = g \times h \quad \Leftrightarrow \quad f(n) = \sum_{d|n} g(d) h\left(\frac{n}{d}\right)

也就是狄利克雷卷积，其中 h 必须为积性函数。

然后考虑从 h 入手，先计算 h(p^k)。

在 n = p^e 时，f(p^e) = \sum_{i=0}^e g(p^i) h(p^{e-i})，这是卷积的关系。
特别的。如果 e = 1（素数），则 f(p) = g(p)h(1) + g(1)h(p), 带入条件有：g(p) = g(p) + h(p) \quad \Rightarrow \quad h(p) = 0

一条关键的结论：h(p) 为 0 对于所有素数 p 成立。质数部分就算完了。

然后考虑 h 非 0 值和 PN 的关系，由于 h 为积性函数且 h(p) = 0 在 p 为质数下。那么对于包含一次质因子的 n 都有 h(n) = 0，这个非常好证明。

然后我们就要捣鼓前缀和了，前方高能！

由 f = g \times h：

f(n) = \sum_{d|n} g(d) h\left(\frac{n}{d}\right)

带入求前缀和：

S_f(n) = \sum_{i=1}^n f(i) = \sum_{i=1}^n \sum_{d|i} g(d) h\left(\frac{i}{d}\right)

令 k = \frac{i}{d} 则 i = dk 且 d \mid i，交换求和：

S_f(n) = \sum_{d=1}^n g(d) \sum_{k=1}^{\lfloor \frac{n}{d} \rfloor} h(k)

由于 h(k) 仅在 k 是 PN 时非 0，则我们可以尝试枚举所有可能的 k，从原始的二层求和展开：

S_f(n) = \sum_{i=1}^n \sum_{d|i} g(d) h\left(\frac{i}{d}\right)

固定 k 以后：

S_f(n) = \sum_{k=1}^n h(k) \sum_{d=1}^{\lfloor \frac{n}{k} \rfloor} g(d) = \sum_{k=1}^n h(k) \cdot G\left(\left\lfloor \frac{n}{k} \right\rfloor\right)

根据我们上面的性质有：

S_f(n) = \sum_{\substack{k=1 \\ k \text{ is PN}}}^n h(k) \cdot G\left(\left\lfloor \frac{n}{k} \right\rfloor\right)

现在开始，PN 已经被我们推导完了核心部分！

计算办法很简单，由卷积关系在 n = p^e 处：

f(p^e) = \sum_{i=0}^e g(p^i) h(p^{e-i})

分离 i = 0 项：

f(p^e) = g(1)h(p^e) + \sum_{i=1}^e g(p^i) h(p^{e-i})

由于 g(1) = 1，得：

h(p^e) = f(p^e) - \sum_{i=1}^e g(p^i) h(p^{e-i})

所以算法的步骤就是：

找出辅助函数，满足上面要求。
筛出所有 \le \sqrt{n} 的质数，预处理一下 G(m) 在 m = \lfloor \frac{n}{k} \rfloor。
最后计算每一个 h，最后对 h 进行求和就完了。

又来到了时间复杂度分析的部分，但是由于这个过于简单，就随便写一下，你只需要知道 PN 个数是 \sqrt{n} 的大致就知道时间复杂度了。

References:

[1].《Min_25 筛学习笔记》， Luogu 用户 @hzoi_liuchang

[2]. Oi-Wiki (Min_25)(Powerful Number)

[3].《一些特殊的数论函数求和问题》（非原文）