P7884 【模板】Meissel–Lehmer 算法（亚线性求素数个数）

题目链接：点击这里

题目大意：
给定整数 $n$ ，求出 $\pi (n)$ 的值。
其中 $\pi (n)$ 表示 $1\sim n$ 的整数中质数的个数。

题目分析：
朴素的 $min25$ 筛在求素数个数上被此算法完爆（此题 $min25$ 求 $10^{13}$ 的数据会 $t$ 飞）
原理不明，存个模板
$\text{Meissel–Lehmer}$ 算法时间复杂度 $10^{13}$ 以内为 $O(可过)$
PS：树状数组优化 $min25$ 筛可通过此题，待补…

具体细节见代码：

//#pragma GCC optimize(2)
//#pragma GCC optimize("Ofast","inline","-ffast-math")
//#pragma GCC target("avx,sse2,sse3,sse4,mmx")
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<vector>
#include<set>
#include<map>
#include<stack>
#include<queue>
#define ll long long
#define inf 0x3f3f3f3f
#define Inf 0x3f3f3f3f3f3f3f3f
#define endl '\n'
#define IOS ios::sync_with_stdio(0); cin.tie(0); cout.tie(0)
using namespace std;
int read()
{
	int res = 0,flag = 1;
	char ch = getchar();
	while(ch<'0' || ch>'9')
	{
		if(ch == '-') flag = -1;
		ch = getchar();
	}
	while(ch>='0' && ch<='9')
	{
		res = (res<<3)+(res<<1)+(ch^48);//res*10+ch-'0';
		ch = getchar();
	}
	return res*flag;
}

const int N = 3e7 + 7, M = 7e6 + 7, K = 50, P = 2 * 3 * 5 * 7 * 11 * 13 * 17, Q = 7, R = 229;
// const int N = 20781012 + 7, M = 6324554 + 7, K = 50, P = 2 * 3 * 5 * 7 * 11 * 13 * 17, Q = 7, R = 229;
const double eps = 1e-18;
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];
double inv[N];
bool p[N], vis[M];

inline ll sqrt(ll n)
{
	ll ans = sqrt((double)n);
	while (ans * ans <= n) ans++;
	return ans - 1;
}

inline ll max(ll a, ll b)
{
	return a > b ? a : b;
}

inline void init(ll n){
	int cnt = 0, id = 0;
	sqrt_n = sqrt(n);
	limit = max(sqrt_n * sqrt(log2(n)), R);
	p[0] = p[1] = true;
	pi_list[1] = 0;
	for (register ll i = 2; i <= limit; i++)
	{
		if (!p[i])
		{
			cnt++;
			prime[cnt] = i;
			inv[cnt] = 1.0 / i + eps;
		}
		pi_list[i] = cnt;
		for (register int j = 1; j <= cnt && i * prime[j] <= limit; j++)
		{
			p[i * prime[j]] = true;
			if (i % prime[j] == 0) break;
		}
	}
	for (register 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 (register int i = 1; i <= P; i++)
		pre_phi[0][i] = i;
	for (register int i = 1; i <= K; i++)
		for (register int 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 (register int i = 1; i <= Q; i++) 
		product[i] = product[i - 1] * prime[i];
}

inline int get_id(ll n, ll m)
{
	return n <= sqrt_n ? id1[n] : id2[m / n];
}

inline ll cbrt(ll n)
{
	ll ans = cbrt((double)n);
	while (ans * ans * ans <= n) ans++;
	return ans - 1;
}

ll pi(ll n, ll m);

inline ll sum1(ll n)
{
	return n * (n + 1) / 2;
}

inline ll p2(ll n, ll m, ll k)
{
	ll a = pi(sqrt(n), k), ans = sum1(m - 1) - sum1(a - 1);
	for (register ll 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 <= 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] = true;
	return big_pi_list[id] = phi(n, a, m) - p2(n, a, m) + a - 1;
}

int main()
{
	ll n;
	scanf("%lld", &n);
	init(n);
	printf("%lld", pi(n, n));
	return 0;
}