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?

1. bool np[maxn];  
2. int prime[maxn],pi[maxn];  
3.   
4. int getprime()  
5. {  
6.     int cnt=0;  
7.     np[0]=np[1]=true;  
8.     pi[0]=pi[1]=0;  
9.     for(int i=2; i<maxn; ++i)  
10.     {  
11.         if(!np[i]) prime[++cnt]=i;  
12.         pi[i]=cnt;  
13.         for(int j=1; j<=cnt&&i*prime[j]<maxn; ++j)  
14.         {  
15.             np[i*prime[j]]=true;  
16.             if(i%prime[j]==0) break;  
17.         }  
18.     }  
19.     return cnt;  
20. }  
21. const int M=7;  
22. const int PM=2*3*5*7*11*13*17;  
23. int phi[PM+1][M+1],sz[M+1];  
24. void init()  
25. {  
26.     getprime();  
27.     sz[0]=1;  
28.     for(int i=0; i<=PM; ++i) phi[i][0]=i;  
29.     for(int i=1; i<=M; ++i)  
30.     {  
31.         sz[i]=prime[i]*sz[i-1];  
32.         for(int j=1; j<=PM; ++j)  
33.         {  
34.             phi[j][i]=phi[j][i-1]-phi[j/prime[i]][i-1];  
35.         }  
36.     }  
37. }  
38. int sqrt2(ll x)  
39. {  
40.     ll r=(ll)sqrt(x-0.1);  
41.     while(r*r<=x) ++r;  
42.     return int(r-1);  
43. }  
44. int sqrt3(ll x)  
45. {  
46.     ll r=(ll)cbrt(x-0.1);  
47.     while(r*r*r<=x) ++r;  
48.     return int(r-1);  
49. }  
50. ll getphi(ll x,int s)  
51. {  
52.     if(s==0) return x;  
53.     if(s<=M) return phi[x%sz[s]][s]+(x/sz[s])*phi[sz[s]][s];  
54.     if(x<=prime[s]*prime[s]) return pi[x]-s+1;  
55.     if(x<=prime[s]*prime[s]*prime[s]&&x<maxn)  
56.     {  
57.         int s2x=pi[sqrt2(x)];  
58.         ll ans=pi[x]-(s2x+s-2)*(s2x-s+1)/2;  
59.         for(int i=s+1; i<=s2x; ++i)  
60.         {  
61.             ans+=pi[x/prime[i]];  
62.         }  
63.         return ans;  
64.     }  
65.     return getphi(x,s-1)-getphi(x/prime[s],s-1);  
66. }  
67. ll getpi(ll x)  
68. {  
69.     if(x<maxn) return pi[x];  
70.     ll ans=getphi(x,pi[sqrt3(x)])+pi[sqrt3(x)]-1;  
71.     for(int i=pi[sqrt3(x)]+1,ed=pi[sqrt2(x)]; i<=ed; ++i)  
72.     {  
73.         ans-=getpi(x/prime[i])-i+1;  
74.     }  
75.     return ans;  
76. }  
77. ll lehmer_pi(ll x)  
78. {  
79.     if(x<maxn) return pi[x];  
80.     int a=(int)lehmer_pi(sqrt2(sqrt2(x)));  
81.     int b=(int)lehmer_pi(sqrt2(x));  
82.     int c=(int)lehmer_pi(sqrt3(x));  
83.     ll sum=getphi(x,a)+ll(b+a-2)*(b-a+1)/2;  
84.     for(int i=a+1; i<=b; i++)  
85.     {  
86.         ll w=x/prime[i];  
87.         sum-=lehmer_pi(w);  
88.         if(i>c) continue;  
89.         ll lim=lehmer_pi(sqrt2(w));  
90.         for(int j=i; j<=lim; j++)  
91.         {  
92.             sum-=lehmer_pi(w/prime[j])-(j-1);  
93.         }  
94.     }  
95.     return sum;  
96. }  

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

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

[html]  view plain  copy

  print ?

1. int main()  
2. {  
3.     init();  
4.     ll n;  
5.     while(scanf("%I64d",&n))  
6.     {  
7.         printf("%I64d\n",lehmer_pi(n))  
8.     }  
9.     return 0;  
10. }  

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

[html]  view plain  copy

 print ?

