筛法

埃筛

求 $\pi(n)$，$n \le 10^8$。

对于每一个质数进行一轮筛。

时间复杂度 $O(n \ln \ln n)$。

const int N = 1e8 + 5;
int n; bitset<N> ip;
void solve() {
	cin >> n;
	ip[1] = 1;
	FOR(i, 2, n) if(!ip[i]) 
		for(int j = i << 1; j <= n; j += i) 
			ip[j] = 1;
	cout << n - ip.count() << endl;
}

线性筛

求 $\pi(n)$，$n \le 10^8$。

将每个数筛出其最小的质因数。

时间复杂度 $O(n)$。

const int N = 1e8 + 5;
int n, p[N / 10], cnt;
bitset<N> ip;
void solve() {
	cin >> n;
	ip[1] = 1;
	FOR(i, 2, n) {
		if(!ip[i]) {
			p[++cnt] = i;
		}
		for(int j = 1; j <= cnt && p[j] * i <= n; j++) {
			ip[p[j] * i] = 1;
			if(i % p[j] == 0) break;
		}
	}
	cout << cnt << endl;
}

由于线性筛是筛出其最小的质因数，所以只要 $f(p^k)$ 能快速求出，线性筛就可以在 $O(n)$ 筛出任何积性函数。

筛 $\varphi(n)$ 和 $\mu(n)$，$n \le 10^7$。

只要知道 $f(p^k)$ 即可：

$\varphi(p)=p-1$
$\varphi(p^k)=(p-1)p^{k-1}=\varphi(p^{k-1})p~(k>1)$

$\mu(p)=-1$
$\mu(p^k)=0~(k>1)$

时间复杂度 $O(n)$。

const int N = 1e7 + 5;
int n, p[N], ip[N], phi[N], mu[N], cnt;
void solve() {
	cin >> n;
	mu[1] = phi[1] = ip[1] = 1;
	FOR(i, 2, n) {
		if(!ip[i]) {
			p[++cnt] = i;
			phi[i] = i - 1;
			mu[i] = -1;
		}
		for(int j = 1; j <= cnt && p[j] * i <= n; j++) {
			ip[p[j] * i] = 1;
			if(i % p[j] == 0) {
				phi[p[j] * i] = phi[i] * p[j];
				break;	
			}
			else {
				phi[p[j] * i] = phi[i] * (p[j] - 1);
				mu[p[j] * i] = - mu[i];
			}
		}
	}
	FOR(i, 1, n) cout << i << " " << phi[i] << " " << mu[i] << endl;
}

杜教筛

求 $\sum_{i=1}^{n}\varphi(i)$ 和 $\sum_{i=1}^{n}\mu(i)$，$n \le 10^{10}$

首先考虑积性函数 $f(n)$ 的前缀和函数 $s(n)$，构造函数 $g(n)$，考虑 $(f \ast g)$ 的前缀和：

$\sum_{i=1}^{n}(f \ast g)(i)$
$=\sum_{i=1}^{n}\sum_{d|i}g(d)f(\frac{i}{d})$
$=\sum_{i=1}^{n}g(i)\sum_{ij\le n}f(j)$
$=\sum_{i=1}^{n}g(i)s(\left \lfloor \frac{n}{i} \right \rfloor )$

把 $i=1$ 差出来：

$g(1)s(n)$
$=\sum_{i=1}^{n}g(i)s(\left \lfloor \frac{n}{i} \right \rfloor )-\sum_{i=2}^{n}g(i)s(\left \lfloor \frac{n}{i} \right \rfloor )$
$=\sum_{i=1}^{n}(f \ast g)(i)-\sum_{i=2}^{n}g(i)s(\left \lfloor \frac{n}{i} \right \rfloor )$

对后面的式子进行整除分块和记忆化即可，可以做到 $O(n^\frac{3}{4})$。

若线性筛预处理到 $n^\frac{2}{3}$，则可以做到 $O(n^\frac{2}{3})$。

