【学习笔记】Min_25 筛
学长来讲 Min_25 筛了。
然后我咕了快一周了，终于写了！！
前情提要：没有复杂度证明。
Min_25 筛干啥的？
问 Min_25 啊
Min_25 筛是用来在 \(O(\frac{n^{\frac{3}{4}}}{\log n})\) 的复杂度内解决一些积性函数前缀和的问题。
先摆出这个积性函数的要求：
\(p\) 是质数，那么要求有：
- \(f(p)\) 是一个关于 \(p\) 的多项式。
- \(f(p^k)\) 便于计算。
符号规定
- \(p_i\) 代表第 \(i\) 个质数
- \(\mathrm{minp}(i)\) 代表 \(i\) 的最小质因子
Part I
首先，我们定义一个函数 \(S(x,j)\)：
既：\([2,x]\) 中最小质因子大于等于第 \(j\) 个质数的所有数的函数值的和。
我们先不考虑 \(1\)，那么发现答案就是 \(S(x,1) + 1\)。
然后把这东西展开：
其实就是把所有数分成了这三种类型：
- \(p_k\)：这个计算方法后面会讲。
- \(p_k^e\times \cdots\)：由于是积性函数，直接乘即可。剩下的数最小质因子就必须大于 \(p_k\) 了，所以是 \(S\left(\left\lfloor\dfrac{x}{p_k^e}\right\rfloor, k + 1\right)\)。
- \(p_k^e (e\ge 2)\)：直接计算即可。式子中用 \(e+1\) 的原因是防止重复计算 \(e=1\) 的情况。
然后大概讲一下中间两个 \(\sum\) 的限制条件：
- \(p_k^2 \le n\)：如果 \(p_k^2 > n\)，那么首先 \(p_k^{e+1} > n\)，这一项就不存在了；并且 \(p_k \times \cdots > p_k^2 > n\)，所以第一项也不存在，所以只需要计算 \(p_k^2 \le n\) 的情况就可以了。
- \(p_k^{e+1}\le n\)：当 \(p_k^{e+1}> n\) 时，形如 \(p_k^e \times \cdots\) 的数也一定大于 \(n\)，所以当 \(p_k^{e+1}\le n\) 时，两个式子都会卡到上界。
于是这一大堆式子就是 \(O(\frac{n^{\frac{3}{4}}}{\log n})\) 的。
我不会证
Part II
那么如何计算 \(\sum_{i\ge j,p_i\le n} f(p_i)\)？
首先我们可以转换一下，改成前缀和相减的形式。\(j\) 比较小，\(\sum_{\substack{i<j}} f(p_i)\) 是可以预处理出来的，所以我们只需要求出 \(\sum_{i\ge 1,p_i\le n} f(p_i)\) 即可。
我们再设一个函数 \(g(x,j)\)：
那么发现：当 \(j\rightarrow +\infty\) 时，就不存在 \(i\) 满足 \(\mathrm{minp}(i)>p_j\) 了，也就是：
由于 \(f(p)\) 是一个多项式，那么我们只需要把若干个 \(g_k(x,+\infty)\) 加起来，就是我们要求的答案。
我们可以考虑递归地求这个东西。
当 \(p_j^k > n\) 时，发现此时就不存在满足 \(\mathrm{minp}(i)>p_j\) 且不是质数的数了，所以： $$g_k(x,j)=g_k(x,j-1)$$
当 \(p_j^k \le n\) 时，我们考虑 \(\mathrm{minp}(i)>p_j\) 在 \(\mathrm{minp}(i)>p_{j-1}\) 的基础上多筛掉了哪些数。首先肯定 \(\mathrm{minp}(i)=p_j\)，那么因为 \(i^k\) 是完全积性函数，所以我们直接乘上 \(p_j^k\)，然后得到这样的式子：
因为 \(p_j\) 的指数不一定只有 \(1\)，所以这里就体现出完全积性函数的好处了：可以只提出来一个 \(p_j\)，后面可以随意乘，这样得到的数一定有 \(p_j\) 而且指数不一定。但是我们需要减去所有质数处的值，因为前后的质数都需要计算。
于是这个式子可以直接进行递推求，然后我们再带回原来的式子就好了。
Part III
代码实现：
#include <bits/stdc++.h>
using namespace std;
const int MAXN = 1000005, P = 1000000007;
long long n;
int sqr;
int pri[MAXN], pcnt, m;
bool vis[MAXN];
long long w[MAXN];
int g1[MAXN], g2[MAXN];
int p1[MAXN], p2[MAXN];
int mp[2][MAXN];
int qpow(int a, int b) {
int ans = 1;
while (b) {
if (b & 1) ans = 1ll * ans * a % P;
a = 1ll * a * a % P;
b >>= 1;
}
return ans;
}
void seive(int n) {
for (int i = 2; i <= n; i++) {
if (!vis[i]) {
pri[++pcnt] = i;
}
for (int j = 1; j <= pcnt && i * pri[j] <= n; j++) {
vis[i * pri[j]] = 1;
if (i % pri[j] == 0) break;
}
}
for (int i = 1; i <= pcnt; i++) p1[i] = (p1[i - 1] + pri[i]) % P;
for (int i = 1; i <= pcnt; i++) p2[i] = (p2[i - 1] + 1ll * pri[i] * pri[i] % P) % P;
}
int getId(long long x) { // 离散化
return x <= sqr ? mp[0][x] : mp[1][n / x];
}
const int INV2 = qpow(2, P - 2), INV6 = qpow(6, P - 2);
void getG(long long n) {
for (long long l = 1, r; l <= n; l = r + 1) { // 离散化
r = n / (n / l);
w[++m] = n / l;
if (w[m] <= sqr) mp[0][w[m]] = m;
else mp[1][n / w[m]] = m;
g1[m] = (w[m] % P * ((w[m] + 1) % P) % P * INV2 % P + P - 1) % P;
g2[m] = (w[m] % P * ((w[m] + 1) % P) % P * ((2 * w[m] % P + 1) % P) % P * INV6 % P + P - 1) % P;
}
for (int j = 1; j <= pcnt; j++) {
for (int i = 1; i <= m && 1ll * pri[j] * pri[j] <= w[i]; i++) {
int id = getId(w[i] / pri[j]);
g1[i] = (g1[i] - 1ll * pri[j] * (g1[id] - p1[j - 1] + P) % P + P) % P;
g2[i] = (g2[i] - 1ll * pri[j] * pri[j] % P * (g2[id] - p2[j - 1] + P) % P + P) % P;
}
}
}
int S(long long x, int j) {
if (x <= 1 || pri[j] > x) return 0;
int id = getId(x), ans = ((1ll * (g2[id] - g1[id]) - (p2[j - 1] - p1[j - 1])) % P + P) % P;
for (int k = j; k <= pcnt && 1ll * pri[k] * pri[k] <= x; k++) {
long long p1 = pri[k], p2 = p1 * p1;
for (int e = 1; p2 <= x; e++, p1 = p2, p2 *= pri[k]) {
ans = (ans + p1 % P * (p1 % P - 1) % P * S(x / p1, k + 1) % P + p2 % P * (p2 % P - 1) % P) % P;
}
}
return ans;
}
int main() {
scanf("%lld", &n);
sqr = sqrt(n);
seive(sqr);
getG(n);
printf("%d\n", (S(n, 1) + 1) % P);
return 0;
}
