前言
这里假设你已经会了以下东西：
- 狄利克雷卷积
- 莫比乌斯反演
- 数论分块
正文
一切的开始
先来看一个题目（Bzoj 4805）：
给你一个整数$n$，求以下东西：
求$sum_{i=1}^n varphi_i$
可以线性筛吧...但是如果$nleq 2e9$呢？
杜教筛是啥？
设$mathbf S(n)=sum_{i=1}^nmathbf f(i)$，其中$mathbf f$是积性函数。
我们不妨设一个$mathbf g$，不管它具体是啥，做一个狄利克雷卷积：
$$ (mathbf gast mathbf f)(i)=sum_{d|i}mathbf g(d)mathbf f(frac id) $$
然后再求一个和：
$$ sum_{i=1}^nsum_{d|i}mathbf g(d)mathbf f(frac id) \=sum_{d=1}^nmathbf g(d)sum_{d|i}mathbf f(frac id) \=sum_{d=1}^nmathbf g(d)sum_{i=1}^{n/d}mathbf f(i) \=sum_{d=1}^nmathbf g(d)mathbf S(frac nd) $$
然后还有：
$$ mathbf g(1)mathbf S(n)=sum_{i=1}^nmathbf g(i)mathbf S(frac ni)- sum_{i=2}^nmathbf g(i)mathbf S(frac ni) $$
前面减数就是我们之前推的那个东西，后面被减数可以数论分块递归算：
$$ mathbf g(1)mathbf S(n)=sum_{i=1}^n(mathbf gast mathbf f)(i)- sum_{i=2}^nmathbf g(i)mathbf S(frac ni) $$
比如要求$sum_{i=1}^nmu(i)$
$$ mathbf g(1)mathbf S(n)=sum_{i=1}^n(mathbf gast mu)(i)- sum_{i=2}^nmathbf g(i)mathbf S(frac ni) $$
找一个积性函数让他们的狄利克雷卷积很好算：
我们知道$mu$的逆为$1$，也就是说$(1astμ)=epsilon$
那么单位元的前缀和是1，舒服啊，所以取$mathbf g(n)$为$mathbf 1(n)=1$，
所以原式可化为：
$$ mathbf S(n)=1-sum_{i=2}^nmathbf S(frac ni) $$
那么回归到例题
$$ mathbf S(n)=sum_{i=1}^n varphi_i \mathbf g(1)mathbf S(n)=sum_{i=1}^n(mathbf gast varphi)(i)- sum_{i=2}^nmathbf g(i)mathbf S(frac ni) $$
这时我们要找一个$mathbf g$使得$mathbf g$和$mathbf g ast varphi$的前缀和都很好算
这时我们想到了$varphi=mu*mathbf{id}$，所以啊：$varphi*mathbf 1=mathbf{id}$，所以同样取$mathbf g(n)$为$mathbf 1(n)=1$
$$ \mathbf S(n)=frac{n imes (n+1)}{2}- sum_{i=2}^nmathbf S(frac ni) $$
所以，杜教筛的重点就在于选取一个有利于计算的$mathbf g$...
下面给出这题的代码：
#include <cstdio>
#include <cstring>
#include <algorithm>
using std::min; using std::max;
using std::swap; using std::sort;
typedef long long ll;
const int N = 2147483647, M = 2e6 + 10, K = 1.3e3 + 10;
ll phi[M], s1[M], s2[K]; bool vis[K], notprime[M];
int _n, prime[M], tot;
template<typename T>
void read(T &x) {
int flag = 1; x = 0; char ch = getchar();
while(ch < '0' || ch > '9') { if(ch == '-') flag = -flag; ch = getchar(); }
while(ch >= '0' && ch <= '9') x = x * 10 + ch - '0', ch = getchar(); x *= flag;
}
ll S(int n) {
if(n < M) return s1[n];
int x = _n / n;
if(vis[x]) return s2[x];
vis[x] = true; ll &ans = s2[x];
ans = 1ll * n * (n + 1) / 2;
for(int i = 2, j; i <= n; i = j + 1)
j = n / (n / i), ans -= (j - i + 1) * S(n / i);
return ans;
}
void get_phi() {
phi[1] = 1;
for(int i = 2; i < M; ++i) {
if(!notprime[i]) prime[++tot] = i, phi[i] = i - 1;
for(int j = 1; j <= tot && i * prime[j] < M; ++j) {
notprime[i * prime[j]] = true;
if(!(i % prime[j])) { phi[i * prime[j]] = phi[i] * prime[j]; break; }
else phi[i * prime[j]] = phi[i] * (prime[j] - 1);
}
}
}
int main () {
read(_n), get_phi();
for(int i = 1; i < M; ++i) s1[i] = s1[i - 1] + phi[i];
printf("%lld
", S(_n));
return 0;
}
一些例题：
- Bzoj3512 DZY Loves Math IV
- LuoguP3768 简单的数学题