先构造 $g(n)$。因为 $I \ast \varphi=id$ 和 $I \ast \mu=\epsilon $，所以直接构造 $g=I$。

时间复杂度 $O(n^\frac{2}{3})$。

const int N = 5e6 + 5;
int p[N], ip[N], phi[N], mu[N], cnt;
ll sphi[N]; int smu[N];
map<ll, ll> mphi;
map<ll, int> mmu;
void init() {
	mu[1] = phi[1] = ip[1] = 1;
	FOR(i, 2, N - 1) {
		if(!ip[i]) {
			p[++cnt] = i;
			phi[i] = i - 1;
			mu[i] = -1;
		}
		for(int j = 1; j <= cnt && p[j] * i < N; j++) {
			ip[p[j] * i] = 1;
			if(i % p[j] == 0) {
				phi[p[j] * i] = phi[i] * p[j];
				break;	
			}
			else {
				phi[p[j] * i] = phi[i] * (p[j] - 1);
				mu[p[j] * i] = - mu[i];
			}
		}
	}
	FOR(i, 1, N - 1) sphi[i] = sphi[i - 1] + phi[i];
	FOR(i, 1, N - 1) smu[i] = smu[i - 1] + mu[i];
}
ll get_phi(ll n) {
	if(n < N) return sphi[n];
	if(mphi.count(n)) return mphi[n];
	ll res = n * (n + 1) / 2;
	for(ll l = 2, r; l <= n; l = r + 1) {
		r = n / (n / l);
		res -= (r - l + 1) * get_phi(n / l);
	}
	return mphi[n] = res;
}
int get_mu(ll n) {
	if(n < N) return smu[n];
	if(mmu.count(n)) return mmu[n];
	int res = 1;
	for(ll l = 2, r; l <= n; l = r + 1) {
		r = n / (n / l);
		res -= (r - l + 1) * get_mu(n / l);
	}
	return mmu[n] = res;
}
void solve() {
	ll n; cin >> n;
	cout << get_phi(n) << " " << get_mu(n) << endl;
}

Meissel-Lehmer

摘自 ケロシの题解：P7884 【模板】Meissel-Lehmer

求 $\pi(n)$，$n \le 10^{13}$。

一道 Meissel-Lehmer 算法的模板题洛谷题解里没有 Meissel-Lehmer，大悲啊。

这里介绍 Meissel-Lehmer 算法。

先线性筛算出前大于等于 $n^\frac{1}{2}$ 个数内的素数和 $\pi(n)$。

然后对于求 $\pi(n)$，考虑对其产生贡献的素数分治。 考虑对于 $n^\frac{1}{3}$ 分治，先算出 $(n^\frac{1}{3},n]$ 中对于 $n^\frac{1}{3}$ 以内的素数筛不掉的数，暂定一个 dp，$dp_{i,j}$ 表示 $[1,i]$ 中不会被前 $j$ 个素数筛掉的个数。

接下来不难发现，现在筛剩下的 $(n^\frac{1}{3},n]$ 中的数，如果是合数，则其必定为两个大于 $n^\frac{1}{3}$ 的素数的乘积，其他即为素数。

计算这类由两个素数相乘的数的个数很简单，直接考虑枚举小的素数，去计算另外一个素数的方案数即可，假设当前选取了小素数 $p$，则其筛掉的数就是 $[p,\frac{n}{p}]$ 中的素数的个数，即为 $\pi(\frac{n}{p})-\pi(p)+1$。

结合一下上述，可得出 $\pi(n)$ 的递归公式：
$$ \pi(n)=dp_{n,\pi(n^\frac{1}{3})}-1+\pi(n^\frac{1}{3}) -\sum_{p\in(n^\frac{1}{3},n^\frac{1}{2}]} \pi(\frac{n}{p})-\pi(p)+1 $$
公式前面减一是因为 $1$ 是筛不掉的，其中 $p$ 为素数。

