hdu 5901 Count primes (meisell-Lehmer)
Count primes
Time Limit: 12000/6000 MS (Java/Others) Memory Limit: 65536/65536 K (Java/Others)
Total Submission(s): 2625 Accepted Submission(s): 1337
Problem Description
Easy question! Calculate how many primes between [1...n]!
Input
Each line contain one integer n(1 <= n <= 1e11).Process to end of file.
Output
For each case, output the number of primes in interval [1...n]
Sample Input
2
3
10
Sample Output
1
2
4
1
2
4
C/C++:
1 #include <iostream> 2 #include <cstdio> 3 #include <cstring> 4 #include <algorithm> 5 #include <cmath> 6 using namespace std; 7 typedef long long ll; 8 const int N=5e6+2; 9 bool np[N]; 10 int prime[N],pi[N]; 11 int getprime() 12 { 13 int cnt=0; 14 np[0]=np[1]=true; 15 pi[0]=pi[1]=0; 16 for(int i=2; i<N; ++i) 17 { 18 if(!np[i]) 19 prime[++cnt]=i; 20 pi[i]=cnt; 21 for(int j=1; j<=cnt&&i*prime[j]<N; ++j) 22 { 23 np[i*prime[j]]=true; 24 if(i%prime[j]==0) break; 25 } 26 } 27 return cnt; 28 } 29 const int M=7; 30 const int PM=2*3*5*7*11*13*17; 31 int phi[PM+1][M+1],sz[M+1]; 32 void init() 33 { 34 getprime(); 35 sz[0]=1; 36 for(int i=0; i<=PM; ++i) 37 phi[i][0]=i; 38 for(int i=1; i<=M; ++i) 39 { 40 sz[i]=prime[i]*sz[i-1]; 41 for(int j=1; j<=PM; ++j) 42 phi[j][i]=phi[j][i-1]-phi[j/prime[i]][i-1]; 43 } 44 } 45 int sqrt2(ll x) 46 { 47 ll r=(ll)sqrt(x-0.1); 48 while(r*r<=x) ++r; 49 return int(r-1); 50 } 51 int sqrt3(ll x) 52 { 53 ll r=(ll)cbrt(x-0.1); 54 while(r*r*r<=x) ++r; 55 return int(r-1); 56 } 57 ll getphi(ll x,int s) 58 { 59 if(s==0) 60 return x; 61 if(s<=M) 62 return phi[x%sz[s]][s]+(x/sz[s])*phi[sz[s]][s]; 63 if(x<=prime[s]*prime[s]) 64 return pi[x]-s+1; 65 if(x<=prime[s]*prime[s]*prime[s]&&x<N) 66 { 67 int s2x=pi[sqrt2(x)]; 68 ll ans=pi[x]-(s2x+s-2)*(s2x-s+1)/2; 69 for(int i=s+1; i<=s2x; ++i) 70 ans+=pi[x/prime[i]]; 71 return ans; 72 } 73 return getphi(x,s-1)-getphi(x/prime[s],s-1); 74 } 75 ll getpi(ll x) 76 { 77 if(x<N) return pi[x]; 78 ll ans=getphi(x,pi[sqrt3(x)])+pi[sqrt3(x)]-1; 79 for(int i=pi[sqrt3(x)]+1,ed=pi[sqrt2(x)]; i<=ed; ++i) 80 ans-=getpi(x/prime[i])-i+1; 81 return ans; 82 } 83 ll lehmer_pi(ll x) 84 { 85 if(x<N) return pi[x]; 86 int a=(int)lehmer_pi(sqrt2(sqrt2(x))); 87 int b=(int)lehmer_pi(sqrt2(x)); 88 int c=(int)lehmer_pi(sqrt3(x)); 89 ll sum=getphi(x,a)+ll(b+a-2)*(b-a+1)/2; 90 for(int i=a+1; i<=b; i++) 91 { 92 ll w=x/prime[i]; 93 sum-=lehmer_pi(w); 94 if(i>c) 95 continue; 96 ll lim=lehmer_pi(sqrt2(w)); 97 for(int j=i; j<=lim; j++) 98 sum-=lehmer_pi(w/prime[j])-(j-1); 99 } 100 return sum; 101 } 102 int main() 103 { 104 init(); 105 ll n; 106 while(cin>>n) 107 cout<<lehmer_pi(n)<<endl; 108 return 0; 109 }