min_25筛学习笔记
用途
筛一个积性函数\(f(n)\)的前缀和\(S(n)=\sum_{i=1}^n f(n)\)，这个函数还需要满足的性质是\(f(p)\)是一个关于\(p\)的多项式、\(f(p^k)\)可以快速计算。
思想
min_25筛的核心思想就是每个合数的最小质因子不超过\(Sqr=\sqrt{n}\)，所以考虑质数单独提出来算，合数在它最小的质数处统计。
第一步
在这一步当中，我们要求出
也就是所有质数的前缀和。
考虑积性函数的性质还不够优秀，不如把多项式的各项拆开来分别算，最后再合在一起。
于是在这一步中假定\(f(p)=p^t\)，它就完全积性了。
记\(P_j\)为从小到大第\(j\)个质数，设
显然，要求所有质数的\(f\)值和，就是\(g(n,|P|)\)。
而\(g\)的初值比较麻烦，但考虑所有不为质数的值最终都会被筛掉，那么一开始可以把所有\(f(i)\)当成质数来算，即\(f(i)=i^t\)。
从小到大枚举\(j\)来筛掉所有的合数（思想类似于埃氏筛），有如下转移式：
第一个式子：\(P_j^2>n\)了，所以不会有多的数被筛掉。
第二个式子：把最小质因子是\(P_j\)的筛掉。
（注意现在这个函数的性质非常好，所以当成质数来算还是可以有完全积性函数的性质，可以直接乘在一起）
于是可以\(O(n\times |P|)\)搞出\(n\)个数的\(g(n,|P|)=\sum_{i=1}^n [i\in Prime]f(i)\)，但这样还是太多了。在第二步中我们会发现其实根本不需要求出所有\(g(n,|P|)\)。
第二步
在这一步中，我们要利用积性函数的性质，把所有合数的\(f\)值加回来。
请注意在这一步中\(f(n)\)不再完全积性了。
记
那么可以先把质数的答案加上，然后枚举最小的质因子，又会有递推式：
（这里面\(f(P_k^e)S(\frac{n}{P_k^e},k+1)\)可以直接乘起来的原因是\(S(n/P_k^e,k+1)\)中所有元素都与\(P_k\)互质）
注意\(S(n,j)\)只调用了\(g(n,j)\)和\(S(\frac{n}{P_k^e},k+1)\)，所以第一步中只需要算出所有\(g(n/i,j)\)就够了，而这是\(O(\sqrt{n}\times |P|)\)的。
这东西直接算就好了，不用记忆化。
最后总复杂度可以被证明是\(O(\frac{n^{3/4}}{\log n})\)？反正我不会证明，直接用就好了。
代码
这里以洛谷P5325 【模板】Min_25筛为例：
#include<bits/stdc++.h>
clock_t t=clock();
namespace my_std{
using namespace std;
#define pii pair<int,int>
#define fir first
#define sec second
#define MP make_pair
#define rep(i,x,y) for (int i=(x);i<=(y);i++)
#define drep(i,x,y) for (int i=(x);i>=(y);i--)
#define go(x) for (int i=head[x];i;i=edge[i].nxt)
#define templ template<typename T>
#define sz 2020200
#define mod 1000000007ll
typedef long long ll;
typedef double db;
mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());
templ inline T rnd(T l,T r) {return uniform_int_distribution<T>(l,r)(rng);}
templ inline bool chkmax(T &x,T y){return x<y?x=y,1:0;}
templ inline bool chkmin(T &x,T y){return x>y?x=y,1:0;}
templ inline void read(T& t)
{
t=0;char f=0,ch=getchar();double d=0.1;
while(ch>'9'||ch<'0') f|=(ch=='-'),ch=getchar();
while(ch<='9'&&ch>='0') t=t*10+ch-48,ch=getchar();
if(ch=='.'){ch=getchar();while(ch<='9'&&ch>='0') t+=d*(ch^48),d*=0.1,ch=getchar();}
t=(f?-t:t);
}
template<typename T,typename... Args>inline void read(T& t,Args&... args){read(t); read(args...);}
char __sr[1<<21],__z[20];int __C=-1,__zz=0;
inline void Ot(){fwrite(__sr,1,__C+1,stdout),__C=-1;}
inline void print(register int x)
{
if(__C>1<<20)Ot();if(x<0)__sr[++__C]='-',x=-x;
while(__z[++__zz]=x%10+48,x/=10);
while(__sr[++__C]=__z[__zz],--__zz);__sr[++__C]='\n';
}
void file()
{
#ifdef NTFOrz
freopen("a.in","r",stdin);
#endif
}
inline void chktime()
{
#ifndef ONLINE_JUDGE
cout<<(clock()-t)/1000.0<<'\n';
#endif
}
#ifdef mod
ll ksm(ll x,int y){ll ret=1;for (;y;y>>=1,x=x*x%mod) if (y&1) ret=ret*x%mod;return ret;}
ll inv(ll x){return ksm(x,mod-2);}
#else
ll ksm(ll x,int y){ll ret=1;for (;y;y>>=1,x=x*x) if (y&1) ret=ret*x;return ret;}
#endif
// inline ll mul(ll a,ll b){ll d=(ll)(a*(double)b/mod+0.5);ll ret=a*b-d*mod;if (ret<0) ret+=mod;return ret;}
templ inline T MOD(T x) {return x - ((mod-x-1) >> 31 & mod);}
}
using namespace my_std;
const ll inv6=inv(6);
ll n,Sqr;
inline ll f(register ll x){x%=mod;return x*(x-1)%mod;}
int pri[sz],cnt;
ll sum1[sz],sum2[sz];
bool npri[sz];
void init()
{
#define x (i*pri[j])
rep(i,2,sz-1)
{
if (!npri[i])
pri[++cnt]=i,
sum1[cnt]=(sum1[cnt-1]+i)%mod,
sum2[cnt]=(sum2[cnt-1]+1ll*i*i%mod)%mod;
for (int j=1;j<=cnt&&x<sz;j++)
{
npri[x]=1;
if (i%pri[j]==0) break;
}
}
#undef x
}
int id1[sz],id2[sz],m;
ll w[sz];
ll g1[sz],g2[sz];
ll S(const ll &N,const int &j)
{
if (N<=1) return 0;
int id=(N>=Sqr?id2[n/N]:id1[N]);
ll ret=(g2[id]-g1[id]-(sum2[j-1]-sum1[j-1])+mod+mod)%mod;
for (int k=j;1ll*pri[k]*pri[k]<=N;++k)
for (ll P=pri[k];P*pri[k]<=N;P=P*pri[k])
ret=MOD(ret+MOD(f(P)*S(N/P,k+1)%mod+f(P*pri[k])));
return ret;
}
int main()
{
file();
read(n);Sqr=sqrt(n);
init();
for (ll i=1,j;i<=n;i=j+1)
{
ll x=w[++m]=n/i;
j=n/x;
if (x<Sqr) id1[x]=m;
else id2[j]=m;
x%=mod;
g1[m]=(x*(x+1)/2)%mod-1;
g2[m]=x*(x+1)%mod*(x+x+1)%mod*inv6%mod-1;
}
for (int i=1;pri[i]<=Sqr;i++) rep(N,1,m)
{
if (1ll*pri[i]*pri[i]>w[N]) break;
ll x=w[N]/pri[i];
int id=(x>=Sqr?id2[n/x]:id1[x]);
g1[N]=MOD(g1[N]-1ll*pri[i]*MOD(g1[id]-sum1[i-1]+mod)%mod+mod);
g2[N]=MOD(g2[N]-1ll*pri[i]*pri[i]%mod*MOD(g2[id]-sum2[i-1]+mod)%mod+mod);
}
cout<<(S(n,1)+1)%mod<<'\n';
return 0;
}
更新
2021.7.25 来更新一下这篇远古 blog ，权当攒 rp 。
首先补一个复杂度证明：在第一步中，设 \(m=\lfloor {n/i}\rfloor\) ，则 \(m\) 要用到的质数只有 \(O({\sqrt m\over \log m})\) 个，所以所有 \(m\) 加起来就大约是
第二步中用了递归且没有记忆化，而且还枚举了 \(m\) 的最小质数是哪个，复杂度就看起来比较迷惑了。
更新当然不能只补一个复杂度证明，所以还要介绍从麦老大那边学过来的一个更好写的做法，并且能把所有 \(n/i\) 的答案一起算出来。
第一步的做法维持不变，代码也不变。
第二步中，定义 \(S(m,k)=\sum_{i=1}^m [i\in \text{prime}\lor \text{minp}(i)\ge k] f(i)\) 。初值为 \(S(m,+\infty)=g(m)\) 。
从大到小枚举 \(k\) ，并枚举 \(m\ge p_k^2\) ，把 \(\text{minp}(i)=k\) 的 \(i\) 给加上。这只需要枚举 \(p_k\) 的次数，从 \(S(m/p_k^e,k+1)\) 转移过来即可。
除了枚举 \(p_k\) 的次数以外复杂度就和第一步完全一样了。
以下是洛谷模板的代码：好像也没短多少
int pri[sz],npri[sz],cc;
ll S1[sz],S2[sz];
void init(){rep(i,2,sz-1) if (!npri[i]) { pri[++cc]=i; S1[cc]=(S1[cc-1]+i)%mod,S2[cc]=(S2[cc-1]+1ll*i*i)%mod; for (int j=i+i;j<sz;j+=i) npri[j]=1; } }
ll n,B;
int id1[sz],id2[sz],m; ll w[sz];
int& id(ll x){return x<=B?id1[x]:id2[n/x];}
ll g1[sz],g2[sz];
ll S[sz];
int main()
{
file();
init();
read(n); B=sqrt(n);
for (ll l=1,r;l<=n;l=r+1)
{
r=n/(n/l); ll x=n/l;
w[++m]=x; id(x)=m;
x%=mod;
g1[m]=x*(x+1)/2%mod-1; g2[m]=x*(x+1)%mod*(x+x+1)%mod*inv(6)%mod-1;
}
rep(i,1,cc) rep(j,1,m)
if (1ll*pri[i]*pri[i]<=w[j])
(g1[j]+=mod-pri[i]*(g1[id(w[j]/pri[i])]-S1[i-1])%mod)%=mod,
(g2[j]+=mod-1ll*pri[i]*pri[i]%mod*(g2[id(w[j]/pri[i])]-S2[i-1])%mod)%=mod;
else break;
rep(i,1,m) S[i]=(g2[i]-g1[i]+mod)%mod;
drep(i,cc,1) rep(j,1,m)
if (1ll*pri[i]*pri[i]<=w[j])
for (ll k=pri[i];k*pri[i]<=w[j];k*=pri[i],S[j]+=(k%mod)*(k%mod-1)%mod)
(S[j]+=(k%mod)*(k%mod-1)%mod*(S[id(w[j]/k)]-S2[i]+S1[i]+mod)%mod)%=mod;
else break;
chktime();
cout<<(S[1]+1)%mod;
return 0;
}