然后考虑如何求 dp 数组，考虑转移。考虑计算最后一个素数的贡献，发现其贡献是其倍数中没有被前面素数筛掉的数的个数，发现可以再利用 dp 数组，所以转移为 $dp_{i,j}=dp_{i,j-1}-dp_{\frac{i}{p_j},j-1}$，其中 $p_j$ 表示第 $j$ 个素数。

考虑边界情况，显然的是 $dp_{i,0}=i$。但这样还不够优，考虑 $p_j^2\ge i$ 的情况，这时 $[1,i]$ 中的合数全被筛过，只剩下素数，所以值即为 $\max(0,\pi(i)-j+1)$，加一是为了留下 $1$，与前面同步。

这样一个基础的 Meissel-Lehmer 算法原理讲解就完成了。

接下来是实现，对于 dp 数组，可以预处理出一部分，例如 $i\le 1.8 \times 10^6,j\le 60$ 的情况。然后由于这道题空间有限，前面公式里的 $\pi(\frac{n}{p})$ 会超出预处理的范围，这是要继续递归。

时间复杂度应该在 $O(n^\frac{2}{3})$ 级别，足以通过。

const int N = 8e6 + 10;
const int MI = 1.8e6;
const int MJ = 60;
ll n;
int f[MI + 5][MJ + 5];
int p[N / 10], ip[N], g[N], cnt;
void init() {
	FOR(i, 2, N - 1) {
		if(!ip[i]) {
			p[++cnt] = i;
		}
		for(int j = 1; j <= cnt && p[j] * i < N; j++) {
			ip[i * p[j]] = 1;
			if(i % p[j] == 0) {
				break;
			}
		}
	}
	FOR(i, 1, N - 1) g[i] = g[i - 1] + !ip[i];
	FOR(i, 1, MI) f[i][0] = i;
	FOR(i, 1, MI) FOR(j, 1, MJ)
		f[i][j] = f[i][j - 1] - f[i / p[j]][j - 1];
}
ll dp(ll i, ll j) {
	if(i <= MI && j <= MJ) return f[i][j];
	if(!i || !j) return i;
	if(i < N && 1ll * p[j] * p[j] >= i) return max(0ll, g[i] - j);
	return dp(i, j - 1) - dp(i / p[j], j - 1);
}
ll pi(ll n) {
	if(n < N) return g[n] - 1;
	ll sn = pow(n, 1. / 3);
	ll m = g[sn] - 1;
	ll res = dp(n, m) + m - 1;
	for(m++; 1ll * p[m] * p[m] <= n; m++)
		res -= pi(n / p[m]) - pi(p[m]) + 1;
	return res;
}

Min_25 筛

定义积性函数 $f(x)$，且 $f(p ^ k) = p ^ k(p ^ k - 1)$（$p$ 是一个质数），求 $$\sum_{i = 1} ^ n f(i)$$ 对 $10 ^ 9 + 7$ 取模。 $n \le 10^{10}$ 。

考虑拆解质因数，设 $n$ 的最小质因子为 $mp(n)$，

设 $S(n,j)$ = $\sum_{i=2}^{n}f(i)[mp(i)>j]$。

然后考虑转移，转移考虑质数和合数处的值。

设 $f(i)$ 质数处的前缀和为 $q(n)=\sum_{i=2}^{n}f(i)[i \in p]$。

那么质数处的和即为 $q(n)-\sum_{i=1}^{j}f(p_i)$。

然后对于合数部分，拆出最小质因数以及其次数，得出：

$\sum_{k>j}\sum_{c\le 1}f(p_k^c)(S(\frac{n}{p_k^c}, k)+[c>1])$

然后考虑 $q(n)$ 怎么求，为了删掉合数的贡献，考虑设 $g(n,j)=\sum_{i=2}^{n}f(i)[i \in p \mid mp(i)>j]$。

考虑 $g(n,j-1)$ 到 $g(n,j)$ 的转移。

当 $p_j^2>n$ 时，没有符合条件的合数，直接转移。

