Meisell-Lehmer算法(大素数模板)

来源：互联网发布：php截取html字符串编辑：程序博客网时间：2024/06/05 12:00

比赛的时候搜到了这个模板，居然没细看……看到其他博主有说还有更快的叫Deleglise Rivat的算法，并没有搜到……

努力看明白……

网上搜的模板：

#include <iostream>#include <cstdio>#include <cstring>#include <algorithm>#include <cmath>#define LL long longusing namespace std; const int N=5e6+2;bool np[N];int prime[N],pi[N];int getprime(){    int cnt=0;    np[0]=np[1]=true;    pi[0]=pi[1]=0;    for(int i=2;i<N;++i){        if(!np[i]) prime[++cnt]=i;        pi[i]=cnt;        for(int j=1;j<=cnt&&i*prime[j]<N;++j){            np[i*prime[j]]=true;            if(i%prime[j]==0) break;        }    }    return cnt;}const int M=7;const int PM=2*3*5*7*11*13*17;int phi[PM+1][M+1],sz[M+1];void init(){    getprime();    sz[0]=1;    for(int i=0;i<=PM;++i) phi[i][0]=i;    for(int i=1;i<=M;++i){        sz[i]=prime[i]*sz[i-1];        for(int j=1;j<=PM;++j){            phi[j][i]=phi[j][i-1]-phi[j/prime[i]][i-1];        }    }}int sqrt2(LL x){    LL r=(LL)sqrt(x-0.1);    while(r*r<=x) ++r;    return int(r-1);}int sqrt3(LL x) {    LL r=(LL)cbrt(x-0.1);    while(r*r*r<=x) ++r;    return int(r-1);}LL getphi(LL x,int s){    if(s==0) return x;    if(s<=M) return phi[x%sz[s]][s]+(x/sz[s])*phi[sz[s]][s];    if(x<=prime[s]*prime[s]) return pi[x]-s+1;    if(x<=prime[s]*prime[s]*prime[s]&&x<N){        int s2x=pi[sqrt2(x)];        LL ans=pi[x]-(s2x+s-2)*(s2x-s+1)/2;        for(int i=s+1;i<=s2x;++i){            ans+=pi[x/prime[i]];        }        return ans;    }    return getphi(x,s-1)-getphi(x/prime[s],s-1);}LL getpi(LL x){    if(x<N) return pi[x];    LL ans=getphi(x,pi[sqrt3(x)])+pi[sqrt3(x)]-1;    for(int i=pi[sqrt3(x)]+1,ed=pi[sqrt2(x)];i<=ed;++i){        ans-=getpi(x/prime[i])-i+1;    }    return ans;}LL lehmer_pi(LL x){    if(x<N) return pi[x];    int a=(int)lehmer_pi(sqrt2(sqrt2(x)));    int b=(int)lehmer_pi(sqrt2(x));    int c=(int)lehmer_pi(sqrt3(x));    LL sum=getphi(x,a)+LL(b+a-2)*(b-a+1)/2;    for(int i=a+1;i<=b;i++){        LL w=x/prime[i];        sum-=lehmer_pi(w);        if(i>c) continue;        LL lim=lehmer_pi(sqrt2(w));        for(int j=i;j<=lim;j++){            sum-=lehmer_pi(w/prime[j])-(j-1);        }    }    return sum;}int main() {    init();    LL n;    ios::sync_with_stdio(false);    while(cin>>n){        cout<<lehmer_pi(n)<<endl;    }    return 0;}

0 0