min25筛 和 洲阁筛
min_25 筛 和 洲阁筛
min25 筛
已知一积性函数 \(f(x)\)，\(f(p)\) 为关于 \(p\) 的多项式，且 \(f(p^k)\) 可以快速算出来，求 \(\sum_{i=1}^nf(i)\)。
例：求 \(\sum_{i=1}^n \sigma_0(i^k)\)，其中 \(n = 10^{10}\)
step1 : 求质数的积性函数前缀和
求 \(1 ... n\) 以内的质数个数。
我们设 \(f(n,m)\) 表示 \(2...n\) 中的所有质数或者不含最小质因子为 \(p_1...p_m\) （\(p_i\) 表示第 \(i\) 个质数）的数的个数（即最小质因子大于 \(p_m\)）
那么有：
意义为：考虑从 \(f(n,m-1)\) 处转移过来，那么我们只需要筛掉最小质因子恰好是 \(p_m\) 的数。那么我们把所有要筛的数提出来个 \(p_m\)，剩下的数的最小质因子为 \(p_m,p_{m+1},...\)。注意删掉“质数"部分。
注意，我们并不想筛掉 \(p_m\)，但是注意到 \(f\) 是不包括 \(1\) 的，所以后面的部分并没有把 \(1\) 筛掉，符合要求。
更普遍的情况。
类似地，我们设 \(f(n,m)\) 表示 \(2...n\) 中的所有质数或者不含最小质因子为 \(p_1...p_m\) （\(p_i\) 表示第 \(i\) 个质数）的数的函数和。这时候我们需要认为合数也有一个假的和质数一样的函数和。那么有：
（函数重名了，用 \(F(n)\) 代替了积性函数；\(psum\) 为预处理出来的质数函数前缀和）
由于是多项式，我们可以拆成 \(F(p) = p^k\) 的形式，最终拼起来。这样，假函数成为完全积性函数，直接提出来 \(F(p_m)\) 就是对的了。
- 复杂度证明
见那篇博客。
发现 \(p_m^2 \le n\)
可以积分证明，复杂度为 \(O(n^{0.75} / \log n)\)
step2 : 求所有数的积性函数和
我们设 \(g(n,m)\) 表示 \(1...n\) 中最小质因子大于 \(p_m\) 的数的（真）积性函数和。
那么有：
意义为：分质数，合数考虑贡献。质数的贡献我们已经搞好了，为 \(f(n,\infty) - psum(m)\)，合数的话我们需要枚举最小质因子及其次数，即枚举 \(k > m\) 的 \(p_k\)，再枚举次数 \(e\)。由于是积性函数，我们可以直接将 \(F(p_k^e)\) 提出来。值得注意的是，形如 \(p^k,k > 1\) 的数我们既没有在质数处枚举到，也没有在合数处提出来 \(p_k^e\) 后枚举到，所以需要加上个 \([e \not= 1]\)。
复杂度：不会证。状态数和上面一样是 \(O(n^{0.75} / \log n)\) 的，但是转移复杂度比较玄学，据说总复杂度是 \(O(n^{0.75} / \log n)\) 的，也有的说是 \(O(n^{1 - \epsilon})\) 的。
洲阁筛
功能基本相同，但是跑的没 min25 筛快。不过可以一次性求出所有 \(\lfloor \frac{n}i \rfloor\) 的前缀和的值。
（好像和正常的洲阁筛不太一样，这个是 hs-black 教我的）
step1:求质数的积性函数前缀和
和 min25 筛完全一致，直接把 \(f(n,m)\) 拿过来即可。
step2 : 求所有数的积性函数和
状态不太一样，有点像 step1 中的那个 \(f\) 数组。
定义 \(g(n,m)\) 表示 \(2...n\) 中质数或最小质因子大于 \(p_m\) 的所有数的真的积性函数和。那么有：
意义为：
对于不包含 \(p_m\) 的那部分，直接从 \(g(n,m)\) 那里转移过来；对于包含 \(p_m\) 的那部分，枚举指数 \(c\)，如果形如 \(p_m^c\)，那么答案就是 \(F(p_m^c)\)；否则还会有若干比 \(p_m\) 大的质因子，答案为 \(F(p_m^c)(g(\lfloor \frac{n}{p_m^c} \rfloor,m) - psum(m))\)。
注意到当 \(p_m^2 > n\) 的时候 \(g(n,m) = f(n,m)\)，于是将 \(g\) 初始化设为 \(f\)，按照 \(m\) 从大到小递推即可。最终答案为 \(g(n,0)\)。好处是推完后我们得到的是一个数组，包含了所有 \(g(\lfloor \frac{n}{i} \rfloor,0)\)，于是再求 \(\lfloor \frac{n}{i} \rfloor\) 的时候就不用重新求了。
代码实现为 \(g(n,m-1) = g(n,m) + \sum_c F(p_m^c)(g(\lfloor \frac{n}{p_{m}^c} \rfloor,m) - psum(m)) + F(p_m^{c+1})\)，主要是为了方便减少枚举上界。
- 复杂度
不知道，大家说还是 \(O(n^{0.75} / \ln n)\) 的，但是总是感觉多枚举个 \(c\) 可能会增大复杂度。但是可以证明复杂度是不超过 \(O(n^{0.75})\) 的。
应用
最常见的当属各种积性函数求和题，不过基本上只是在套板子。
还有一些考察对 min25 筛的灵活应用。
质数计数
求 \(\le n\) 的质数中模 \(mod\) 等于 \(r\) 的数有几个。\(n \le 10^{10},r \le 13\)
对min25筛第一步的考察。
可以在第一个DP数组中新加一维，\(f(n,m,r)\) 表示 \(2...n\) 中不包含前 \(m\) 个质数或者本身就是质数的数个数，转移还是去除掉 \(p_m\):
次大质因子求和
求 \(\sum_{i \le n} f(i)\)，其中 \(f(i)\) 表示 \(i\) 的非严格次大质因子（除掉最大质因子后的最大质因子）。如果没有则为 \(0\)。\(n \le 10^{11}\)。
对 min25 筛第二步的考察。
考虑 min25 筛第二步实质上是在对每个数的各个质因子从大到小筛，其中最大的纯质数单独处理。我们想要筛出次大质因子，就在筛的时候钦定当前考虑的这个因子就是次大质因子（不是的话可以继承），剩下的比它还大的数就只能是单个的质数了，快速质数计数即可。
代码
min25 筛
//divcntk
int block, itot, pos[N];
ull id[N], n, K, f[N];
inline int get_id(ull x) {
return pos[x <= block ? x : block + n / x];
}
inline void get_f() {
for (ull l = 1, r; l <= n; l = r + 1) {
ull zhi = n / l;
r = n / zhi;
id[++itot] = zhi;
pos[zhi <= block ? zhi : block + n / zhi] = itot;
f[itot] = zhi - 1;
}
for (int j = 1; 1ull * pri[j] * pri[j] <= n; ++j) {
for (int i = 1; id[i] >= 1ull * pri[j] * pri[j]; ++i) {
f[i] = f[i] - f[get_id(id[i] / pri[j])] + j - 1;
}
}
}
ull get_g(ull n, int j) {
if (pri[j] >= n) return 0;
ull res = (f[get_id(n)] - j) * (K + 1);//bug
for (int k = j + 1; 1ull * pri[k] * pri[k] <= n; ++k) {
ull tmp = pri[k];
for (int e = 1; tmp <= n; ++e, tmp = tmp * pri[k]) {
res += (K * e + 1) * (get_g(n / tmp, k) + (e > 1));
}
}
return res;
}
inline void work() {
read(n), read(K); block = sqrt(n);//bug
get_f();
ull res = get_g(n, 0) + 1;
printf("%llu\n", res);
}
洲阁筛
//divcnt3
inline void init(const int up) //筛素数
inline int Id(ll x) { return x <= S ? x : tot - n / x + 1; }
ll r[N], tot, g[N];
inline ll calc(ll k) { return 3 * k + 1; }//h(p^k)
inline void work() {
read(n);
S = sqrt(n); init(S + 200);
for (ll l = 1; l <= n; l = r[tot] + 1)
r[++tot] = n / (n / l), g[tot] = r[tot] - 1;
int ptr = 1;
for (int j = 1; pri[j] * pri[j] <= n; ++j) {
while (pri[j] * pri[j] > r[ptr]) ++ptr;
for (int i = tot; i >= ptr; --i) {
g[i] -= g[Id(r[i] / pri[j])] - j + 1;
}
}
for (int i = 1; i <= tot; ++i) g[i] *= 4;
ptr = tot + 1;
for (int j = pcnt; j; --j) {
while (pri[j] * pri[j] <= r[ptr - 1]) --ptr;
for (int i = tot; i >= ptr; --i) {
ll nw = pri[j], nww = nw * nw;
for (int c = 1; nww <= r[i]; ++c, nw = nww, nww *= pri[j]) {
g[i] += calc(c) * (g[Id(r[i] / nw)] - 4 * j) + calc(c + 1);//Attention! 略有修改
}
}
}
ll res = g[tot] + 1;
printf("%lld\n", res);
}
- 注意：
分清 i 和 r[i]！