否则删去质因子含有 $p_j$ 的合数的贡献：

$g(n,j)=g(n,j-1)-f(p_j)(S(\frac{n}{p_j},j-1)-\sum_{i=1}^{j-1}f(p_i))$

不难发现这个转移中的 $f$ 只需要用到质数处的值，所以可以直接用一个完全积性函数替换掉。求 $g(n,0)$ 需要快速求出这个完全积性函数的前缀和，如果不好算可以进行多项式拆解，例如例题中的 $f(p)=p(p-1)=p^2-p$，分别计算 $f(p)=p^2$ 和 $f(p)=p$ 的值，然后减一减即可。

然后也不难发现 $g(i,j)$ 有用的 $i$ 就是所有 $\left \lfloor \frac{n}{k} \right \rfloor $，一共只有 $O(\sqrt{n})$ 种。

那么 $q(n)=g(n,|p|)$，答案等于 $S(n,0)+f(1)$。

时间复杂度 $O(\frac{n^\frac{3}{4}}{\log n})$。

const int N = 2e5 + 5;
const int P = 1e9 + 7;
inline int add(int x, int y) { return (x + y < P ? x + y : x + y - P); }
inline void Add(int &x, int y) { x = (x + y < P ? x + y : x + y - P); }
inline int sub(int x, int y) { return (x < y ? x - y + P : x - y); }
inline void Sub(int &x, int y) { x = (x < y ? x - y + P : x - y); }
inline int mul(int x, int y) { return (1ll * x * y) % P; }
inline void Mul(int &x, int y) { x = (1ll * x * y) % P; }
ll n; int sn, m;
int ip[N], p[N], cnt;
int s1[N], s2[N], g1[N], g2[N], id1[N], id2[N];
ll w[N];
void init() {
	ip[1] = 1;
	FOR(i, 2, sn) {
		if(!ip[i]) {
			p[++cnt] = i;
		}
		for(int j = 1; j <= cnt && p[j] * i <= sn; j++) {
			ip[p[j] * i] = 1;
			if(i % p[j] == 0) break;
		}
	}
	FOR(i, 1, cnt) {
		s1[i] = add(s1[i - 1], p[i]);
		s2[i] = add(s2[i - 1], mul(p[i], p[i]));
	}
}
int f1(int n) {
	return mul(mul(n, n + 1), (P + 1) / 2);
}
int f2(int n) {
	return mul(mul(mul(n, n + 1), n * 2 + 1), (P + 1) / 6);
}
int get(ll x) {
	if(x <= sn) return id1[x];
	return id2[n / x];
}
int S(ll n, int j) {
	if(p[j] > n) return 0;
	int res = sub(sub(g2[get(n)], g1[get(n)]), sub(s2[j], s1[j]));
	for(int i = j + 1; i <= cnt && 1ll * p[i] * p[i] <= n; i++) 
		for(ll c = 1, v = p[i]; v <= n; v *= p[i], c++)
			Add(res, mul(mul(v % P, v % P - 1), S(n / v, i) + (c > 1)));
	return res;
}
void solve() {
	cin >> n; sn = sqrt(n); init();
	for(ll l = 1, r; l <= n; l = r + 1) {
		r = n / (n / l);
		w[++m] = n / l;
		g1[m] = f1(w[m] % P) - 1;
		g2[m] = f2(w[m] % P) - 1;
		if(w[m] <= sn) id1[w[m]] = m;
		else id2[n / w[m]] = m;
	}
	FOR(i, 1, cnt) {
		for(int j = 1; j <= m && 1ll * p[i] * p[i] <= w[j]; j++) {
			Sub(g1[j], mul(p[i], sub(g1[get(w[j] / p[i])], s1[i - 1])));
			Sub(g2[j], mul(mul(p[i], p[i]), sub(g2[get(w[j] / p[i])], s2[i - 1])));
		}
	}
	cout << add(S(n, 0), 1) << endl;
}