ZJCOJ L先生与质数V3/V4 (Meisell-Lehmer算法)

Problem L: L先生与质数V3/V4（应各位菊苣要求）

Time Limit:1
Sec Memory Limit: 128/16 MB

Submit:298 Solved:65

[Submit][Status][WebBoard]

Description

在解决了上一个质数问题之后，L先生依然不甘心，他还想计算下更多范围内的质数，你能帮助他吗？（没错这题题面和V3一毛一样）

Input

有多组测试例。（测试例数量<70)

每个测试例一行，输入一个数字n（0<n<=3000000），输入0表示结束。

Output

输出测试例编号和第N个质数。

Case X: Y

Sample Input

1

2

3

4

10

100

0

Sample Output

Case 1: 2

Case 2: 3

Case 3: 5

Case 4: 7

Case 5: 29

Case 6: 541

（如果要看对小范围的数据进行素数筛法的操作，请看这里：传送门）

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

[html] view
plain copy

print?

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;

}

其中maxn值的大小取为sqrt（maxn（n））

求得1~n中素数个数的操作即为

[html] view
plain copy

print?

int main()

{

init();

ll n;

while(scanf("%I64d",&n))

{

printf("%I64d\n",lehmer_pi(n))

}

return 0;

}

有了上面知识的铺垫后，容易发现我们要完成上面那道题只要用Meisell-Lehmer+二分搜索答案即可

[html] view
plain copy

print?

#include <cstdio>

#include <cstring>

#include <cmath>

#include <iostream>

using namespace std;

#define mst(a,b) memset((a),(b),sizeof(a))

#define f(i,a,b) for(int i=(a);i<(b);++i)

typedef long long ll;

const int maxn= 10005;

const int mod = 10007;

const ll INF = 0x3f3f3f3f;

const double eps = 1e-6;

#define rush() int T;scanf("%d",&T);while(T--)

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;

}

int main()

{

init();

ll n;

int cas=1;

while(~scanf("%lld",&n)&&n)

{

ll l=2,r=5e7+10; //r为估算的最大值

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;

}

（注：由于phi数组消耗内存很大，依然过不了V4（好像可以用米勒罗宾随机算法，学好后补上），但学到了Meisell-Lehmer算法+二分的方法）







贴上官方题解

[cpp] view
plain copy

print?

#include <bits/stdtr1c++.h>

#define MAXN 100

#define MAXM 10001

#define MAXP 40000

#define MAX 400000

#define clr(ar) memset(ar, 0, sizeof(ar))

#define read() freopen("lol.txt", "r", stdin)

#define dbg(x) cout << #x << " = " << x << endl

#define chkbit(ar, i) (((ar[(i) >> 6]) & (1 << (((i) >> 1) & 31))))

#define setbit(ar, i) (((ar[(i) >> 6]) |= (1 << (((i) >> 1) & 31))))

#define isprime(x) (( (x) && ((x)&1) && (!chkbit(ar, (x)))) || ((x) == 2))

using namespace std;

namespace pcf{

long long dp[MAXN][MAXM];

unsigned int ar[(MAX >> 6) + 5] = {0};

int len = 0, primes[MAXP], counter[MAX];

void Sieve(){

setbit(ar, 0), setbit(ar, 1);

for (int i = 3; (i * i) < MAX; i++, i++){

if (!chkbit(ar, i)){

int k = i << 1;

for (int j = (i * i); j < MAX; j += k) setbit(ar, j);

}

}

for (int i = 1; i < MAX; i++){

counter[i] = counter[i - 1];

if (isprime(i)) primes[len++] = i, counter[i]++;

}

}

void init(){

Sieve();

for (int n = 0; n < MAXN; n++){

for (int m = 0; m < MAXM; m++){

if (!n) dp
[m] = m;

else dp
[m] = dp[n - 1][m] - dp[n - 1][m / primes[n - 1]];

}

}

}

long long phi(long long m, int n){

if (n == 0) return m;

if (primes[n - 1] >= m) return 1;

if (m < MAXM && n < MAXN) return dp
[m];

return phi(m, n - 1) - phi(m / primes[n - 1], n - 1);

}

long long Lehmer(long long m){

if (m < MAX) return counter[m];

long long w, res = 0;

int i, a, s, c, x, y;

s = sqrt(0.9 + m), y = c = cbrt(0.9 + m);

a = counter[y], res = phi(m, a) + a - 1;

for (i = a; primes[i] <= s; i++) res = res - Lehmer(m / primes[i]) + Lehmer(primes[i]) - 1;

return res;

}

}

long long solve(long long n){

int i, j, k, l;

long long x, y, res = 0;

for (i = 0; i < pcf::len; i++){

x = pcf::primes[i], y = n / x;

if ((x * x) > n) break;

res += (pcf::Lehmer(y) - pcf::Lehmer(x));

}

for (i = 0; i < pcf::len; i++){

x = pcf::primes[i];

if ((x * x * x) > n) break;

res++;

}

return res;

}

int main(){

pcf::init();

long long n, res,L,R,M,ca=1;

while (scanf("%lld", &n) != EOF){

if(n==0) break;

L=2;R=1e8;

while(L<R){

M=(L+R)/2;

res=pcf::Lehmer(M);

if(res>=n) R=M;

else L=M+1;

}

printf("Case %lld: %lld\n",ca++,L);

}

return 0;

}

顺便给出可以HDOJ上可以用Meisell-Lehmer算法过的一道题的链接：HDOJ
5901 Count primes