Min_25筛
upd on 11.20：好像不少锅，懒得修了。
一个神奇的筛法，复杂度是神奇而且常数贼小的\(O(\frac {n^{\frac 34}}{\log n})\)。当然是拿来求积性函数前缀和的东西。一般拿来解决\(10^{10}\)级别的前缀和处理。
首先对于我们\(O(n^{\frac 23})\)的杜教筛，这个东西的泛用性非常狭窄，也就是必须找两个合适的函数卷积起来才能搞，如果你找不着函数的话那就没治，比如说\(\pi(n)\)（\(n\)以内素数个数）这种东西（当然这个有个\(O(\frac {n^{\frac 23}}{\log^2n})\)的玄妙做法但我看不懂）。所以这时候就需要一些别的亚线性复杂度的筛法（min25，pn，洲阁）来搞。
然后是Min_25筛的适用条件：
- 是积性函数
- 在素数处的取值是个项数较少的多项式
- \(f(p^k)\)的值能快速取得
就这仨（一般来说是的，但是还有二般，比如LOJ572这个简单莫反+可能还行的杜教筛+sbMin_25筛的题）。下面拿洛谷板子举例子。
题面：有积性函数\(f(x)\),满足：\(f(1)=1,f(p^k)=p^k(p^k-1)\)，求\(\sum_{i=1}^nf(i) \mod 1000000007\)的值，\(n\le 10^{10}\)。
首先第一步我们把质数取值和合数取值分开。
分开看两个部分。
第一部分
于是我们现在要求所有质数取值处的前缀和。发现正常搞绝对不行，所以考虑一个人类智慧的dp。为了方便，我们把\(f(i)\)看成一个完全积性函数，要求这个完全积性函数在质数处的取值和原来的函数一样。如果原来函数是完全积性函数的话直接搞，如果不是就拆开成若干完全积性函数然后加起来。比如这个题就是\(f(p^k)=(p^k)^2-p^k\)。要求质数取值一样，直接设个\(f_1(i)=i,f_2(i)=i^2\)。之后一大段这个完全积性函数图方便直接用\(f\)表示了。
设
可能有一大群看不懂的表示，解释一下。\(P\)是\(n\)以内的素数集，\(minn(i)\)是\(i\)最小的质因子，\(p_j\)是第\(j\)个质数。于是我们找到最大的小于\(\sqrt n\)的质数\(p_x\)，我们想要的就是\(g(n,x)\)。
考虑这个东西怎么转移。第一维一会再说，带个\(n\)不太好预处理。先看第二维，这个最大就是\(\sqrt n\)的。考虑如何从\(g(n,j-1)\)转移到\(g(n,j)\)。我们发现，如果这个\(j\)增大了那会减少一些数，这些数就是最小质因子是\(p_j\)的合数。
考虑如何把这些数用\(g\)表示出来。发现它们就是所有的数减掉被拿出来的质数。第一部分我们把所有\(p_j\)的倍数除掉一个\(p_j\)因数，剩下的就是\(\lfloor \frac n{p_j} \rfloor\)，这些数中最小质因数大于\(p_{j-1}\)的就是要去掉的数(因为去掉了一个最小质因数\(p_j\)，如果仍然存在比它小的质因数那么最小质因数就不是\(p_j\))，说人话就是\(g(\lfloor \frac n{p_j} \rfloor,j-1)\)。然后被拿出来的质数就是小于\(p_{j-1}\)的所有质数，就是\(g(p_{j-1},j-1)\)。由于我们去掉了一个\(p_j\)因子的贡献，而\(f\)又是完全积性函数，所以我们可以把\(f(p_j)\)提出来。结果就是下面这一个大柿子：
然后考虑第一维怎么搞，好像只能预处理那就预处理一下。先看前边半拉\(g(\lfloor \frac n{p_j} \rfloor,j-1)\)。显然这个一定会落到\(\lfloor \frac nx \rfloor\)的范围内，\(x\)是随便一个数。而这样的数根据整除分块只有\(\sqrt n\)个，也就是说我们所有要预处理的状态是\(\sqrt n\)级别的。然后看后面的\(p_{j-1}\)，由于我们每个合数的最小质因子绝对不超过\(\sqrt n\)，所以我们在处理前面的\(\sqrt n\)个状态的时候已经把所有\(\sqrt n\)以内的\(p\)处理过了（具体证明可以看洛谷日报）。所以说我们预处理所有\(\sqrt n\)以内的\(g(i,0)\)就行了。放这题就是一个\(1-n\)的整数和还有一个\(1-n\)的整数平方和。就没有任何限制吗。然后\(dp\)暴力就行，防止MLE压掉一维。还有带着个既不是质数又不是合数的\(1\)有点麻烦，去掉。还有看不懂的可以看注释。
long long find(long long x){
return x<sq?pos1[x]:pos2[n/x];//就是把离散化后的值找回来
}
scanf("%lld",&n);
sq=sqrt(n);
get(sq);//筛素数
for(long long l=1,r;l<=n;l=r+1){
r=n/(n/l);//整除分块处理出根号个状态
w[++m]=n/l;
if(w[m]<sq)pos1[w[m]]=m;
else pos2[n/w[m]]=m;
//这里显然你要开n的数组是存不下的 但是观察到根号前和根号后的取值均不超过根号n所以考虑离散化到根号范围内的数组
g1[m]=(1ll*w[m]*(w[m]+1)%mod*inv2%mod-1+mod)%mod;//1-w[m]的和 即使第二维是0
g2[m]=(1ll*w[m]*(w[m]+1)%mod*(2ll*w[m]+1)%mod*inv6%mod-1+mod)%mod;//1-w[m]的平方和
}
for(int j=1;j<=cnt;j++){//dp第二维
for(int i=1;i<=m&&1ll*p[j]*p[j]<=w[i];i++){
g1[i]=(g1[i]-1ll*p[j]*(g1[find(w[i]/p[j])]-g1[find(p[j-1])])%mod+mod)%mod;//按照定义dp即可
g2[i]=(g2[i]-1ll*p[j]*p[j]%mod*(g2[find(w[i]/p[j])]-g2[find(p[j-1])])%mod+mod)%mod;
}
}
第二部分
之后直接把\(g(n,x)\)（\(x\)是小于\(\sqrt n\)的最大的素数）写成\(g(n)\)。
既然已经dp了那就dp到底吧。设
注意这里的\(f(i)\)是原函数，不是什么完全积性函数。也就是说我们不需要什么完全积性的性质。于是我们最后的答案就是\(s(n,0)\)。
仍然考虑把质数和合数的贡献分开。Min_25筛有另一个名字叫“扩展埃氏筛”，就是因为这一步里要像埃氏筛一样枚举质数的倍数转移。质数的贡献就是所有大于\(p_j\)的质数，就是\(g(n)-g(p_j)\)。然后考虑合数。我们枚举剩下的所有质数的倍数来转移前缀和。具体地就像这样：
来解释一下右面一堆东西是什么牛马。类似埃氏筛，我们用\(k\)枚举所有比\(j\)大的质数，用\(e\)枚举质数的次幂。然后后面的\(s\)是这个要从所有\(p_k^e\)的倍数里最小质因数比\(k\)大的转移状态（显然这可以枚举全部状态）。最后的\(e\neq 1\)是去掉质数\(p_k\)的贡献（因为前面已经算过了）。
这个甚至不需要记忆化，直接暴力递归就可以（根据joke3579对板子题的记忆化\(10^{11}\)的实验，扫到记忆化的次数是0）。没看懂的仍然看注释。
long long s(long long x,long long y){
if(p[y]>=x)return 0;
long long ret=(g2[find(x)]-g1[find(x)]-g2[find(p[y])]+g1[find(p[y])]+mod+mod)%mod;
//看着多 其实f=g2-g1显然 这个就是把括号拆开
for(int i=y+1;i<=cnt&&1ll*p[i]*p[i]<=x;i++){//枚举第y项后的所有质数
for(long long j=p[i];j*p[i]<=x;j=j*p[i]){//枚举质数的倍数
ret=(ret+f(j)*s(x/j,i)%mod+f(j*p[i])%mod)%mod;//直接暴力递归转移
//这个f(j*p[i])就是+1来跳过次数是1的部分 所以要使次数+1之后仍然小于x
//j=p[i]时这个是p[i]^2 依此类推
}
}
return ret;
}
printf("%lld",(s(n,0)+1)%mod);//最后别忘了加上没算的f(1)
顺便地，补一下复杂度证明。两步是类似的，都采用了类似的枚举质数倍数的方式（实际上Min_25筛还有一种非递归写法，第二步和第一步枚举质数的方法差不多）。于是复杂度就可以变成：
解释一下高等数学相关吧（虽然这是Min_25筛的学习笔记）：（参考了这个）
第一个等号就是根据定义列柿子。第二个也差不多，就是根据素数定理拆开。然后后面的显然大于前面的，舍掉前面然后拆成积分。\(\ln\)和\(\log\)差不多大所以整个常用一点的。然后关于这个积分怎么积：
首先我们看到下面真没法积，所以考虑一个近似做法，把下面变成\(\log n\)。反正带着\(\log\)怎么变也无所谓，常数区别。然后把它带回去看上面：
然后空间复杂度显然\(O(\sqrt n)\)。
于是洛谷的板子也就可以整出来了。贴个完整代码。
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
const int mod=1e9+7,inv2=500000004,inv6=166666668;
long long n,m,sq,cnt,g1[200010],g2[200010],w[200010];
int p[200010],pos1[200010],pos2[200010];
bool v[200010];
long long s1(long long x){
x%=mod;
return x*(x+1)%mod*inv2%mod;
}
long long s2(long long x){
x%=mod;
return x*(x+1)%mod*(2*x+1)%mod*inv6%mod;
}
long long Sq(long long x){
x%=mod;
return x*x%mod;
}
long long f(long long x){
x%=mod;
return (Sq(x)-x+mod)%mod;
}
long long find(long long x){
return x<sq?pos1[x]:pos2[n/x];
}
void get(int n){
for(int i=2;i<=n;i++){
if(!v[i])p[++cnt]=i;
for(int j=1;j<=cnt&&i*p[j]<=n;j++){
v[i*p[j]]=true;
if(i%p[j]==0)break;
}
}
}
long long s(long long x,long long y){
if(p[y]>=x)return 0;
long long ret=(g2[find(x)]-g1[find(x)]-g2[find(p[y])]+g1[find(p[y])]+mod+mod)%mod;
for(int i=y+1;i<=cnt&&1ll*p[i]*p[i]<=x;i++){
for(long long j=p[i];j*p[i]<=x;j=j*p[i]){
ret=(ret+f(j)*s(x/j,i)%mod+f(j*p[i])%mod)%mod;
}
}
return ret;
}
int main(){
scanf("%lld",&n);
sq=sqrt(n);
get(sq);
for(long long l=1,r;l<=n;l=r+1){
r=n/(n/l);
w[++m]=n/l;
if(w[m]<sq)pos1[w[m]]=m;
else pos2[n/w[m]]=m;
g1[m]=(s1(w[m])-1+mod)%mod;
g2[m]=(s2(w[m])-1+mod)%mod;
}
for(int j=1;j<=cnt;j++){
for(int i=1;i<=m&&1ll*p[j]*p[j]<=w[i];i++){
g1[i]=(g1[i]-1ll*p[j]*(g1[find(w[i]/p[j])]-g1[find(p[j-1])])%mod+mod)%mod;
g2[i]=(g2[i]-Sq(p[j])*(g2[find(w[i]/p[j])]-g2[find(p[j-1])])%mod+mod)%mod;
}
}
printf("%lld",(s(n,0)+1)%mod);
return 0;
}
（vscode编辑markdown这么好看的吗没搞过的一定要试试虽然开着预览打字的时候会一跳一跳的）
第二步有一个类似的方法，把递归改成递推：
意义显然。最后减一个 \(g\left(\min\left(\left\lfloor\dfrac n{p_{j+1}^e}\right\rfloor,p_{j+1}\right)\right)\) 是由于 \(s\left(\left\lfloor\dfrac n{p_{j+1}^e}\right\rfloor,j+1\right)\) 带上了 \(\le p_{j+1}\) 的质数的贡献，应当减掉。但是这个 \(\min\) 看着很难受，所以改一下写法：
这个跑的慢点，但是可以计算出所有的 \(s\left(\left\lfloor\dfrac nx\right\rfloor\right)\) 。同样上个代码。
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#define int long long
using namespace std;
const int mod=1e9+7,inv2=500000004,inv6=166666668;
long long n,m,sq,cnt,g1[200010],g2[200010],w[200010],s[200010];
int p[200010],pos1[200010],pos2[200010];
bool v[200010];
long long s1(long long x){
x%=mod;
return x*(x+1)%mod*inv2%mod;
}
long long s2(long long x){
x%=mod;
return x*(x+1)%mod*(2*x+1)%mod*inv6%mod;
}
long long Sq(long long x){
x%=mod;
return x*x%mod;
}
long long f(long long x){
x%=mod;
return (Sq(x)-x+mod)%mod;
}
long long find(long long x){
return x<sq?pos1[x]:pos2[n/x];
}
void get(int n){
for(int i=2;i<=n;i++){
if(!v[i])p[++cnt]=i;
for(int j=1;j<=cnt&&i*p[j]<=n;j++){
v[i*p[j]]=true;
if(i%p[j]==0)break;
}
}
}
signed main(){
scanf("%lld",&n);
sq=sqrt(n);
get(sq);
for(long long l=1,r;l<=n;l=r+1){
r=n/(n/l);
w[++m]=n/l;
if(w[m]<sq)pos1[w[m]]=m;
else pos2[n/w[m]]=m;
g1[m]=(s1(w[m])-1+mod)%mod;
g2[m]=(s2(w[m])-1+mod)%mod;
}
for(int j=1;j<=cnt;j++){
for(int i=1;i<=m&&p[j]*p[j]<=w[i];i++){
g1[i]=(g1[i]-p[j]*(g1[find(w[i]/p[j])]-g1[find(p[j-1])])%mod+mod)%mod;
g2[i]=(g2[i]-Sq(p[j])*(g2[find(w[i]/p[j])]-g2[find(p[j-1])])%mod+mod)%mod;
}
}
//以下
for(int i=1;i<=m;i++)s[i]=g1[i]=(g2[i]-g1[i]+mod)%mod;
for(int j=cnt;j>=1;j--){
for(int i=1;i<=m&&p[j]*p[j]<=w[i];i++){
for(int tmp=p[j],e=1;p[j]*tmp<=w[i];e++,tmp*=p[j]){
s[i]=(s[i]+f(tmp)*(s[find(w[i]/tmp)]-g1[find(p[j])]+mod)%mod+f(tmp*p[j]))%mod;
}
}
}
printf("%lld\n",(s[1]+1)%mod);
return 0;
}
例题
- 求区间素数个数。
首先我们随手差分一下变成\(1-n\)的素数个数。然后一眼就是Min_25筛的第一步。就是求个\(f(i)=1\)在质数取值的和就可以了。显然这个东西完全积性。代码就不上了，反正贼水。
题面：有个积性函数\(f\)，满足：\(f(1)=1,f(p^c)=p\oplus c\)（\(\oplus\)是异或）。求它的前缀和，数据范围\(10^{10}\)。
这个也还可以。观察质数取值有：\(f(2)=2+1,f(p)=p-1\)。于是我们把它们全看成\(f(p)=p-1\)，然后直接把它当做这个筛，最后特判\(2\)（就是结果加个\(2\)）。对于这个函数，我们可以把它拆成两个完全积性函数\(f(i)=i\)和\(f(i)=1\)(注意不能是\(-1\),因为\((-1)\times (-1)=1\))然后第一步乱搞就可以。第二步直接套板子。
上个代码。其实就是大板子。
long long s(long long x,long long y){
if(p[y]>=x)return 0;
long long ret=g1[find(x)]-g2[find(x)]-g1[find(p[y])]+g2[find(p[y])];
for(int i=y+1;i<=cnt&&1ll*p[i]*p[i]<=x;i++){
for(long long e=1,j=p[i];j*p[i]<=x;e++,j*=p[i]){
ret+=(p[i]^e)*s(x/j,i)+(p[i]^(e+1));
}
}
if(y==0)ret+=2;
return ret;
}
int main(){
scanf("%lld",&n);
sq=sqrt(n);
get(sq);
for(long long l=1,r;l<=n;l=r+1){
r=n/(n/l);
w[++m]=n/l;
if(w[m]<=sq)pos1[w[m]]=m;
else pos2[n/w[m]]=m;
g1[m]=w[m]*(w[m]+1)/2-1;
g2[m]=w[m]-1;
}
for(int j=1;j<=cnt;j++){
for(int i=1;i<=m&&p[j]*p[j]<=w[i];i++){
g1[i]=g1[i]-p[j]*(g1[find(w[i]/p[j])]-g1[find(p[j-1])]);
g2[i]=g2[i]-(g2[find(w[i]/p[j])]-g2[find(p[j-1])]);
}
}
printf("%lld",s(n,0)+1);
return 0;
}
- 求莫比乌斯函数\(\mu\)和欧拉函数\(\varphi\)的前缀和。
显然这个可以杜教筛（废话这是杜教筛板子）。但是我们现在讲的是Min_25筛。
先看\(\mu\)。看它的质数取值全是\(-1\)，那第一步就解决了。第二步发现对于\(\forall k>1,\mu(p^k)=0\)。于是我们甚至连质数次幂都不需要枚举，反正都是\(0\)。
然后看\(\varphi\)。我们有\(\varphi(p)=p-1,\varphi(p^k)=p^{k-1}(p-1)\)。这个就直接套板子筛就行了。代码不上了，懒了。
题意：求
数据范围\(10^{10}\)。
首先我们看三个\(\sum\)有那么亿点点难受，消掉一个。把\(j+k\)看做\(i\)的一个因数\(l\)，那么答案就是：
由于我们同样选两个因数，一次\(j\)更大，一次\(l\)更大，算了两次，所以要除以\(2\)。然后根据组合意义，这就是从\(i\)的因数里随便选两个出来有多少选法，也就是：
其中\(d(i)\)表示 \(i\) 有多少因数。然后随手把组合数拆开变成：
Min_25筛板子。
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
const int mod=998244353;
long long n,m,sq,cnt,sum,g[200010],w[200010];
int p[200010],pos1[200010],pos2[200010];
bool v[200010];
long long find(long long x){
return x<sq?pos1[x]:pos2[n/x];
}
void get(int n){
for(int i=2;i<=n;i++){
if(!v[i])p[++cnt]=i;
for(int j=1;j<=cnt&&i*p[j]<=n;j++){
v[i*p[j]]=true;
if(i%p[j]==0)break;
}
}
}
long long f(long long x){
return 1ll*(x+1)*(x+1)%mod;
}
long long s(long long x,long long y){
if(p[y]>=x)return 0;
long long ret=4ll*(g[find(x)]-g[find(p[y])]+mod)%mod;
for(int i=y+1;i<=cnt&&1ll*p[i]*p[i]<=x;i++){
for(long long j=p[i],e=1;j<=x;j=j*p[i],e++){
ret=(ret+1ll*(e+1)*(e+1)*(s(x/j,i)+(e>1)))%mod;
}
}
return ret;
}
int main(){
scanf("%lld",&n);
sq=sqrt(n);
get(sq);
for(long long l=1,r;l<=n;l=r+1){
r=n/(n/l);
w[++m]=n/l;
sum=(sum+1ll*w[m]*(r-l+1))%mod;
if(w[m]<sq)pos1[w[m]]=m;
else pos2[n/w[m]]=m;
g[m]=(w[m]-1)%mod;
}
for(int j=1;j<=cnt;j++){
for(int i=1;i<=m&&1ll*p[j]*p[j]<=w[i];i++){
g[i]=(g[i]-1ll*(g[find(w[i]/p[j])]-g[find(p[j-1])])%mod+mod)%mod;
}
}
long long ans=1ll*(s(n,0)+1-sum+mod)%mod*((mod+1)>>1)%mod;
printf("%lld\n",ans);
return 0;
}
- 来一个压轴题：开头的LOJ572 Misaka Network 与求和
题意：求
其中\(f\)表示次大质因数，重复的质因数算多次，比如\(f(4)=2\)。规定\(f(p)=1,f(1)=0\)。数据范围\(n\le 2\times 10^9\)。
首先随手套路反演一波（反演再不会的看我之前博客）（都写这么多行了真懒得一步一步化了）：
然后看后面那一堆东西，我们要求\(\mu * f^k\)的前缀和。看见卷积形式，考虑杜教筛化简一波（好怪啊我杜教筛怎么都忘了）。看到卷的是\(\mu\)，我们卷一个\(I\)上去。于是由\(\mu * I=\epsilon\)，我们把这个\(\mu\)去掉了，现在的问题就是求\(f^k\)的前缀和。确实不好搞，而且\(2\times 10^{9}\)更难搞了。
但考虑我们Min_25筛的过程，第二步中我们在枚举每一个可能成为最小质因子的质数并计算它的贡献。我们换一个角度思考，这一个步骤可以是枚举次大质因子的步骤。于是我们可以改变第二步\(s(n,j)\)的定义，是所有\(n\)以内的次大质因子\(\ge p_j\)的\(f^k\)之和。
考虑如何对状态进行转移。每次我们扫描一个次大质因子的时候，我们可以从次大质因子比它大的数转移，还要加上当前次大质因子的贡献，还有所有之后的质数。一步一步来。当前次大质因子的贡献就可以直接找有多少质数比它大来计算，就是\(p_j^k(\sum[i\in P\ and\ i\le p_j])\)。解释一下就是在该数为最大质因子的数之后乘一个更大的质因子，该数就变成了次大质因子，直接计算更大的质数的贡献。然后后边的合数可以类似地递归，整个的就是
具体解释一下。枚举质因子\(p_i\)时，它有贡献则这个数除\(p_i\)为\(1\)或质数。如果不是的话可以除掉所有\(p_i\)然后递归继续枚举下一个质因子。之后就长的和Min_25筛板子差不多了。复杂度因为杜教筛没法预处理所以是\(O(n^{\frac 34})\)。请大力卡常。代码没有因为不会（
终于结束了。
