1. 程式人生 > >Min_25篩

Min_25篩

最小 long lin 並不是 isp clas pan 素數個數 導致

聽說這個東西能給予人力量

那就來學一學吧

功能就是篩一個積性函數\(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)\)

表示\(i\)的質因子,\(P_j\)表示第\(j\)個質數

看起來可能沒什麽含義,但是我們可以理解為正在進行一個埃篩,\(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\)

,最小質因子自然就是\(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篩