Min_25篩
聽說這個東西能給予人力量
那就來學一學吧
功能就是篩一個積性函數\(f(i)\)的前綴和
Min_25篩好像是最近才流行起來的篩法,復雜度是非常神奇的\(O(\frac{n^{\frac{3}{4}}}{logn})\)
和杜教篩一樣,使用這個篩法的也有一定要求,就是\(f(p^c)\)需要在\(O(1)\)求出
來看看這個非常力量的篩法
我們要求的東西是
\[\sum_{i=1}^nf(i)\]
我們先定義一個二元函數\(g(n,j)\)
\[g(n,j)=\sum_{i=1}^n[i\in P\ or\ minp(i)>P_j]f(i)\]
\(P\)表示質數集,\(minp(i)\)
看起來可能沒什麽含義,但是我們可以理解為正在進行一個埃篩,\(P_j\)這個質數已經篩完了,這個時候所有被判定為質數的位置和所有沒被篩到的位置的\(f\)函數的和,就表示\(g(n,j)\)
考慮如何推\(g(n,j)\)
如果\(P_j^2>n\),那麽\(P_j\)這個數非常可憐,一個數都沒有被它篩到,因為\(P_j\)篩到的數最小也得是\(P_j^2\),所以就有
\[g(n,j)=g(n,j-1)\]
如果\(P_j^2<=n\),那麽我們就需要考慮\(P_j\)和\(g(n,j-1)\)篩到了什麽,之後減掉這些貢獻
篩掉的都是那些\(minp(i)=P_j\)的數吧,根據直覺
\[g(n,j)=g(n,j-1)-f(P_j)(g(\left \lfloor \frac{n}{P_j} \right \rfloor,j-1)-\sum_{i=1}^{j-1}f(P_i))\]
來解釋一下為什麽吧
首先那個\(\left \lfloor \frac{n}{P_j} \right \rfloor\)能保證這個範圍裏的數乘上\(P_j\)都不大於\(n\),之後這個範圍裏的數都被\(minp\)小於等於\(P_{j-1}\)的都被幹掉了,所以這個範圍內剩下的數的\(minp\)都大於等於\(P_j\),之後乘上\(P_j\)
但是我們這樣把質數也給算了上去,所以並不是很可行,所以後面還得把\(1\)到\(j-1\)的質數加回來
同時這裏有一個問題,我們將\(\left \lfloor \frac{n}{P_j} \right \rfloor\)範圍的數還原的時候直接乘上\(f(P_j)\)並不科學,因為這裏顯然會有不互質的情況出現,所以這裏是將\(f\)當做完全積性函數來做的
綜上
\[g(n,j)=\begin{cases} g(n,j-1)&P_j^2\gt n\\ g(n,j-1)-f(P_j)[g(\left \lfloor \frac{n}{P_j} \right \rfloor,j-1)-\sum_{i=1}^{j-1}f(P_i)]&P_j^2\le n\end{cases}\]
求出這個\(g\)什麽用都沒有,我們還需要再設
\[S(n,j)=\sum_{i=1}^n[i\in P\ or\ minp(i)>=P_j]f(i)\]
這裏就不在把\(f(i)\)當成完全積性了
這個時候我們的答案自然就是\(S(n,1)+1\)
我們把\(S(n,j)\)分成質數貢獻的一部分和合數貢獻的一部分
質數哪一部分的貢獻自然就是
\[g(n,|P|)-\sum_{i=1}^{j-1}f(P_j)\]
後面那個是減掉小於\(j\)的質數的貢獻
而合數的部分我們通過枚舉最小質因子和最小質因子出現的次數
\[\sum_{k=j}^{P_k^2<=n}\sum_{e=1}^{P_k^e<=n}f(P_k^e)(S(\left \lfloor \frac{n}{P_k^e} \right \rfloor,k+1)+[e>1])\]
這裏跟上面類似,由於後面是\(j+1\)所以會導致\(S(\left \lfloor \frac{n}{P_k^e} \right \rfloor,k+1)\)內算貢獻的數都最小質因子都大於等於\(P_{k+1}\),直接乘上\(P_k^e\)不可能產生不互質的情況
由於\(1\)的貢獻是最後算的,\(S(\left \lfloor \frac{n}{P_k^e} \right \rfloor,k+1)\)並沒有算到\(1\)也就是還原回來的\(P_k^e\)所以在最後需要加\(1\)
特殊的是\(e=1\)本身就是質數,之前算過貢獻這裏就不需要加\(1\)
代碼實現這個篩法就非常復雜了
來一道簡單的例題
loj#6235. 區間素數個數
顯然這個東西不是積性函數,但是考慮一下Min_25篩的第一步,就是只算質數部分的貢獻
我們這裏把\(f(i)=[i>1]\)就好了
代碼
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#define maxn 1000005
#define re register
#define LL long long
const LL mod=1000000007;
LL n,w[maxn],Sqr,tot;
LL p[maxn],pre[maxn],id1[maxn],id2[maxn],h[maxn],g[maxn],m;
int f[maxn];
int main()
{
scanf("%lld",&n);Sqr=std::sqrt(n);
f[1]=1,pre[1]=1;
for(re int i=2;i<=Sqr;i++) {
if(!f[i]) p[++tot]=i,pre[tot]=(pre[tot-1]+(i^1))%mod;
for(re int j=1;j<=tot&&p[j]*i<=Sqr;j++) {
f[p[j]*i]=1;if(i%p[j]==0) break;
}
}//預處理根號範圍內的質數
for(re LL l=1,r;l<=n;l=r+1) {
r=n/(n/l);w[++m]=n/l;
if(w[m]<=Sqr) id1[w[m]]=m;
else id2[n/w[m]]=m;
g[m]=w[m]-1;//這裏的g[m]=g(w[m],0),所以初值是除去1以外的數的個數
}
for(re int j=1;j<=tot;j++)
for(re int i=1;i<=m&&p[j]*p[j]<=w[i];i++) {
int k=(w[i]/p[j]<=Sqr)?id1[w[i]/p[j]]:id2[n/(w[i]/p[j])];
//求w[i]/p[j]在哪一個塊裏
g[i]=g[i]-g[k]+j-1;//這裏實際上就是從g(i,j-1)推向g(i,j)
}
printf("%lld\n",g[1]);
return 0;
}
loj#6053. 簡單的函數
首先註意到除去唯一的偶質數\(2\),所有的\(p\)都滿足\(p⊕1=p-1\)
而\(2⊕1=3\)
先設\(f(p)=p\),我們再求一個區間內素數個數求篩的時候減掉這個區間素數個數
如果有\(2\)這個質數記得把答案加\(2\)
之後也沒什麽了
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#define maxn 300005
#define re register
#define LL long long
const LL mod=1000000007;
const LL inv2=500000004;
LL n,w[maxn],Sqr,tot;
LL p[maxn],pre[maxn],id1[maxn],id2[maxn],h[maxn],g[maxn],m;
int f[maxn];
inline LL S(LL x,int y) {
if(x<=1||p[y]>x) return 0;
int k=(x<=Sqr)?id1[x]:id2[n/x];
LL ans=(h[k]-g[k]-pre[y-1]+y-1+((y==1)?2ll:0))%mod;ans=(ans+mod)%mod;
for(re int k=y;k<=tot&&p[k]*p[k]<=x;k++) {
LL p1=p[k];
for(re int e=1;p1<=x;e++,p1*=p[k])
ans+=(1ll*(S(x/p1,k+1)+((e>1)?1ll:0))*(p[k]^e)%mod),ans%=mod;
}
return ans;
}
int main()
{
scanf("%lld",&n);Sqr=std::sqrt(n)+1;
f[1]=1;
for(re int i=2;i<=Sqr;i++) {
if(!f[i]) p[++tot]=i,pre[tot]=(pre[tot-1]+i)%mod;
for(re int j=1;j<=tot&&p[j]*i<=Sqr;j++) {
f[p[j]*i]=1;if(i%p[j]==0) break;
}
}
for(re LL l=1,r;l<=n;l=r+1) {
r=n/(n/l);w[++m]=n/l;
if(w[m]<=Sqr) id1[w[m]]=m;
else id2[n/w[m]]=m;
g[m]=(w[m]-1)%mod;
h[m]=((w[m]+2ll)%mod)*((w[m]-1ll)%mod)%mod;h[m]=(h[m]*inv2)%mod;
}
for(re int j=1;j<=tot;j++)
for(re int i=1;i<=m&&p[j]*p[j]<=w[i];i++) {
int k=(w[i]/p[j]<=Sqr)?id1[w[i]/p[j]]:id2[n/(w[i]/p[j])];
g[i]=g[i]-g[k]+j-1ll; g[i]=(g[i]+mod)%mod;
h[i]=(h[i]-(h[k]-pre[j-1]+mod)%mod*p[j]%mod+mod)%mod;
}
printf("%lld\n",1+S(n,1));
return 0;
}
Min_25篩