前言
由于其由 Min_25 发明并最早开始使用，故称「Min_25 筛」。
- 从此种筛法的思想方法来说，其又被称为「Extended Eratosthenes Sieve」。
其可以在 \(O(\frac{n^{\frac{3}{4}}}{logn})\) 的时间复杂度下解决一类 积性函数 的前缀和问题。
要求： \(f(p)\) 是一个关于 \(p\) 的项数较少的多项式或可以快速求值； \(f(p^{c})\) 可以快速求值。
摘自 OI wiki
解析
例题
积性函数 \(f(p^k) = p^k(p^k - 1)\)
求：\(\sum_{i = 1} ^ n f(i)\)
准备 \(1\)
规定 \(p\) 代表质数, \(minp(i)\) 代表 \(i\) 的最小质因子（其中 \(minp(p) = p\)），则：
\(\sum_{i = 1} ^ n f(i) = \sum_{p \leq n} f(p) + \sum_{i \leq n \ and \ i \notin p} f(i)\)
这一步表示将答案拆成质数的答案和合数的答案。
那么对于合数部分：
根据积性函数的性质，则有：
\(f(i) = \prod_{k = 1} ^ m f(p_k ^ {c_k})\)。
提尽 \(i\) 的最小质因子 \(p_e\)，将 \(i\) 分成两部分：
\(p_e ^ {c_e} \times \frac{i}{p_e ^ {c_e}}\)
则上式可变为：
\(f(i) = f(p_e ^ {c_e}) \times f(\frac{i}{p_e ^ {c_e})})\)
那么合数部分的答案可转化为：
\(\sum_{p ^ e \leq n} f(p ^ e) \sum_{1 \leq i \leq n \ and \ p ^ e \mid i \ and \ minp(\frac{i}{p ^ e}) > p} f(\frac{i}{p ^ e})\)
代表的意思是：
枚举质数 \(p\) 及它的指数 \(e\)，
在可以提的合数 \(i\) 中提出 \(p ^ e\)，将合数 \(i\) 转化为 \(p ^ e \times \frac{i}{p ^ e}\)，即上式乘法分配律展开的结果。
然后再判断如下条件：
- \(p\) 是否是 \(i\) 的最小质因子
- \(p\) 是否提尽
这两个条件均可通过 \(minp(\frac{i}{p ^ e})\) 是否大于 \(p\) 来判断：
对于 \(minp(\frac{i}{p ^ e}) < p\) 的 \(i\)，那么不满足条件 \(1\)。
对于 \(minp(\frac{i}{p ^ e}) = p\) 的 \(i\)，那么不满足条件 \(2\)。
这些均可通过对 \(i\) 分解质因数来理解。
那么对于 \(f(p ^ e)\) 由题意可以直接求，对于 \(f(\frac{i}{p ^ e})\) 则是同样的求法，直到 \(\frac{i}{p ^ e}\) 是 \(p_k ^ {c_k}\) 可以直接求。
而上述的约束则保证了不会多求漏求。
准备 \(2\)
所以我们现在需要求出所有 \(f(p ^ k)\) 的值。
因为 \(n\) 的范围很大，所以我们没法线性筛求出 \(f\) 的前缀和。
考虑 DP。。。
设 \(h_k(i) = i ^ k\)，显然 \(h_k\) 是一个完全积性函数。
设 \(g_k(n,j) = \sum_{i = 1} ^ n [i \in Prime \ or \ minp(i) > p_j] h_k(i)\)，其中 \(p_j\) 表示第 \(j\) 个质数。
换句话说 \(g_k(n,j)\) 也就是 \(1\) ~ \(n\) 中质数或者最小质因数大于 \(p_j\) 的 \(i\) 的 \(h_k(i)\) 的前缀和。
\(g_k\) 的作用在下文，但也可以先抛开问题本身，看懂 \(g_k\) 的求法。
考虑转移。。。
对于 \(g_k(n, j - 1)\) 根据定义，其中包含两部分：
- \(g_k(n, j)\)
- 最小质因数等于 \(p_j\) 的 \(h_k(i)\) 的和。
所以可由 \(g_k(n, j - 1)\) 减去多余的部分得到 \(g_k(n, j)\)。
所以有状态转移方程如下：
\(g_k(n,j) = g_k(n, j - 1) - p_{j} ^ k (g(\frac{n}{p_{j} ^ k}, j - 1) - g(p_j, j - 1))\)
从状态转移方程可以看出：
最小质因数等于 \(p_j\) 的 \(h_k(i)\) 的和即为 \(g_k(n, j - 1) - p_{j} ^ k (g(\frac{n}{p_j}, j - 1) - g(p_{j - 1}, j - 1))\)。
原因如下：
根据定义不难得到 \(g_k(p_j, j - 1)\) 即为前 \(j - 1\) 个质数的 \(k\) 次方和。
而对于 \(g_k(\frac{n}{p_j}, j - 1)\)，包含了两部分：
- \(1\) ~ \(p_{j - 1}\) 的质数和
- 最小质因数大于等于 \(p_j\) 的合数和，及 \(p_{j}\) ~ \(\frac{n}{p_j}\) 的质数和
注意第二点，如果我们将这些最小质因数大于等于 \(p_j\) 的合数或者后一部分的质数同时乘上 \(p_j\) 那么我们得到的是不是都是最小质因数等于 \(p_j\) 的合数，这点应该很好理解。
那也就是说得到了我们想要筛去的合数。
因为我们只知道上述两部分的和 \(g_k(\frac{n}{p_j}, j - 1)\)，所以在乘的过程中第一部分也会被乘上一个 \(p_j\)，那么我们应该把这一部分质数减去。
也就是 \(g_k(p_j, j - 1)\)。
又因为 \(h_k\) 是完全积性函数，而我们谈论的和是 \(h_k\) 的累加，根据乘法分配律，我们可以对 \(g_k\) 运用完全积性函数的性质（\(f(ab) = f(a)f(b)\)）。
那么要得到最小质因数等于 \(p_j\) 的 \(h_k(i)\) 的和，我们可以直接给我们的第二部分乘上 \(h_k(p_j)\) 的值，也就是 \(p_j ^ k\)，得到我们要求的函数值。
所以关于状态转移方程的解释就结束了。
关于初始化：
显然对于 \(g_k(n, 0) = \sum_{i = 1} i ^ k\)
求解
而我们得到的 \(g_k\) 有什么用呢？
设 \(g(n, j) = \sum_{i = 1} ^ n [i \in Prime \ or \ minp(i) > p_j] f(i)\)。
因为 \(f(p) = p(p - 1) = p ^ 2 - p\)。
那么将满足条件的 \(f(i)\) 累加，就得到两部分：
- \(i ^ 2\) 的和
- \(i\) 的和的相反数
如果令上文的 \(k\) 分别等于 \(1, 2\)。
那么就有： \(g(n, j) = g_2(n, j) - g_1(n, j)\)。
如果设 \(p_x\) 为小于等于 \(\sqrt{n}\) 的最大的质数，那么 \(1\) ~ \(n\) 的质数的 \(f\) 的和 \(g(n, x)\) 记为 \(g(n)\)。
设 \(S(n, j) = \sum_{minp(i) > p_j} f(i)\)。
显然我们也需要 DP 。。。
同样的，考虑 \(S(n, j)\) 的组成：
- 质数部分：大于 \(p_j\) 小于等于 \(n\) 的质数
- 合数部分
设 \(sum_k(x)\) 表示小于等于 \(p_x\) 的质数的 \(k\) 次方的和。
那么质数部分就等于：
\(g(n) - （sum_2(j) - sum_1(j)）\)
对于合数部分：
我们模仿 \(g_k\) 的求法
因为 \(f\) 不是完全积性函数，而是积性函数，我们不能像求 \(g_k\) 那样只提一个 \(p_j\)。
所以我们要把 \(p_j\) 提尽。
得到：
\(\sum_{k > j \ and \ p_k ^ e \leq n} f(p_k ^ e)(S(\frac{n}{p_k ^ e} + [e \neq 1]))\)
利用 \(S(\frac{n}{p_k ^ e}\) 的思想和求 \(g_k\) 的时候完全一样，这里就不在赘述了，实在不懂的同学可以手推一下。
后面加上一个 \([e \neq 1]\) 是因为如果 \(e > 1\)，那么\(p ^ e\) 在我们提 \(p_k ^ e\) 时没有算进去，而 \(e = 1\) 时，在质数部分就算了，就不用了算进去。
下标的离散化
因为 \(n = 1e10\)，直接开肯定不行。
由性质 \(\lfloor \frac{\lfloor \frac{x}{a} \rfloor}{b} \rfloor = \lfloor \frac{x}{ab} \rfloor\)，在递归的时候，无论我们把 \(n\) 除以几，得到的都是 \(n\) 除以某个数，所以我们可以直接只保存 \(n / x\) 的值。
来自 wucstdio 大佬
代码
#include<cstdio>
#include<cmath>
#include<cstring>
using namespace std;
typedef long long LL;
const LL mod = 1e9 + 7, inv2 = 500000004, inv3 = 333333336;
const int N = 1e5 + 5;
LL n, Prime[N + 5], num, sp1[N + 5], sp2[N + 5], tot = 0, w[3 * N], g1[3 * N], g2[3 * N], sqr;
int id1[N + 5], id2[N + 5];
bool notPrime[N + 5];
void init() {
memset(notPrime, false, sizeof(notPrime));
num = 0;
for(int i = 2; i <= N; i++) {
if(!notPrime[i]) {
Prime[++num] = i;
sp1[num] = (sp1[num - 1] + i) % mod;
sp2[num] = (sp2[num - 1] + 1ll * i * i % mod) % mod;
}
for(int j = 1; j <= num && Prime[j] * i <= N; j++) {
notPrime[i * Prime[j]] = true;
if(i % Prime[j] == 0) break;
}
}
}
LL s(LL x, int y) {
if(Prime[y] >= x) return 0;
LL k = x <= sqr ? id1[x] : id2[n / x];
LL ans = (g2[k] - g1[k] + mod - (sp2[y] - sp1[y]) + mod) % mod;
for(int i = y + 1; i <= num && Prime[i] * Prime[i] <= x; i++) {
LL P = Prime[i];
for(int e = 1; P <= x; e++, P = P * Prime[i]) {
LL xx = P % mod;
ans = (ans + xx * (xx - 1) % mod * (s(x / P, i) + (e != 1))) % mod;
}
}
return ans % mod;
}
int main() {
init();
scanf("%lld", &n);
sqr = sqrt(n);
for(LL l = 1, r; l <= n; l = r + 1) {
r = (n / (n / l));
w[++tot] = n / l;
g1[tot] = w[tot] % mod;
g2[tot] = g1[tot] * (g1[tot] + 1) / 2 % mod * (2 * g1[tot] + 1) % mod * inv3 % mod - 1;
if(g2[tot] < 0) g2[tot] += mod;
g1[tot] = g1[tot] * (g1[tot] + 1) / 2 % mod - 1;
if(g1[tot] < 0) g1[tot] += mod;
if(w[tot] <= sqr) id1[w[tot]] = tot;
else id2[n / w[tot]] = tot;
}
for(int i = 1; i <= num; i++)
for(int j = 1; j <= tot && Prime[i] * Prime[i] <= w[j]; j++) {
LL k = w[j] / Prime[i];
k = k <= sqr ? id1[k] : id2[n / k];
g1[j] -= Prime[i] * (g1[k] - sp1[i - 1] + mod) % mod;
g2[j] -= Prime[i] * Prime[i] % mod * (g2[k] - sp2[i - 1] + mod) % mod;
g1[j] %= mod, g2[j] %= mod;
if(g1[j] < 0) g1[j] += mod;
if(g2[j] < 0) g2[j] += mod;
}
printf("%lld", (s(n, 0) + 1) % mod);
return 0;
}