素数问题——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。主要可以通过这个题学会在用二分进行优化的方法。这个模板处理大范围数据还是很有用的