$Min\_25$篩學習筆記

時間 2019-11-09

標籤 min 學習筆記简体版

原文原文鏈接

$Min\_25$篩學習筆記

這種神仙東西不寫點東西一下就忘了QAQ

資料和代碼出處
 資料2
資料3
打死我也不認可參考了yyb的html

$Min\_25$篩能夠幹嗎？下文中未特殊說明$P$均指質數集合,$p_i$或$p$指某個具體質數。函數

求一類積性函數$f(x)$的前綴和，須要知足$f(p)$能夠寫成多項式的形式，或者操做一下能夠寫成多項式（如例題），且$f(p^k)$能快速求出。學習

講真學這個東西比我什麼都不會的時候學$FFT$都累。spa

Round 1

先求質數的貢獻。咱們要求code

\[\sum_{i=1}^x[i是質數]f(i)\]。htm

咱們能夠設blog

\[g(n,j)=\sum_{i=1}^{n}[i \in P \ or\ \min(p)>P_j]f(i)\]get

注意到當$p_j^2>x$時，$g(n,j)=g(n,j-1)$由於沒有新的貢獻了，因此咱們只要篩到$\sqrt n$的質數就能夠了，這也是下降複雜度的關鍵。string

咱們還注意到$g(n,|P|)=\sum_{i=1}^x[i是質數]f(i)$，$|P|$是$\sqrt n$內的質數集合大小。it

關於$g(n,j)$的轉移，咱們有：

\[g(n,j)=g(n,j-1)-f(P_j)[g(\frac{n}{P_j},j-1)-\sum_{i=1}^{j-1}f(P_i)],P_j^2\le n\]

意思大概就是減掉全部最小質因子爲$p_j$的貢獻，但因爲$g(\frac n{p_j},j-1)$裏包含了質數，而$<p_j$的質數是不能算的，因此要減掉。

因爲咱們只須要$g(n,|P|)$，因此下面的代碼是一維的，用遞推實現。注意到過程當中咱們只須要$g(\lfloor\frac ni\rfloor,|P|)$，因此最多隻有$2\sqrt n$種取值

至於實現，因爲我參考了$gsy$的實現，痛不欲生，因而決定把他的代碼蒯走。這份代碼是篩$f(p)=1$的。

//這兩個鬼id就是你在杜教篩中碰到的卡常卡空間技巧，這份代碼你理解了這個就能看懂
//至於YL，是機房衆人mo的巨佬，因此是模數
//Sq是根號n
for (int i=1,j;i<=n;i=j+1)
{
    j=n/(n/i);w[++m]=n/i;g[m]=(w[m]-1)%YL;//除法分塊，根號n求出全部有用的值
    if(w[m]<=sq)id1[w[m]]=m;else id2[n/w[m]]=m;
}
for (int j=1;j<=tot;++j)
    for (int i=1;i<=m&&pri[j]*pri[j]<=w[i];++i)//i再往上就是全部的質數，會被後面抵消
    {
        int k=(w[i]/pri[j]<=sq)?id1[w[i]/pri[j]]:id2[n/(w[i]/pri[j])];
        g[i]=(g[i]-g[k]+(j-1))%YL;g[i]=(g[i]+YL)%YL;
    }

Round 2

然而敵人並無這麼容易就被打倒，咱們還有合數沒算呢。那麼咱們鼓搗一個$S(n,i)$出來

\[S(n,j)=\sum_{i=1}^n[\min(p)\ge P_j]f(i)\]

注意$S(n,1)$沒算到$f(1)$。答案就是$S(n,1)+f(1)$。

遞推式來了

\[S(n,j)=(質數的貢獻)+(合數的貢獻)\]

有沒有感受被騙了$QwQ$，咱們繼續

\[(質數的貢獻)=g(n,|P|)-\sum_{i=1}^{j-1}f(P_i)\]

\[(合數的貢獻)=\sum_{k=j}^{P_k^2\le n}\sum_{e=1}^{P_k^{e+1}\le n}S(\frac{n}{P_k^e},k+1)\times f(P_k^e)+f(P_k^{e+1})\]

通過$YL$的指點，我能夠口胡一下了，每一個合數要在最小質因子處被篩到。當$\frac n{p_k^e}<p_{k+1}$時，確定是$0$就不必繼續了。

舉個栗子，形如$t*p_k^3$（其中$t$的最小質因子大於$p_k$）的數會被$S(\frac n{p_k^3},4)$包含。

因爲$S$這個函數在任什麼時候候都不包含$f(1)$因此咱們要手動加上$f(p_k^2),f(p_k^3),...,f(p_k^{e+1})$其中$p_k^{e+2}>n$。什麼，你問我$f(p_k)$去哪了，這不是個質數麼。

完結撒花

放上LOJ6053簡單的函數做爲例題，這裏面的$p$^$c$對於質數來講除了$2$之外都是$p-1$，像個多項式。

代碼中$g(x,|P|)=\sum_{i=1}^x[i是質數]i$，$h(x,|P|)=\sum_{i=1}^x[i是質數]1$

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define gt getchar()
#define ll long long
#define File(s) freopen(s".in","r",stdin),freopen(s".out","w",stdout)
inline ll in()
{
    ll k=0;char ch=gt;
    while(ch<'-')ch=gt;
    while(ch>'-')k=k*10+ch-'0',ch=gt;
    return k;
}
const int N=1e6+5,YL=1e9+7,mod=YL;
inline int MO(const int &x){return x>=YL?x-YL:x;}
ll n,w[N];
int np[N],pr[N],tot,h[N],g[N],sp[N],id1[N],id2[N],m,sq;
void seive(int n)
{
    for(int i=2;i<=n;++i)
    {
        if(!np[i])pr[++tot]=i,sp[tot]=MO(sp[tot-1]+i);
        for(int j=1;i*pr[j]<=n;++j)
        {np[i*pr[j]]=1;if(i%pr[j]==0)break;}
    }
}

int S(ll x,int y)
{
    if(x<=1||pr[y]>x)return 0;
    int k=(x<=sq?id1[x]:id2[n/x]);
    int res=MO(((ll)g[k]-h[k]-sp[y-1]+y-1)%YL+YL);
    if(y==1)res+=2;
    for(int i=y;i<=tot&&1ll*pr[i]*pr[i]<=x;++i)
    {
        ll p1=pr[i],p2=p1*pr[i];
        for(int e=1;p2<=x;++e,p1=p2,p2*=pr[i])
            res=MO(res+(1ll*S(x/p1,i+1)*(pr[i]^e)+(pr[i]^e+1))%YL);
    }
    return res;
}

int main()
{
    n=in();sq=sqrt(n);seive(sq);ll t;
    for(ll i=1,j;i<=n;i=j+1)
    {
        j=n/(n/i),w[++m]=n/i;
        if(w[m]<=sq)id1[w[m]]=m;else id2[n/w[m]]=m;
        h[m]=(w[m]-1)%YL;g[m]=((w[m]+2)%YL)*((w[m]-1)%YL)%YL;
        if(g[m]&1)g[m]+=YL;g[m]>>=1;
    }
    
    for(int j=1;j<=tot;++j)
        for(int i=1;i<=m&&1ll*pr[j]*pr[j]<=w[i];++i)
        {
            t=w[i]/pr[j];int k=(t<=sq?id1[t]:id2[n/t]);
            h[i]=MO((h[i]-h[k]+j-1)%YL+YL);
            g[i]=MO(MO(g[i]-1ll*pr[j]*(g[k]-sp[j-1])%YL)+YL);
        }
    printf("%d\n",S(n,1)+1);
    return 0;
}