Min_25 筛
前置知识：积性函数，数论分块。
1 概述
前面我们学习过杜教筛，它可以在亚线性的复杂度内求出一些数论函数的前缀和。而 Min_25 筛同样可以在亚线性复杂度内求出一些积性函数的前缀和。
考虑我们朴素的筛法求积性函数值无法低于线性复杂度的原因就是其必须枚举区间内的每一个数，而 Min_25 筛则是将区间内的所有数分为质数和合数两类进行处理，然后相加得出答案。于是 Min_25 筛的使用对于求解的积性函数 \(f(x)\) 有一定要求：对于质数 \(p\)，\(f(p)\) 可以看作一个关于 \(p\) 的低阶多项式，并且它的值在质数的幂次方处可以快速求出。
2 算法思想
2.1 分类
上面说过，Min_25 筛的重点在于将所有数分为质数和合数两类，也就是：
然后进一步的，我们枚举后面合数的最小质因子以及质因子的次数，可以得到：
其中 \(\text{minp}(i)\) 表示 \(i\) 的最小质因子。注意一个合数的最小质因子一定小于 \(\sqrt n\)，所以我们的 \(p\) 只需要枚举到 \(\sqrt n\) 即可。
这个分类思想是 Min_25 筛中一个重要的思想，在最后求解答案时会再次用到。
2.2 质数求解
由于 Min_25 筛求解的 \(n\) 很大，所以肯定不能直接筛出 \(f\) 的和。 我们假设对于质数 \(p\)，\(f(p)\) 可以写成 \(f'(p)=\sum a_i p^i\) 的多项式形式，那么我们现在的目标就是统计出 \(\sum p^k\) 的值。
考虑利用埃氏筛的思路，同时考虑一个 dp，设 \(g(n,j)\) 表示 \(n\) 范围内经过 \(j\) 轮埃氏筛后剩下数字的 \(k\) 次方之和。再说直白点就是所有质数或者最小质因子大于第 \(j\) 个质数 \(p_j\) 的所有数的 \(k\) 次方之和。则最后要求的就是 \(g(n,j)\)，其中 \(p_j\) 是最后一个 \(\le \sqrt n\) 的质数。
考虑怎样转移，我们考虑从 \(g(n,j-1)\) 转移到 \(g(n,j)\)，此时需要减去最小质因子是 \(p_j\) 的所有合数的贡献。
既然最小质因子是 \(p_j\)，那么我们直接提出一个 \(p_j\)，剩下数的范围就是 \(\lfloor\tfrac{n}{p_j}\rfloor\)，为了保证最小质因子是 \(p_j\)，我们需要保证除掉 \(p_j\) 后的这些数的最小质因子大于 \(p_{j-1}\)。于是剩下的部分的贡献就是 \(g(\lfloor\tfrac n{p_j}\rfloor,j-1)\)，但是此时我们可能会把 \([1,p_{j-1}]\) 范围内的质数的贡献也减掉，所以我们需要把 \(g(p_{j-1},j-1)\) 加回来。于是可以得到转移方程：
这里我们能直接提出 \(p_j\) 的原因就在于幂函数是一个完全积性函数，所以不必在意是否互质，直接相乘即可。
不过现在直接暴力求解的复杂度仍然没有降下去。考虑中间的一项 \(g(\lfloor\tfrac{n}{p_j}\rfloor,j-1)\)。由数论分块的结论 \(\left\lfloor\dfrac{\left\lfloor\tfrac{a}{b}\right\rfloor}{c}\right\rfloor=\left\lfloor\dfrac{a}{bc}\right\rfloor\) 容易发现，上述式子只会用到所有形如 \(\lfloor\tfrac nx\rfloor,x\le n\) 位置上的 dp 值，而这些值最多只有 \(2\lfloor\sqrt n\rfloor\) 种。所以只需要求出它们的值，然后直接 dp 即可。
不过式子中还有一项是 \(g(p_{j-1},j-1)\)，也就是质数处的取值。这一部分可以直接线性筛预处理出来，不过实际上没有必要，因为用于转移的 \(p_j\) 一定 \(<\lfloor\sqrt n\rfloor\) ，所以它们一定可以表示为 \(\lfloor\tfrac nx\rfloor,x\le n\) 的形式，因此可以直接调用之前求出的 \(g\) 值。
此处有一个命题：\(\forall k< \lfloor\sqrt n\rfloor,\exists x,k=\lfloor \tfrac nx\rfloor\)。
考虑反证，设存在一个 \(i\)，满足 \(\lfloor\tfrac{n}{i+1}\rfloor<k<\lfloor\tfrac ni\rfloor\)，则 \(k\) 不可被表示。于是有 \(\tfrac {n}{i+1}<k\)，即 \(n<k(i+1)\)。代回去后可得 \(k<\lfloor\tfrac ni\rfloor<\lfloor\tfrac {k(i+1)}i\rfloor=k+\lfloor\tfrac ki\rfloor\)。
另一方面来讲，由于 \(\tfrac{n}{i+1}<k<\sqrt n\)，所以 \(\tfrac n{i+1}<\sqrt n\)，即 \(i+1>\sqrt n\)。由于 \(k<\lfloor\sqrt n \rfloor\)，所以 \(k\le \sqrt n-1\)，于是 \(i>\sqrt n -1\ge k\)，所以 \(\lfloor\tfrac ki\rfloor=0\)。
综合可得 \(k<k\)，显然不成立，于是原命题成立。
注意这里的 \(g\) 和下文中的 \(S\) 都没有计算 \(1\) 的贡献，不过这对于 \(g\) 没有影响。
2.3 合数求解
2.3.1 方法一
和上面类似的，设 \(S(n,j)\) 表示 \(n\) 范围内最小质因子大于 \(p_j\) 的所有数的 \(f\) 值之和。接下来令 \(G(n)\) 表示 \(n\) 范围内所有质数的 \(f\) 值之和，显然这可以用上面的 \(g(n,j)\) 简单求出。
考虑按照最开始分类时的思路，将这个贡献拆成质数的贡献和合数的贡献。质数贡献容易计算，就是 \(G(n)-G(p_j)\)，合数的贡献依然考虑枚举最小质因子和次数。于是有以下转移式：
注意枚举的 \(p_k\le \lfloor\sqrt n\rfloor,p_k^e\le n\)。上文说过，\(S\) 是不计算 \(1\) 的贡献的，那么 \(p_k^e\) 的贡献我们就统计不上了，需要另外加上去。不过当 \(e=1\) 的时候 \(p_k\) 就是个质数，而质数贡献在前面是统计过的，因此不用额外加上去。
求解的话直接暴力递归求解即可。
2.3.2 方法二
我们可以直接套用求解质数的状态，设 \(S(n,j)\) 表示 \(n\) 范围内所有质数或最小质因子 \(>p_j\) 的数的 \(f\) 值之和。我们从 \(S(n,j+1)\) 推到 \(S(n,j)\)。思路和求解质数是完全一致的，考虑最小质因子为 \(p_{j+1}\) 的数字个数即可。不过由于此时的 \(f\) 是积性函数，所以要一次性将所有 \(p_{j+1}\) 全部提出来才行。于是有以下转移：
枚举时依然需要保证 \(p_{j+1}^e \le n\)。这里加 \([e\neq 1]\) 的理由和上文一致，而要加上取 \(\min\) 操作的原因是我们无法保证 \(p_{j+1}^{e+1}\) 与 \(n\) 的大小关系，而求解质数时我们一定有 \(p_j^2\le n\)，所以不必取 \(\min\)。
但是再细想一下就会发现不对，如果 \(p_{j+1}^{e+1} >n\)，那么 \(\lfloor\tfrac{n}{p_{j+1}^e}\rfloor\) 就一定小于 \(p_{j+1}\)，于是求出的 \(S\) 中就不存在合数了；而剩下的质数的函数值之和又恰好是 \(G\) 的值，两者正好抵消，没有贡献。于是这一项的贡献就只剩下了一个 \(f(p_{j+1}^{e})\)。而如果我们将 \(f(p_{j+1}^e)[e\neq 1]\) 的贡献提前到 \(e-1\) 处计算后，这一项的贡献就彻底没了。于是我们可以简化成：
此时需要保证 \(p_{j+1}^{e+1}\le n\)。这样的话设初值 \(S(n,+\infty)=G(n)\)，我们就可以直接递推求解了。
不管采用什么方法，可以证明，Min_25 筛的总复杂度是 \(O(\dfrac{n^{0.75}}{\log n})\) 的。不过值得注意的是，方法一常数较小；而方法二不仅能求出 \(S(n)\)，也能求出所有的 \(S\left(\lfloor\tfrac{n}{i}\rfloor\right)\)，在某些题目中比较有用。
2.4 代码
模板题：【模板】Min_25 筛。题目中已经告诉我们 \(f(p^k)=p^k(p^k-1)=(p^k)^2-(p^k)\)，也就是一个二次多项式。直接套模板即可。采用方法一代码如下：
#include <bits/stdc++.h>
#define int long long
using namespace std;
const int Maxn = 2e5 + 5;
const int Inf = 2e9;
const int Mod = 1e9 + 7;
const int Inv2 = 500000004;
const int Inv6 = 166666668;
int n, m;
int prim[Maxn], cnt, vis[Maxn];
void init(int N) {//预处理质数
for(int i = 2; i <= N; i++) {
if(!vis[i]) {
prim[++cnt] = i;
}
for(int j = 1; prim[j] * i <= N; j++) {
vis[i * prim[j]] = 1;
if(i % prim[j] == 0) break;
}
}
}
int ind1[Maxn], ind2[Maxn];
//分别存储 x 和 n/x 在 dp 数组的编号，这里 x <= sqrt(n)
int g1[Maxn], g2[Maxn], val[Maxn], tot;
//一次项和二次项之和 所有 n/x 的值
int getid(int x) {return x <= m ? ind1[x] : ind2[n / x];}//获取当前值的编号
int F(int x) {x %= Mod; return (x * x % Mod - x + Mod) % Mod;}
int S(int n, int j) {//递归求解 S(n,j)
if(prim[j] > n) return 0;//递归边界，直接返回
int id1 = getid(n), id2 = getid(prim[j]);
int res = (g2[id1] - g1[id1]) - (g2[id2] - g1[id2]);
res = (res % Mod + Mod) % Mod;//算出 G(n)-G(p[j])
for(int k = j + 1; k <= cnt && prim[k] <= n / prim[k]; k++) {
int num = prim[k];
for(int e = 1; num <= n; e++, num = num * prim[k]) {
res = (res + F(num) * (S(n / num, k) + (e > 1)) % Mod) % Mod;
//按照递归式直接递归即可
}
}
return res;
}
signed main() {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> n;
m = sqrt(n);
init(m);
int l = 1, r;
while(l <= n) {//数论分块求出所有值
r = n / (n / l);
val[++tot] = n / l;
if(val[tot] <= m) ind1[val[tot]] = tot;
else ind2[n / val[tot]] = tot;
int num = val[tot] % Mod;//注意要先取模
g1[tot] = num * (num + 1) % Mod * Inv2 % Mod - 1;
g2[tot] = num * (num + 1) % Mod * (2 * num + 1) % Mod * Inv6 % Mod - 1;
//先求出所有 g(x,0)，显然就是前缀和
//减一是为了减掉 1 的贡献
l = r + 1;
}
for(int j = 1; j <= cnt; j++) {//先枚举第二维
for(int i = 1; i <= tot && prim[j] <= val[i] / prim[j]; i++) {//枚举第一维
int id1 = getid(val[i] / prim[j]), id2 = getid(prim[j - 1]);
g1[i] -= prim[j] * (g1[id1] - g1[id2]) % Mod;
g2[i] -= prim[j] * prim[j] % Mod * (g2[id1] - g2[id2]) % Mod;
g1[i] = (g1[i] % Mod + Mod) % Mod;
g2[i] = (g2[i] % Mod + Mod) % Mod;//按照式子求解
}
}
cout << (S(n, 0) + 1) % Mod << '\n';//求出答案之后记得加一
return 0;
}
采用方法二代码如下：
#include <bits/stdc++.h>
#define int long long
using namespace std;
const int Maxn = 2e5 + 5;
const int Inf = 2e9;
const int Mod = 1e9 + 7;
const int Inv2 = 500000004;
const int Inv6 = 166666668;
int n, m;
int prim[Maxn], cnt, vis[Maxn];
void init(int N) {
for(int i = 2; i <= N; i++) {
if(!vis[i]) prim[++cnt] = i;
for(int j = 1; prim[j] * i <= N; j++) {
vis[i * prim[j]] = 1;
if(i % prim[j] == 0) break;
}
}
}
int ind1[Maxn], ind2[Maxn];
int val[Maxn], tot, g1[Maxn], g2[Maxn];
int S[Maxn], G[Maxn];
int getid(int x) {return x <= m ? ind1[x] : ind2[n / x];}
int F(int x) {x %= Mod; return (x * x % Mod - x + Mod) % Mod;}
signed main() {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> n;
m = sqrt(n);
init(m);
int l = 1, r;
while(l <= n) {
r = n / (n / l);
val[++tot] = n / l;
if(val[tot] <= m) ind1[val[tot]] = tot;
else ind2[n / val[tot]] = tot;
int num = val[tot] % Mod;
g1[tot] = num * (num + 1) % Mod * Inv2 % Mod - 1;
g2[tot] = num * (num + 1) % Mod * (2 * num + 1) % Mod * Inv6 % Mod - 1;
l = r + 1;
}
for(int j = 1; j <= cnt; j++) {
for(int i = 1; i <= tot && prim[j] * prim[j] <= val[i]; i++) {
int id1 = getid(val[i] / prim[j]), id2 = getid(prim[j - 1]);
g1[i] -= prim[j] * (g1[id1] - g1[id2]) % Mod;
g2[i] -= prim[j] * prim[j] % Mod * (g2[id1] - g2[id2]) % Mod;
g1[i] = (g1[i] % Mod + Mod) % Mod, g2[i] = (g2[i] % Mod + Mod) % Mod;
}
}
//到这里都与方法一一致
for(int i = 1; i <= tot; i++) {
S[i] = G[i] = (g2[i] - g1[i] + Mod) % Mod;
}
for(int j = cnt; j >= 1; j--) {
for(int i = 1; i <= tot && prim[j] * prim[j] <= val[i]; i++) {
int num = prim[j];
for(int e = 1; num <= val[i] / prim[j]; e++, num = num * prim[j]) {
int id1 = getid(val[i] / num), id2 = getid(prim[j]);
(S[i] += F(num) * (S[id1] - G[id2] + Mod) % Mod + F(num * prim[j])) %= Mod;//按照递推式直接求 S
}
}
}
cout << (S[1] + 1) % Mod << '\n';
return 0;
}
3 例题
例 1 [LOJ6235] 区间素数个数
很简单，我们只需要求出 \(n\) 范围内所有质数的 \(0\) 次方和即可，只用做第一步质数求解即可得出最后答案。
例 2 [LOJ6053] 简单的函数
实际上还有 [LOJ6783 ~ 6785] 简单的函数 加强版，\(\text{Link}\)。
这个函数很有意思，首先它满足在质数的幂处可以简单求解，不过在 \(p\) 处的值是 \(p\oplus 1\)，看上去不像一个多项式。不过我们发现，除了 \(2\) 以外的质数都是奇数，也就是除了 \(2\) 以外的 \(p\) 的 \(f(p)\) 都是 \(p-1\)，而 \(f(2)=2+1\)。
不过我们完全可以将 \(f(2)\) 也看成 \(2-1\)，最后特判一下再加上 \(2\) 即可。然后这个函数在 \(p\) 处就可以转化为多项式 \(p-1\)，维护一次项和与零次项和即可得出答案。
后面三道加强版的题就需要用方法二求出所有 \(S\) 的值了，最后一个题还需要卡常。
浙公网安备 33010602011771号