算法学习笔记：Min_25 筛
P2441M
·
·
算法·理论
一些记号
下文中令:
Min_25 筛
Min_25 筛采用素数拟合的思路，我们分别考虑质数和合数处的 f 值之和。先来考虑
\sum_{i=1}^n[i\in\mathbb{P}]f(i)
我们构造完全积性函数 F，使得其在质数处的取值和 f 相等。那上式显然变成了
\sum_{i=1}^n[i\in\mathbb{P}]F(i)
由此考虑 DP。令
$$
g(n,x)=\sum_{i=1}^n[i\in\mathbb{P}\lor\operatorname{lpf}(i)>p_x]F(i)
$$
考虑如何从 $g(n,x-1)$ 转移到 $g(n,x)$，显然此时限制条件变强。不难发现去掉的恰好是那些 $\operatorname{lpf}$ 为 $p_x$ 的数，即
$$
\Delta=g(n,x-1)-g(n,x)=\sum_{i=1}^{n}[\operatorname{lpf}(i)=p_x]F(i)
$$
尝试刻画 $\operatorname{lpf}(i)=p_x$ 的条件。注意到 $F$ 是完全积性函数，我们可以提一个 $F(p_x)$ 出来，这样子剩余部分应该满足 $\operatorname{lpf}\geq p_x$。考察 $g\left(\left\lfloor\frac{n}{p_x}\right\rfloor,x-1\right)$，我们发现这仅仅多算了那些 $<p_x$ 的质数处的 $F$ 值之和。由于我们只考虑 $\leq\lfloor\sqrt{n}\rfloor$ 的质数，所以不难预处理出 $sp_x=\sum_{i=1}^xF(p_i)$。于是我们得到了递推式：
$$
g(n,x)=g(n,x-1)-F(p_x)\left(g\left(\left\lfloor\frac{n}{p_x}\right\rfloor,x-1\right)-sp_{x-1}\right)
$$
根据杜教筛那一套理论，我们只需要 $g\left(\left\lfloor\frac{n}{i}\right\rfloor,x\right)(p_x^2\leq \left\lfloor\frac{n}{i}\right\rfloor)$ 处的值，可以证明对这些点值 DP，时间复杂度是 $\mathcal{O}\left(\frac{n^{\frac{3}{4}}}{\log{n}}\right)$ 的。
实现时，我们需要给 $\left\lfloor\frac{n}{i}\right\rfloor$ 编号。这个挺典的，我们开两个 $\mathcal{O}(\sqrt{n})$ 的索引数组 $id_1$ 和 $id_2$。若 $\left\lfloor\frac{n}{i}\right\rfloor\leq\sqrt{n}$，则将编号存在 $id_1,\left\lfloor\frac{n}{i}\right\rfloor$ 中，否则令 $j=\left\lfloor\frac{n}{\left\lfloor\frac{n}{i}\right\rfloor}\right\rfloor$，将编号存在 $id_{2,j}$ 中即可。
注意我们最终只需要用到 $g(n,\pi(\lfloor\sqrt{n}\rfloor))=\sum_{i=1}^n[i\in\mathbb{P}]F(i)$，因此 $g$ 的第二维可以用滚动数组去掉。
回到原题，我们类似地设计 DP，令
$$
S(n,x)=\sum_{i=1}^n[\operatorname{lpf}(i)>p_x]f(i)
$$
还是分成质数和合数两部分考虑。对于前者，不难得出其贡献为
$$
g(n,\pi(\lfloor\sqrt{n}\rfloor))-sp_x
$$
由于不保证 $f$ 是完全积性函数，我们不能只提单个质因子出来，需要一次性全部提出。也就是我们先枚举 $p_i(i>x)$，再枚举其指数 $c$，这样就可以把 $f(p_i^c)$ 提出来了。剩余部分应该满足 $\operatorname{lpf}>p_i$，也就是 $S\left(\left\lfloor\frac{n}{p_i^c}\right\rfloor,i\right)$，不过由于我们不把 $1$ 纳入 $S,g$ 的计算范围内，所以这样会漏掉 $f(p_i^c\times 1)$，当 $c=1$ 时，$f(p_i)$ 我们已经在质数的部分计算过了，不用加回去，而当 $c>1$ 时，我们则需要把 $f(p_i^c)$ 加回去。由此得到递推式：
$$
S(n,x)=g(n,\pi(\lfloor\sqrt{n}\rfloor))-sp_x+\sum_{i=x+1}^{p_i^2\leq n}\sum_{c=1}^{p_i^c\leq n}f(p_i^c)\left(S\left(\left\lfloor\frac{n}{p_i^c}\right\rfloor,i\right)+[c>1]\right)
$$
这部分的时间复杂度被证明是 $\mathcal{O}(n^{1-\varepsilon})$ 的。
最终答案就是 $S(n,0)+f(1)$。
根据上面的过程，不难发现 Min_25 筛要求 $f$ 和我们构造出的 $F$ 具备以下性质：
- $F(p)=f(p)$ 可以快速计算，$F$ 的前缀和可以快速计算。
- $f(p^c)$ 可以快速计算。
这些性质实际上是不难满足的。
注意我们实际上可以把 $f(p)$ 分解成若干个完全积性函数 $F(p)$ 的线性组合，在计算 $S(n,x)$ 用到 $\sum_{i=x+1}^{\pi(\lfloor\sqrt{n}\rfloor)}f(p_i)$ 时，用同样的方式线性组合对应的 $g$ 和 $sp$ 的值即可。
### [Luogu P5325 【模板】Min_25 筛](https://www.luogu.com.cn/problem/P5325)
板子题。本题中 $f(p)=p^2-p$，构造 $F_1(p)=p$，$F_2(p)=p^2$ 即可。
代码：
```cpp
ll n, sn, num[N << 1];
int sz, pr[N / 10], sp1[N / 10], sp2[N / 10];
int cnt, id1[N], id2[N];
int g1[N << 1], g2[N << 1];
bool v[N];
inline void sieve(int n) {
for (int i = 2; i <= n; ++i) {
if (!v[i]) {
pr[++sz] = i;
sp1[sz] = add(sp1[sz - 1], i);
sp2[sz] = add(sp2[sz - 1], 1ll * i * i % MOD);
}
for (int j = 1; j <= sz && i * pr[j] <= n; ++j) {
v[i * pr[j]] = 1;
if (i % pr[j] == 0) break;
}
}
}
inline int calc(ll x, int y) {
if (pr[y] >= x) return 0;
int k = x <= sn ? id1[x] : id2[n / x];
int ans = sub(sub(g2[k], g1[k]), sub(sp2[y], sp1[y]));
for (int i = y + 1; i <= sz && 1ll * pr[i] * pr[i] <= x; ++i)
for (ll j = pr[i], c = 1; j <= x; j *= pr[i], ++c) {
ll j2 = j % MOD;
cadd(ans, j2 * (j2 - 1) % MOD * (calc(x / j, i) + (c > 1)) % MOD);
}
return ans;
}
int main() {
ios::sync_with_stdio(0), cin.tie(0);
cin >> n, sieve(sn = sqrt(n));
for (ll l = 1, r, x; l <= n; l = r + 1) {
num[++cnt] = x = n / l, r = n / x;
ll x2 = x % MOD;
g1[cnt] = sub(((x2 * (x2 + 1)) >> 1) % MOD, 1);
g2[cnt] = sub(x2 * (x2 + 1) % MOD * (x2 * 2 % MOD + 1) % MOD * inv6 % MOD, 1);
if (x <= sn) id1[x] = cnt;
else id2[n / x] = cnt;
}
for (int i = 1; i <= sz; ++i)
for (int j = 1; j <= cnt && 1ll * pr[i] * pr[i] <= num[j]; ++j) {
int k = (num[j] / pr[i]) <= sn ? id1[num[j] / pr[i]] : id2[n / (num[j] / pr[i])];
csub(g1[j], 1ll * pr[i] * sub(g1[k], sp1[i - 1]) % MOD);
csub(g2[j], 1ll * pr[i] * pr[i] % MOD * sub(g2[k], sp2[i - 1]) % MOD);
}
cout << add(calc(n, 0), 1);
return 0;
}
```
### [SPOJ 34096 DIVCNTK - Counting Divisors (general)](https://www.spoj.com/problems/DIVCNTK/)
不难看出 $f(n)=d(n^k)$ 是积性函数，且 $f(p)=k+1$，$f(p^c)=ck+1$，上 Min_25 筛，令 $F=1$ 即可，后面再乘上 $k+1$。
### [Luogu P5493 质数前缀统计](https://www.luogu.com.cn/problem/P5493)
用 Min_25 筛的前半部分做 DP，边界 $g(n,0)$ 是自然数幂和，直接 $\mathcal{O}(k)$ 拉插即可。时间复杂度 $\mathcal{O}\left(\frac{n^{\frac{3}{4}}}{\log{n}}+\sqrt{n}k\right)$。
### [LOJ 6235 区间素数个数](https://loj.ac/p/6235)
令 $F(p)=1$，用 Min_25 筛前半部分直接 DP 即可。时间复杂度 $\mathcal{O}\left(\frac{n^{\frac{3}{4}}}{\log{n}}+\sqrt{n}\right)$。
### [SPOJ 22549 DIVFACT4 - Divisors of factorial (extreme)](https://www.spoj.com/problems/DIVFACT4/)
容易想到分别考虑每个 $[1,n]$ 中的质数 $p$，经典地，$n!$ 中 $p$ 因子的个数是
$$
\sum_{c=1}^{\log_p{n}}\left\lfloor\frac{n}{p^c}\right\rfloor
$$
因此
$$
d(n!)=\prod_{i=1}^{\pi(n)}\left(1+\sum_{c=1}^{\log_{p_i}{n}}\left\lfloor\frac{n}{p_i^c}\right\rfloor\right)
$$
注意到对于 $>\sqrt{n}$ 的 $p$，它们都只有 $1+\left\lfloor\frac{n}{p}\right\rfloor$ 的贡献，可以整除分块，这样对于一个块 $[l,r]$，计算出 $[l,r]$ 内的质数个数 $cnt$，给答案乘上 $\left(1+\left\lfloor\frac{n}{l}\right\rfloor\right)^{cnt}$ 即可。我们询问的端点都是形如 $\left\lfloor\frac{n}{i}\right\rfloor$ 的位置，用 Min_25 筛的前半部分即可解决。这部分时间复杂度为 $\mathcal{O}\left(\frac{n^{\frac{3}{4}}}{\log{n}}+\sqrt{n}\log{n}\right)$。
而对于 $\leq\sqrt{n}$ 的 $p$，我们可以 $\mathcal{O}(\log{n})$ 暴力跳，共有 $\pi\left(\lfloor\sqrt{n}\rfloor\right)=\mathcal{O}\left(\frac{\sqrt{n}}{\log{n}}\right)$ 个这样的 $p$，时间复杂度为 $\mathcal{O}(\sqrt{n})$。
由于模数 $m$ 很大，需要上快速乘。
总时间复杂度为 $\mathcal{O}\left(\frac{n^{\frac{3}{4}}}{\log{n}}+\sqrt{n}\log{n}\right)$。
代码：
```cpp
inline ll qmul(ll a, ll b) {
ull res = 1ull * a * b - (ull)((long double)a / m * b + 0.5L) * m;
return res < m ? res : res + m;
}
inline ll qpow(ll a, ll b) {
a %= m;
ll res = 1;
for (; b; b >>= 1) {
if (b & 1) res = qmul(res, a);
a = qmul(a, a);
}
return res;
}
inline void sieve(int n) {
v[1] = 1;
for (int i = 2; i <= n; ++i) {
if (!v[i]) v[i] = pr[++sz] = i;
for (int j = 1; j <= sz && i * pr[j] <= n; ++j) {
int x = pr[j];
v[i * x] = x;
if (v[i] == x) break;
}
}
}
inline int get_id(ll x) { return x <= sn ? id1[x] : id2[n / x]; }
inline ll calc_cnt(ll n, int p) {
ll res = 0;
for (; n >= p; n /= p) cadd(res, (n / p) % m);
return res;
}
int main() {
ios::sync_with_stdio(0), cin.tie(0);
sieve(N - 5);
cin >> t;
while (t--) {
cin >> n >> m, sn = sqrt(n), cnt = 0, ans = 1 % m;
for (ll l = 1, r, x; l <= n; l = r + 1) {
r = n / (x = n / l);
num[++cnt] = x, g[cnt] = x - 1;
if (x <= sn) id1[x] = cnt;
else id2[r] = cnt;
}
for (int i = 1; i <= sz && pr[i] <= sn; ++i)
for (int j = 1; j <= cnt && 1ll * pr[i] * pr[i] <= num[j]; ++j) {
int pos = get_id(num[j] / pr[i]);
g[j] -= g[pos] - (i - 1);
}
for (int i = 1; i <= sz && pr[i] <= sn; ++i) ans = qmul(ans, add(calc_cnt(n, pr[i]), 1));
for (ll l = sn + 1, r, x, c; l <= n; l = r + 1) {
r = n / (x = n / l);
ans = qmul(ans, qpow(x + 1, g[get_id(r)] - g[get_id(l - 1)]));
}
cout << ans << '\n';
}
return 0;
}
```
### [SPOJ 19975 APS2 - Amazing Prime Sequence (hard)](https://www.spoj.com/problems/APS2/)
对 Min_25 筛不那么裸的运用。
还是分开考虑质数和合数。对于质数，相当于求
$$
\sum_{i=1}^n[i\in\mathbb{P}]i
$$
这个 Min_25 筛随便做，令 $F(p)=p$，做前半部分的 DP 即可。
对于合数，我们考察 Min_25 筛前半部分的 DP：
$$
g(n,x)=g(n,x-1)-F(p_x)\left(g\left(\left\lfloor\frac{n}{p_x}\right\rfloor,x-1\right)-sp_{x-1}\right)
$$
注意到 $g\left(\left\lfloor\frac{n}{p_x}\right\rfloor,x-1\right)-sp_{x-1}$ 计算的就是 $n$ 以内 $\operatorname{lpf}=p_x$ 的数的个数，于是我们在 DP 过程中累加答案即可。
时间复杂度 $\mathcal{O}\left(\frac{n^{\frac{3}{4}}}{\log{n}}+\sqrt{n}\right)$。