Count primes

  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 }