素数问题——Meisell-Lehmer算法

我觉得素数问题不是一两篇文章就可以解释清楚的，在数论中很多相关问题至今都无法解决。这里只是简单的总结在竞赛中用到的一小部分。

HDU5901这道是网络赛的一道素数题，有多种解法，这里推荐一些大佬的博客：http://blog.csdn.net/chaiwenjun000/article/details/52589457

http://blog.csdn.net/angon823/article/details/52577145（有资料链接）

素数筛有很多种方法，比如欧拉筛，埃氏筛，Meisell-Lehmer算法等等

这里主要说一下Meisell-Lehmer算法

Meisell-Lehmer算法模板（计算2~n中的素数个数）

bool np[maxn];  int prime[maxn],pi[maxn];    int getprime()  {      int cnt=0;      np[0]=np[1]=true;      pi[0]=pi[1]=0;      for(int i=2; i<maxn; ++i)      {          if(!np[i]) prime[++cnt]=i;          pi[i]=cnt;          for(int j=1; j<=cnt&&i*prime[j]<maxn; ++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<maxn)      {          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<maxn) 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<maxn) 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;  }  

这种算法在时间和数据量上有很大优势，但是在内存上有一定的不足

改模板的时候，原来的七个数可以取前五个，这里会消耗一些内存

在之江学院的预赛中出了卡内存和时间的素数题，当时过这个题要优化模板还得加上二分

这里有必要研究一下

题目链接http://115.231.222.240:8081/JudgeOnline/problem.php?id=1489

米勒罗宾算法主要是求从1到n内有多少素数的，所以可以用二分优化，上限是估算的最大值

代码：

#include <iostream>#include <cstdio>#include <cstring>#include <cmath>using namespace std;typedef long long ll;const int maxn= 10005;bool np[maxn];int prime[maxn],pi[maxn];int getprime(){    int cnt=0;    np[0]=np[1]=true;    pi[0]=pi[1]=0;    for(int i=2; i<maxn; ++i)    {        if(!np[i]) prime[++cnt]=i;        pi[i]=cnt;        for(int j=1; j<=cnt&&i*prime[j]<maxn; ++j)        {            np[i*prime[j]]=true;            if(i%prime[j]==0) break;        }    }    return cnt;}const int M=5;const int PM=2*3*5*7*11;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<maxn)    {        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<maxn) 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<maxn) 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;    int cas=1;    while(~scanf("%lld",&n)&&n) {        ll l=2,r=5e7+10;         while(l<r){            ll m=(l+r)/2;            ll res=lehmer_pi(m);            if(res>=n)                r=m;            else l=m+1;        }        printf("Case %d: %lld\n",cas++,l);    }    return 0;}

说明：这个题内存限制在16M，但是数据是在三百万之内的，所以可以对模板进行阉割的。杭电上的求10的11次方的素数内存是512M。主要可以通过这个题学会在用二分进行优化的方法。这个模板处理大范围数据还是很有用的