1. #include <cstdio>  
2. #include <cstring>  
3. #include <cmath>  
4. #include <iostream>  
5. using namespace std;  
6. #define mst(a,b) memset((a),(b),sizeof(a))  
7. #define f(i,a,b) for(int i=(a);i<(b);++i)  
8. typedef long long ll;  
9. const int maxn= 10005;  
10. const int mod = 10007;  
11. const ll INF = 0x3f3f3f3f;  
12. const double eps = 1e-6;  
13. #define rush() int T;scanf("%d",&T);while(T--)  
14.   
15. bool np[maxn];  
16. int prime[maxn],pi[maxn];  
17.   
18. int getprime()  
19. {  
20.     int cnt=0;  
21.     np[0]=np[1]=true;  
22.     pi[0]=pi[1]=0;  
23.     for(int i=2; i<maxn; ++i)  
24.     {  
25.         if(!np[i]) prime[++cnt]=i;  
26.         pi[i]=cnt;  
27.         for(int j=1; j<=cnt&&i*prime[j]<maxn; ++j)  
28.         {  
29.             np[i*prime[j]]=true;  
30.             if(i%prime[j]==0) break;  
31.         }  
32.     }  
33.     return cnt;  
34. }  
35. const int M=7;  
36. const int PM=2*3*5*7*11*13*17;  
37. int phi[PM+1][M+1],sz[M+1];  
38. void init()  
39. {  
40.     getprime();  
41.     sz[0]=1;  
42.     for(int i=0; i<=PM; ++i) phi[i][0]=i;  
43.     for(int i=1; i<=M; ++i)  
44.     {  
45.         sz[i]=prime[i]*sz[i-1];  
46.         for(int j=1; j<=PM; ++j)  
47.         {  
48.             phi[j][i]=phi[j][i-1]-phi[j/prime[i]][i-1];  
49.         }  
50.     }  
51. }  
52. int sqrt2(ll x)  
53. {  
54.     ll r=(ll)sqrt(x-0.1);  
55.     while(r*r<=x) ++r;  
56.     return int(r-1);  
57. }  
58. int sqrt3(ll x)  
59. {  
60.     ll r=(ll)cbrt(x-0.1);  
61.     while(r*r*r<=x) ++r;  
62.     return int(r-1);  
63. }  
64. ll getphi(ll x,int s)  
65. {  
66.     if(s==0) return x;  
67.     if(s<=M) return phi[x%sz[s]][s]+(x/sz[s])*phi[sz[s]][s];  
68.     if(x<=prime[s]*prime[s]) return pi[x]-s+1;  
69.     if(x<=prime[s]*prime[s]*prime[s]&&x<maxn)  
70.     {  
71.         int s2x=pi[sqrt2(x)];  
72.         ll ans=pi[x]-(s2x+s-2)*(s2x-s+1)/2;  
73.         for(int i=s+1; i<=s2x; ++i)  
74.         {  
75.             ans+=pi[x/prime[i]];  
76.         }  
77.         return ans;  
78.     }  
79.     return getphi(x,s-1)-getphi(x/prime[s],s-1);  
80. }  
81. ll getpi(ll x)  
82. {  
83.     if(x<maxn) return pi[x];  
84.     ll ans=getphi(x,pi[sqrt3(x)])+pi[sqrt3(x)]-1;  
85.     for(int i=pi[sqrt3(x)]+1,ed=pi[sqrt2(x)]; i<=ed; ++i)  
86.     {  
87.         ans-=getpi(x/prime[i])-i+1;  
88.     }  
89.     return ans;  
90. }  
91. ll lehmer_pi(ll x)  
92. {  
93.     if(x<maxn) return pi[x];  
94.     int a=(int)lehmer_pi(sqrt2(sqrt2(x)));  
95.     int b=(int)lehmer_pi(sqrt2(x));  
96.     int c=(int)lehmer_pi(sqrt3(x));  
97.     ll sum=getphi(x,a)+ll(b+a-2)*(b-a+1)/2;  
98.     for(int i=a+1; i<=b; i++)  
99.     {  
100.         ll w=x/prime[i];  
101.         sum-=lehmer_pi(w);  
102.         if(i>c) continue;  
103.         ll lim=lehmer_pi(sqrt2(w));  
104.         for(int j=i; j<=lim; j++)  
105.         {  
106.             sum-=lehmer_pi(w/prime[j])-(j-1);  
107.         }  
108.     }  
109.     return sum;  
110. }  
111.   
112. int main()  
113. {  
114.     init();  
115.     ll n;  
116.     int cas=1;  
117.     while(~scanf("%lld",&n)&&n)  
118.     {  
119.         ll l=2,r=5e7+10;  //r为估算的最大值  
120.         while(l<r)  
121.         {  
122.             ll m=(l+r)/2;  
123.             ll res=lehmer_pi(m);  
124.             if(res>=n)  
125.                 r=m;  
126.             else l=m+1;  
127.         }  
128.         printf("Case %d: %lld\n",cas++,l);  
129.     }  
130.     return 0;  
131. }  

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

