杜教筛学习笔记

用途&复杂度

杜教筛可用于计算一些积性函数的前缀和： $\sum_{i=1}^nf(i)$ 。

时间复杂度 $n^{2/3}$ ，需要预处理出 $n^{2/3}$ 数据量的前缀和。

原理

令 $S(n)=\sum_{i=1}^nf(i)$ 。

$\sum_{i=1}^n g(i)S(\lfloor\frac{n}{i}\rfloor)=\sum_{i=1}^n (f*g)(i)$

以上就是杜教筛最关键的点。

证明：

$\begin{align} \sum_{i=1}^n(f*g)(i) &= \sum_{i=1}^n\sum_{d|i}g(d)f(\lfloor\frac{i}{d}\rfloor)\\ &= \sum_{d=1}^n g(d)(\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}f(i))\\ &= \sum_{d=1}^ng(d)S(\lfloor\frac{n}{d}\rfloor)\\ \end{align}$

我们要求的是 $S(n)$ ，选择一个积性函数 $g(i)$ ，则 $g(1)=1$ 。

由简单容斥可得

$\begin{align} S(n) &= g(1)\cdot S(n)\\ &= \sum_{i=1}^n g(i)\cdot S(\lfloor\frac{n}{i}\rfloor)-\sum_{i=2}^n g(i)\cdot S(\lfloor\frac{n}{i}\rfloor)\\ \end{align}$

减号前面的部分由上面的关键点可得等于 $f$ 和 $g$ 的卷积的前缀和。

所以有

$\begin{align} S(n) &= \sum_{i=1}^n (f*g)(i)-\sum_{i=2}^n g(i)\cdot S(\lfloor\frac{n}{i}\rfloor)\\ \end{align}$

通过选择合适的积性函数 $g$ ，可以快速计算卷积的前缀和。而减号后面这部分可以通过整除分块计算， $\lfloor\frac{n}{i}\rfloor$ 的取值个数不超过 $\sqrt{n}$ 个。

例题

P4213 【模板】杜教筛（Sum）

求莫比乌斯函数以及欧拉函数的前缀和。

求莫比乌斯函数前缀和：

令 $g(i)=1$ 。

$\begin{align} S(n) &= \sum_{i=1}^n(\mu*I)(i)-\sum_{i=2}^n S(\lfloor\frac{n}{i}\rfloor)\\ &= \sum_{i=1}^n [i=1]-\sum_{i=2}^n S(\lfloor\frac{n}{i}\rfloor)\\ &= 1-\sum_{i=2}^n S(\lfloor\frac{n}{i}\rfloor)\\ \end{align}$

注意先预处理出 $n^{2/3}$ 个前缀和。

后面的部分分块求。

求欧拉函数前缀和：

还是令 $g(i)=1$ .

要知道 $I*\varphi = ID$ 。其中 $ID(x)=x$ 。

$\begin{align} S(n) &= \sum_{i=1}^n(\varphi*I)(i)-\sum_{i=2}^n S(\lfloor\frac{n}{i}\rfloor)\\ &= \sum_{i=1}^n ID(i)-\sum_{i=2}^n S(\lfloor\frac{n}{i}\rfloor)\\ &= \frac{(1+n)n}{2}-\sum_{i=2}^n S(\lfloor\frac{n}{i}\rfloor)\\ \end{align}$

也可以用莫比乌斯反演求。

$\begin{align} S'(n) &= \sum_{i=1}^n \sum_{j=1}^n [gcd(i,j)=1]\\ &= \sum_{i=1}^n\sum_{j=1}^n\sum_{d|gcd(i,j)}\mu(d)\\ &= \sum_{d=1}^n \mu(d)\cdot \lfloor\frac{n}{d}\rfloor^2\\ \end{align}$

由于 $S(n)=\sum_{i=1}^n\sum_{j=1}^i [gcd(i,j)=1]$ ，所以 $S(n)=\frac{S'(n)-1}{2}$ 。

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define debug(x) cout << #x << ":\t" << x << endl;
const int INF = 0x3f3f3f3f;
const ll inf = 0x3f3f3f3f3f3f3f3f;
const int N = 5e6 + 10;
const ll mod = 998244353;
int T;
ll n;
ll mu[N], phi[N], prime[N];
bool vis[N]; int cnt;
void pre(ll n) {
	mu[1] = phi[1] = 1;
	for (ll i = 2; i <= n; i++) {
		if (!vis[i]) {
			prime[cnt++] = i;
			mu[i] = -1;
			phi[i] = i - 1;
		}
		ll d;
		for (int j = 0; j < cnt && (d = i * prime[j]) <= n; j++) {
			vis[d] = 1;
			if (i%prime[j] == 0) {
				mu[d] = 0;
				phi[d] = phi[i] * prime[j];
				break;
			}
			else {
				mu[d] = -mu[i];
				phi[d] = phi[i] * phi[prime[j]];
			}
		}
	}
	for (int i = 1; i <= n; i++)mu[i] += mu[i - 1], phi[i] += phi[i - 1];
}
map<ll, ll>smu, sp;
ll S_mu(ll n) {
	if (n <= 5e6)return mu[n];
	if (smu.count(n))return smu[n];
	ll ans = 0;
	for (ll l = 2, r; l <= n; l = r + 1) {
		r = n / (n / l);
		ans += (r - l + 1)*S_mu(n / l);
	}
	return smu[n] = 1ll - ans;
}
ll S_p(ll n) {
	if (n <= 5e6)return phi[n];
	if (sp.count(n))return sp[n];
	ll ans = 0;
	for (ll l = 2, r; l <= n; l = r + 1) {
		r = n / (n / l);
		ans += (r - l + 1)*S_p(n / l);
	}
	return sp[n] = (1 + n)*n / 2 - ans;
}
int main() {
	pre(5e6);
	scanf("%d", &T);
	while (T--) {
		scanf("%lld", &n);
		printf("%lld %lld\n", S_p(n), S_mu(n));
	}
	return 0;
}