Nth prime & numbers of primes (模板)
都是取的模板,这几天做的素数题挺多的,所以整理了放在这里,感觉有一天回用到的!
SPOJ:Nth Prime: 求第N个素数,N<1e9。
#include<bits/stdc++.h> using namespace std; typedef long long ll; const int N=2000000,P=510510,Q=92160; struct getnthprime { int prime[N+1],pi[N+1],e[P]; void init(){ for(int i=2;i<=N;i++) { if(!prime[i]) prime[++prime[0]]=i,pi[i]=pi[i-1]+1; else pi[i]=pi[i-1]; for(int j=1;j<=prime[0]&&i<=N/prime[j];j++) { prime[i*prime[j]]=1; if(i%prime[j]==0) break; } } for(int i=0;i<P;i++) e[i]=i; for(int i=1;i<=7;i++) { for(int j=P-1;j>=0;j--) e[j]-=e[j/prime[i]]; } } ll get_phi(ll m,int n) { if (n==7) return m/P*Q+e[m%P]; if (m<prime[n]) return 1; if (m<=N&&m<=(ll)prime[n]*prime[n]*prime[n]) { ll ans=pi[m]-n+1; for(int i=n+1,l=pi[(int)sqrt(m+0.1)];i<=l;i++) ans+=pi[m/prime[i]]-i+1; return ans; } return get_phi(m,n-1)-get_phi(m/prime[n],n-1); } ll get_pi(ll m){ if(m<=N) return pi[m]; int n=pi[(int)cbrt(m-0.1)+1]; ll ans=get_phi(m,n)+n-1; for(int i=n+1,l=pi[(int)sqrt(m+0.1)];i<=l;i++) ans-=get_pi(m/prime[i])-i+1; return ans; } bool f[1000010]; ll get_pn(ll n) { if (n<=prime[0]) return prime[n]; ll x=n*(log(n)+log(log(n))-1)+n*(log(log(n))-2)/log(n)-6*n/1000; ll y=n*(log(log(n)))*(log(log(n)))/log(n)/log(n); y=min(y,ll(3500000)); ll l=x,r=x+y,flag = 0; for (int i=0;i<2;i++) { ll m=(l+r)>> 1; ll pm=get_pi(m); if(pm>=n) r=m,flag=0; else l=m+1,flag=pm; } ll count=flag?flag:get_pi(l-1); for(int i=1,li=pi[(int)sqrt(r+0.1)];i<=li;i++) { for(int j=((l-1)/prime[i]+1)*prime[i]-l;j<=r-l+1; j+=prime[i]){ f[j]=true; } } for(int i=0;i<=r-l+1;i++) { if(!f[i]){ count++; if(count==n) return i+l; } } return -1; } }NP; int main() { NP.init(); ll n; scanf("%lld",&n); cout<<NP.get_pn(n)<<endl; return 0; }
HDU5901:Count primes: 求1到N有多少个素数。N<1e11。
#include<bits/stdc++.h> using namespace std; #define ll long long const int N=5e6+2,M=7,PM=2*3*5*7*11*13*17; struct countprimes { bool np[N],did[N]; int prime[N],pi[N],phi[PM+1][M+1],sz[M+1]; vector<ll>v; int getprime() { int cnt=0; np[0]=np[1]=true; pi[0]=pi[1]=0; for(int i=2;i<N;++i){ if(!np[i]) prime[++cnt]=i; pi[i]=cnt; for(int j=1;j<=cnt&&i*prime[j]<N;++j){ np[i*prime[j]]=true; if(i%prime[j]==0) break; } } return cnt; } 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<N) { 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<N) 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<N) 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; } }CP; int main() { CP.init(); ll n,ans=0; while(~scanf("%lld",&n)){ cout<<CP.lehmer_pi(n)<<endl; } }
It is your time to fight!