1e13内的素数个数 Meissel–Lehmer 算法

P7884 【模板】Meissel–Lehmer 算法

时间复杂度是玄学

模板准备

#include <bits/stdc++.h>
#define inf 0x7fffffff
#define ll long long
//#define int long long
//#define double long double
#define re register int
#define void inline void
#define eps 1e-18
//#define mod 1e9+7
#define ls(p) p<<1
#define rs(p) p<<1|1
//#define ls(p) e[p].l
//#define rs(p) e[p].r
//#define pi acos(-1.0)
#define pb push_back
//#define P pair < int , int >
#define mk make_pair
#define fi first
#define se second
using namespace std;
const int mod=998244353;
//const int M=632454+7;
//const int N=1911570+7;//?????????? 4e8.
const int N=20781012+7,M=6324554+7,K=50,P=2*3*5*7*11*13*17,Q=7,R=229;
ll sqrt_n,limit;
int pi_list[N],id1[M],id2[M],pre_phi[K+7][P+7],product[Q+7];
ll prime[N],number[M],big_pi_list[M];
bool p[N],vis[M];
double inv[N];
ll sqrt(ll n)
{
	ll ans=sqrt(1.0*n);
	while(ans*ans<=n)  ans++;
	return ans-1;
}
ll max(ll a, ll b)
{
	return a>b?a:b;
}
void init(ll n)
{
	int cnt=0,id=0;
	sqrt_n=sqrt(n);
	limit=max(1ll*sqrt_n*sqrt(log2(n)),R);
	p[0]=p[1]=1;
	pi_list[1]=0;
	for(ll i=2;i<=limit;i++)
	{
		if(!p[i])  prime[++cnt]=i,inv[cnt]=1.0/i+eps;
		pi_list[i]=cnt;
		for(re j=1;j<=cnt&&i*prime[j]<=limit;j++)
		{
			p[i*prime[j]]=1;
			if(i%prime[j]==0)  break;
		}
	}
	for(ll i=1,j;i<=n;i=j+1)
	{
		ll tn=n/i;
		j=n/tn;
		id++;
		number[id]=tn;
		if(tn<=sqrt_n)  id1[tn]=id;
		else  id2[j]=id;
	}
	for(re i=1;i<=P;i++)  pre_phi[0][i]=i;
	for(re i=1;i<=K;i++)  for(re j=1;j<=P;j++)  pre_phi[i][j]=pre_phi[i-1][j]-pre_phi[i-1][j/prime[i]];
	product[0]=1;
	for(re i=1;i<=Q;i++)  product[i]=product[i-1]*prime[i];
}
int get_id(ll n,ll m)
{
	return n<=sqrt_n?id1[n]:id2[m/n];
}
ll cbrt(ll n)
{
	ll ans=cbrt(1.0*n);
	while(ans*ans*ans<=n)  ans++;
	return ans-1;
}
ll sum1(int n)
{
	return 1ll*n*(n+1)/2;
}
ll pi(ll n,ll m);
ll p2(ll n,ll m,ll k)
{
	ll a=pi(sqrt(n),k),ans=sum1(m-1)-sum1(a-1);
	for(re i=m+1;i<=a;i++)  ans+=pi(n*inv[i],k);
	return ans;
}
ll phi(ll n,ll m,ll k)
{
	if(n==0)  return 0;
	if(m==0)  return n;
	if(n<=prime[m])  return 1;
	if(n<=P&&m<=K)  return pre_phi[m][n];
	if(m<=Q)  return n/product[m]*pre_phi[m][product[m]]+pre_phi[m][n%product[m]];
	if(n<prime[m]*prime[m])  return pi(n,k)-m+1;
	if(n<prime[m]*prime[m]*prime[m])  return pi(n,k)+p2(n,m,k)-m+1;
	return phi(n,m-1,k)-phi(n*inv[m],m-1,k);
}
ll pi(ll n,ll m)
{
	if(n<=limit)  return pi_list[n];
	int id=get_id(n,m);
	if(vis[id])  return big_pi_list[id];
	ll a=pi(cbrt(n),m);
	vis[id]=1;
	return big_pi_list[id]=phi(n,a,m)-p2(n,a,m)+a-1;
}
void solve()
{
	ll n;
	cin>>n;
	init(n);
	cout<<pi(n,n)<<endl;
}
signed main()
{
//	freopen("P1505_1.txt", "r", stdin);
//	freopen("Aout.txt", "w", stdout);
    int T=1;
//    cin>>T;
    for(int index=1;index<=T;index++)
    {
//        printf("Case #%lld: ",index);
        solve();
//        puts("");
    }
    return 0;
}
/*

10 5
hbtngdflmj
1 10 1
2 9 0
3 8 1
4 7 0
5 6 1


*/