贴上官方题解

[cpp] view plain copy

 print?

1. #include <bits/stdtr1c++.h>  
2.    
3. #define MAXN 100  
4. #define MAXM 10001  
5. #define MAXP 40000  
6. #define MAX 400000  
7. #define clr(ar) memset(ar, 0, sizeof(ar))  
8. #define read() freopen("lol.txt", "r", stdin)  
9. #define dbg(x) cout << #x << " = " << x << endl  
10. #define chkbit(ar, i) (((ar[(i) >> 6]) & (1 << (((i) >> 1) & 31))))  
11. #define setbit(ar, i) (((ar[(i) >> 6]) |= (1 << (((i) >> 1) & 31))))  
12. #define isprime(x) (( (x) && ((x)&1) && (!chkbit(ar, (x)))) || ((x) == 2))  
13.    
14. using namespace std;  
15.    
16. namespace pcf{  
17.     long long dp[MAXN][MAXM];  
18.     unsigned int ar[(MAX >> 6) + 5] = {0};  
19.     int len = 0, primes[MAXP], counter[MAX];  
20.    
21.     void Sieve(){  
22.         setbit(ar, 0), setbit(ar, 1);  
23.         for (int i = 3; (i * i) < MAX; i++, i++){  
24.             if (!chkbit(ar, i)){  
25.                 int k = i << 1;  
26.                 for (int j = (i * i); j < MAX; j += k) setbit(ar, j);  
27.             }  
28.         }  
29.    
30.         for (int i = 1; i < MAX; i++){  
31.             counter[i] = counter[i - 1];  
32.             if (isprime(i)) primes[len++] = i, counter[i]++;  
33.         }  
34.     }  
35.    
36.     void init(){  
37.         Sieve();  
38.         for (int n = 0; n < MAXN; n++){  
39.             for (int m = 0; m < MAXM; m++){  
40.                 if (!n) dp[n][m] = m;  
41.                 else dp[n][m] = dp[n - 1][m] - dp[n - 1][m / primes[n - 1]];  
42.             }  
43.         }  
44.     }  
45.    
46.     long long phi(long long m, int n){  
47.         if (n == 0) return m;  
48.         if (primes[n - 1] >= m) return 1;  
49.         if (m < MAXM && n < MAXN) return dp[n][m];  
50.         return phi(m, n - 1) - phi(m / primes[n - 1], n - 1);  
51.     }  
52.    
53.     long long Lehmer(long long m){  
54.         if (m < MAX) return counter[m];  
55.    
56.         long long w, res = 0;  
57.         int i, a, s, c, x, y;  
58.         s = sqrt(0.9 + m), y = c = cbrt(0.9 + m);  
59.         a = counter[y], res = phi(m, a) + a - 1;  
60.         for (i = a; primes[i] <= s; i++) res = res - Lehmer(m / primes[i]) + Lehmer(primes[i]) - 1;  
61.         return res;  
62.     }  
63. }  
64.    
65. long long solve(long long n){  
66.     int i, j, k, l;  
67.     long long x, y, res = 0;  
68.    
69.     for (i = 0; i < pcf::len; i++){  
70.         x = pcf::primes[i], y = n / x;  
71.         if ((x * x) > n) break;  
72.         res += (pcf::Lehmer(y) - pcf::Lehmer(x));  
73.     }  
74.    
75.     for (i = 0; i < pcf::len; i++){  
76.         x = pcf::primes[i];  
77.         if ((x * x * x) > n) break;  
78.         res++;  
79.     }  
80.    
81.     return res;  
82. }  
83.    
84. int main(){  
85.     pcf::init();  
86.     long long n, res,L,R,M,ca=1;  
87.     while (scanf("%lld", &n) != EOF){  
88.         if(n==0) break;  
89.         L=2;R=1e8;  
90.    
91.         while(L<R){  
92.             M=(L+R)/2;  
93.             res=pcf::Lehmer(M);  
94.             if(res>=n) R=M;  
95.             else L=M+1;  
96.         }  
97.         printf("Case %lld: %lld\n",ca++,L);  
98.    
99.     }  
100.     return 0;  
101. }  

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