算法学习笔记（23）：杜教筛

杜教筛

参考来源： OI-Wiki， 网上博客
线性筛可以在线性时间求积性函数前缀和， 而杜教筛可以用低于线性时间求解积性函数前缀和。

我们考虑 \(S(n)\) 就是积性函数的前缀和， 所以我们尝试构造关于 \(\large S(n)\) 关于 \(\large S(\lfloor \frac{n}{i} \rfloor)\) 的递推式。

对于任意一个数论函数 g，必满足：

\(\begin{aligned} \sum_{i=1}^{n}(f * g)(i) & =\sum_{i=1}^{n}\sum_{d \mid i}g(d)f\left(\frac{i}{d}\right) \\ & =\sum_{i=1}^{n}g(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right) \end{aligned} \)

将 \(\large S(n)\) 提出来就可以得到：

\(\begin{aligned} g(1)S(n) & = \sum_{i=1}^n g(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right) - \sum_{i=2}^n g(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right) \\ & = \sum_{i=1}^n (f * g)(i) - \sum_{i=2}^n g(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right) \end{aligned} \)

再将 \(g(1)\) 移过去即可得到 \(S(n)\) 的递推式， 考虑我们可以自己构造 \(g(n)\) 使得计算变快。

当我们构造出这样的 \(g(n)\) 时：

1. 可以快速计算 \(\sum_{i=1}^n(f * g)(i)\)；
2. 可以快速计算 g 的前缀和，以用数论分块求解 \(\sum_{i=2}^ng(i)S\left(\left\lfloor\dfrac{n}{i}\right\rfloor\right)\)。

则我们可以在较短时间内求得 \(g(1)S(n)\)。

例题

P4213 【模板】杜教筛

求 \(\mu\) 和 \(\phi\) 的前缀和， 考虑一个性质， \(\mu * 1 = \varepsilon\)， \(\phi * 1 = id\)

考虑 \(\sum _{i = 1}^n(f * g)(i)\) 也就是 \(\sum id\) 和 \(\sum \varepsilon\)， 可以快速计算。

考虑 \(g(i)\) 也就是 \(1\)， 我们可以快速计算， 所以可以直接数论分块了。

点击查看代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned int uint;
const int N = 2e6 + 10;
unordered_map<int, ll> sphi, smu; 
ll phi[N], mu[N];
int pri[N], tot, n, T;
bool p[N];
void init() {
	phi[1] = 1; mu[1] = 1;
	for(int i = 2; i <= 2e6; i++) {
		if(!p[i]) pri[++tot] = i, phi[i] = i - 1, mu[i] = -1;
		for(int j = 1; j <= tot && i * pri[j] <= 2e6; j++) {
			p[i * pri[j]] = 1;
			if(i % pri[j] == 0) {
				phi[i * pri[j]] = phi[i] * pri[j];
				mu[i * pri[j]] = 0;
				break;
			} 
			else {
				phi[i * pri[j]] = phi[i] * (pri[j] - 1);
				mu[i * pri[j]] = -mu[i];
			}
		}
	}
	for(int i = 2; i <= 2e6; i++) phi[i] += phi[i - 1], mu[i] += mu[i - 1];
}
ll calcphi(int n) {
	if(n <= 2e6) return phi[n];
	if(sphi[n]) return sphi[n];
	ll res = (ll)n * ((ll)n + 1) / 2;
	for(uint l = 2, r; l <= n; l = r + 1) {
		r = n / (n / l);
		res -= 1ll * (r - l + 1) * calcphi(n / l);
	}	
	return sphi[n] = res;
}
ll calcmu(int n) {
	if(n <= 2e6) return mu[n];
	if(smu[n]) return smu[n];
	ll res = 1;
	for(uint l = 2, r; l <= n; l = r + 1) {
		r = n / (n / l);
		res -= 1ll * (r - l + 1) * calcmu(n / l);
	} 
	return smu[n] = res;
}
int main() {
	scanf("%d", &T); init();
	while(T--) {
		scanf("%d", &n);
		printf("%lld %lld\n", calcphi(n), calcmu(n));
	}
	return 0;
}