1. 程式人生 > >Min_25篩學習筆記

Min_25篩學習筆記

感覺好好用啊

Luogu上的杜教篩模版題一發 Min_25搶到了 rank1

 


 

前言

$ Min$_$25$篩可以篩以下一類函式的字首和

$ f(1)=1$,$ f(p_i^{k_i}p_j^{k_j})=f(p_i^{k_i})f(p_j^{k_j})$

$ f(p_i)為一個多項式,f(p_i^{k_i})可以快速計算$

以下部分暫時忽略$ 1$,即只考慮最小質因子$ >=2$的那些數

 


 

先考慮素數貢獻

我們定義$ sp(n)$表示$\sum\limits_{i=1}^n f(p_i)$即前$ n$個素數的積性函式和

這裡我們先假設$ f$對於質數的計算是完全積性函式

$ P_i$表示線篩求出的第$ i$小的質數

令$ g(n,i)$表示$ \sum\limits_{j=2}^n [j的最小質因數>P_i或j是質數]f(j)$

在這裡$ f(j)$表示假設$ j$是質數,以質數方式帶入函式計算的結果

由於合數會被篩掉因而不會影響答案

考慮怎麼計算$ g(n,i)$

類似線性篩的方式每次篩掉一批合數

如果$P_i^2>n$則有$ g(n,i)=g(n,i-1)$

因為第$ i$個質數能篩掉的最小合數是$ P_i^2$

因此篩質數只需要篩到$ \sqrt n$即可

如果$ P_i^2<=n$有$ g(n,i)=g(n,i-1)-f(P_i)*(\ g( \frac{n}{P_i},i-1)-sp_{i-1}\ )$

原理是假設$ P_i$是一個質因數,它能產生的合數貢獻是$ f(P_i)*g( \frac{n}{P_i} ,i-1)$

但是由於$ P_i$不一定是最小質因數,還要加回多減的小質數即$ sp_{i-1}$

由於滿足$ f$是完全積性函式,上面部分還算挺清真的

我們需要求的只是$ g(x,INF)$

注意我們發現我們需要求的$ g(x,INF)$只需要滿足存在$ d$使得$ x=\frac{n}{d} $即可

可以提前整數分塊這樣只需要計算$ \sqrt n$數量級的$ g(x,INF)$即可

可以通過滾動陣列遞推的方式完成這一部分

 


 我們令$ S(n,m)$表示$ \sum\limits_{i=2}^n[i的最小質因數>=P_m]f(i)$

顯然我們要求的是$ S(n,1)$

遞迴求解

貢獻分兩步統計:

質數貢獻:$ g(n,INF)-sp(m-1)$

即去掉較小的質數以外其他質數都會被計算到

合數貢獻:$ \sum\limits_{k=m}^{P_k^2<=n}\sum\limits_{e=1}^{P_k^{e+1}<=n}f(P_k^e)S(\frac{n}{P_k^e},k+1)+f(P_k^{e+1})$

即列舉當前選擇的最小質因數以及數量轉移,同時計算只選擇多於兩個當前因數即不往後轉移的合數情況

這樣直接轉移就好了


栗子:篩$ \sum\limits_{i=1}^n \phi(i)$

發現$ phi(P_i)=P_i-1$即對於質數的計算不是一個完全積性函式

這時候需要拆開計算

令$ g(P_i)=P_i$,$ h(P_i)=1$

這樣分成兩個完全積性函式,分別篩質數求值然後相減即可

推$ S(n,m)$的時候不會有影響

篩$ \sum\limits_{i=1}^n \mu(i)$也沒有本質區別

傳送門: here

$ my \ code:$

#include<ctime>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<queue>
#define rt register unsigned
#define ll long long
#define M 100010
using namespace std;
inline ll read(){
    ll x = 0; char zf = 1; char ch = getchar();
    while (ch != '-' && !isdigit(ch)) ch = getchar();
    if (ch == '-') zf = -1, ch = getchar();
    while (isdigit(ch)) x = x * 10 + ch - '0', ch = getchar(); return x * zf;
}
void write(ll y){if(y<0)putchar('-'),y=-y;if(y>9)write(y/10);putchar(y%10+48);}
void writeln(const ll y){write(y);putchar('\n');}
int i,j,k,m,n,x,y,z,cnt;
ll sp[M],g[M];int h[M],ss[M];bool pri[M];
int id1[M],id2[M],q[M],t,sz;
void init(){
    sz=sqrt(n);
    for(rt i=2;i<=sz;i++){
        if(!pri[i])ss[++cnt]=i,sp[cnt]=sp[cnt-1]+i;
        for(rt j=1;j<=cnt&&i*ss[j]<=sz;j++){
            pri[i*ss[j]]=1;
            if(i%ss[j]==0)break;
        }
    }
    for(rt i=1;i<=n;){
        const unsigned v=n/i;unsigned R=n/v;
        q[++t]=v;if(v<=sz)id1[v]=t;else id2[n/v]=t;
        g[t]=(ll)v*(v+1)/2-1;
        h[t]=v-1;i=R+1;
    }
}
inline int id(int x){return x<=sz?id1[x]:id2[n/x];}
ll S(int n,int m){//phi
    if(n<=1||ss[m]>n)return 0;
    ll ret=g[id(n)]-sp[m-1]+m-1;
    for(rt k=m;ss[k]*ss[k]<=n&&k<=cnt;k++)
    for(rt v=ss[k],p1=ss[k]-1;(ll)v*ss[k]<=n;v=v*ss[k],p1=p1*ss[k])
    ret+=(S(n/v,k+1)+ss[k])*p1;
    return ret;
}

int D(int n,int m){//mu
    if(n<=1||ss[m]>n)return 0;
    int ret=h[id(n)]+(m-1);
    for(rt k=m;ss[k]*ss[k]<=n&&k<=cnt;k++)
    ret-=D(n/ss[k],k+1);
    return ret;
}
int main(){
    for(rt T=read();T;T--){
        n=read();t=cnt=0;init();
        for(rt i=1;i<=cnt;i++){
            const int v=ss[i]*ss[i];
            for(rt j=1;j<=t&&v<=q[j];j++){
                const int k=id(q[j]/ss[i]);
                g[j]-=ss[i]*(g[k]-sp[i-1]);
                h[j]-=h[k]-i+1;
            }
        }
        for(rt i=1;i<=t;i++)g[i]-=h[i],h[i]=-h[i];
        write(S(n,1)+1);putchar(' ');writeln(D(n,1)+1);
    }
    return 0;
}