介绍
我们在 OI 题目中，不难会遇到给定一个积性函数 $f(x)$，求
$$\sum_{i=1}^N f(x)$$
的题目。虽然对于积性函数我们已经有成熟的线性筛方法来说可以做到 $O(n)$，可是有的出题人就是毒瘤，可是有的题目的数据范围要求给出一个复杂度小于线性的做法，这时候神仙就搞出了杜教筛这一东西。
前置芝士
芝士1：线性筛
首先我们知道，基本上所有的积性函数都可以线性筛。
积性函数的特点：
$$\forall x,y,(x,y) = 1 \Rightarrow f(xy) = f(x)f(y)$$
如果我们想线性筛积性函数，首先这个过程一定是建立在欧拉筛素数上的，然后我们只需要知道 $f(p^e)$ 的取值就行了。
我们拿 $\mu$ 举例：
什么你还不知道 $\mu$ 是什么？快去看看莫比乌斯反演吧。
首先发现 $\mu(p) = -1$，$p \in prime$。
然后发现当 $x$ 能不止被一种乘积筛出来的时候，$\mu(x) = 0$。
如果是被最小的质因数筛出的时候，就可以依据积性函数的定义来更新了。
直接放一个线性筛 $\mu$ 和 $\phi$ 的代码吧。
inline void pre(int n=MAXN-5){
phi[1] = mu[1] = 1;
FOR(i,2,n){
if(!p[i]){
prime[++cnt] = i;
mu[i] = -1;phi[i] = i-1;
}
for(Re int j = 1;j <= cnt && i*prime[j] <= n;++j){
p[i*prime[j]] = true;
if(i%prime[j]){
phi[i*prime[j]] = phi[i]*phi[prime[j]];
mu[i*prime[j]] = -mu[i];
}
else{
phi[i*prime[j]] = phi[i]*prime[j];
mu[i*prime[j]] = 0;
break;
}
}
}
}芝士2：常见积性函数
纯粹是列举一下防止后面因为符号不同的问题而看不懂。
$\mu(n)$，莫比乌斯函数。
$\phi(n)$，欧拉函数。
$d(n)$，约数个数，即 $d(n) = \sum_{i=1}^n [i|n]$。
$\sigma(n)$，约数和，即 $\sigma(n) = \sum_{i=1}^n i[i|n]$。
完全积性函数：$\forall x,y \Rightarrow f(xy) = f(x)f(y)$
$I(n)$：常函数 $I(n) = 1$。
$\epsilon(n)$：狄利克雷卷积单位元。
$id(n)$：单位函数 $id(i) = i$。
如果你学过莫比乌斯反演，应该已经掌握了所有的前置芝士了，现在我们进入正题。
杜教筛过程
下文中定义狄利克雷卷积的运算符号为 $ * $
我们现在要对于一个积性函数 $f(x)$，求
$$S(n) = \sum_{i=1}^n f(i)$$
考虑找到一个合适的函数 $g$，使得 $f * g$ 的前缀和易于计算。
我们考虑 $f * g$ 的前缀和，并加以推导，得
$$ \begin{aligned} &\sum_{i=1}^n (f * g)(i) \\ &= \sum_{i=1}^n \sum_{d|i} g(i)f(\frac{n}{i}) \\ &= \sum_{d=1}^n g(d) \sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} f(i) \\ &= \sum_{i=1}^n g(i) S(\lfloor \frac{n}{i}\rfloor) \\ &= g(1)S(n) + \sum_{i=2}^n g(i)S(\lfloor \frac{n}{i}\rfloor) \end{aligned} $$
考虑进行一些小的移项，最后可得我们一般使用的杜教筛的式子：
$$g(1)S(n) = \sum_{i=1}^n (f* g)(i) - \sum_{i=2}^n g(i)S(\lfloor \frac{n}{i} \rfloor)$$
如果我们对于 $(f* g)$ 的前缀和有快速求法（e.g. $O(1)$），那么预处理出 $f$ 函数的前 $n^{\frac{2}{3}}$，后面数论分块+递归求解，就可以做到 $O(n^{\frac{2}{3}})$ 的优秀复杂度内求得 $S(n)$。
现在我们来尝试用杜教筛来解决求一些函数的前缀和问题。
BZOJ 3944 Sum
题目链接
题目要让你求 $\sum_{i=1}^n \mu(i)$ 和 $\sum_{i=1}^n \phi(i)$，其中 $n \leq 2^{31}-1$。
首先考虑 $\mu$：
首先我们需要找到一个函数，使得它与 $\mu$ 的狄利克雷卷积的前缀和是容易计算的。注意到关于 $\mu$ 的一个性质：$\mu * I = \epsilon$。
发现 $\epsilon$ 的前缀和十分好算（就是 $1$），并且 $I$ 非常有利于乘法分块，于是我们取 $f(n) = \mu(n),g(n) = I(n)$，代入上式，得：
$$S(n) = 1-\sum_{i=2}^n S(\lfloor \frac{n}{i} \rfloor)$$
类比我们来考虑 $\phi$：
发现一个性质：$\phi * I = id$
$id$ 的前缀和可以用等差数列求和公式求得，代入上式可以发现：
$$S(n) = \sum_{i=1}^n i - \sum_{i=2}^n S(\lfloor \frac{n}{i} \rfloor)$$
也就是
$$S(n) = \frac{n(n+1)}{2} - \sum_{i=2}^n S(\lfloor \frac{n}{i} \rfloor)$$
现在我们来证明一下这东西的复杂度是什么 。
设 $T(n)$ 表示计算出 $S(n)$ 的复杂度，可以得到：
$$T(n) = O(\sqrt{n}) + \sum_{i=2}^{n}(T(i)+T(\frac{n}{i}))$$
展开一层就可以了，因为往下都是高阶小量可以不考虑。
$$T(n) =O(\sqrt{n})+\sum_{i=2}^n(O(\sqrt{i})+O(\sqrt{\lfloor \frac{n}{i} \rfloor})) = O(n^{\frac{3}{4}})$$
我们可以用筛法处理前 $k$ 个，如果 $k \geq \sqrt{n}$，则
$$T(n) = \sum_{i=1}^{\frac{n}{k}}O(\sqrt{\lfloor \frac{n}{i} \rfloor}) = O(\frac{n}{\sqrt{k}})$$
一般我们 $k$ 取 $n^{\frac{2}{3}}$，可以得到复杂度为 $O(n^{\frac{2}{3}})$ 的算法。
#include <tr1/unordered_map>
#include <algorithm>
#include <iostream>
#include <cstring>
#include <climits>
#include <cstdlib>
#include <cstdio>
#include <vector>
#include <cmath>
#include <ctime>
#include <queue>
#include <stack>
#include <map>
#include <set>
#define fi first
#define se second
#define U unsigned
#define P std::pair
#define Re register
#define LL long long
#define pb push_back
#define MP std::make_pair
#define all(x) x.begin(),x.end()
#define CLR(i,a) memset(i,a,sizeof(i))
#define FOR(i,a,b) for(Re int i = a;i <= b;++i)
#define ROF(i,a,b) for(Re int i = a;i >= b;--i)
#define DEBUG(x) std::cerr << #x << '=' << x << std::endl
using namespace std::tr1;
const int MAXN = 6000000+5;
bool p[MAXN];
int mu[MAXN],phi[MAXN],prime[10001],cnt;
LL sum1[MAXN];int sum2[MAXN];
inline void pre(int n=MAXN-5){
phi[1] = mu[1] = 1;
FOR(i,2,n){
if(!p[i]){
prime[++cnt] = i;
mu[i] = -1;phi[i] = i-1;
}
for(Re int j = 1;j <= cnt && i*prime[j] <= n;++j){
p[i*prime[j]] = true;
if(i%prime[j]){
phi[i*prime[j]] = phi[i]*phi[prime[j]];
mu[i*prime[j]] = -mu[i];
}
else{
phi[i*prime[j]] = phi[i]*prime[j];
mu[i*prime[j]] = 0;
break;
}
}
}
FOR(i,1,n) sum1[i] = sum1[i-1]+phi[i],sum2[i] = sum2[i-1]+mu[i];
}
unordered_map<LL,LL> S1;
unordered_map<int,int> S2;
inline LL calc1(LL x){
if(x <= MAXN-5) return sum1[x];
if(S1[x]) return S1[x];
LL ans = x*(x+1)/2;
for(Re LL l=2,r;l <= x;l = r+1){
r = x/(x/l);
ans -= calc1(x/l)*(r-l+1);
}
return S1[x] = ans;
}
inline int calc2(int x){
if(x <= MAXN-5) return sum2[x];
if(S2[x]) return S2[x];
int ans = 1;
for(Re int l = 2,r;l >= 0 && l <= x;l = r+1){
r = x/(x/l);
ans -= calc2(x/l)*(r-l+1);
}
return S2[x] = ans;
}
inline void Solve(){
int n;scanf("%d",&n);
printf("%lld %d\n",calc1(n),calc2(n));
}
int main(){
pre();
int T;scanf("%d",&T);
while(T--) Solve();
return 0;
}两道小题目
Sum 加强版
$$\sum_{i=1}^n (i \times \phi(i))$$
求上式的值，其中 $n$ 是线性筛跑不过的量级。
我们一眼看不出来 $g$ 用什么好，不妨先假定有了这个函数并代入狄利克雷卷积的形式：$\sum_{d|n} (d * \phi(d)) \times g(\frac{n}{d})$。
这个常数 d 太烦人了，不如我们想办法消去它。我们尝试取 $g = I$，代入化简得到
$$\sum_{d|n} (d* \phi(d))* \frac{n}{d} = n\sum_{d|n} \phi(d) = n^2$$
这样的话前缀和就很好求了，带回原式可得：
$$S(n) = \sum_{i=1}^n i^2 - \sum_{d=2}^n d* S(\lfloor \frac{n}{d}\rfloor)$$
大家都会 $O(1)$ 求 $\sum_{i=1}^n i^2$，也就是
$$\sum_{i=1}^n i^2 = \frac{n(n+1)(2n+1)}{6}$$
Hdu 5608 Function
题目链接
告诉你一个数论函数 $f(d)$，满足
$$\sum_{d|n} f(d) = n^2-3n+2$$
多次询问 $f$ 的前缀和。
$T \leq 500,N \leq 10^9$，只有五组 $N > 10^6$。
长得十分像莫比乌斯反演的经典形式。
我们设 $g(n) = n^2-3n+2$，原式转化为
$$g(n) = \sum_{d|n} f(d)$$
反演一波可以得到：
$$f(n) = \sum_{d|n} g(d)\mu(\frac{n}{d})$$
我们要求的答案变成了：
$$\sum_{i=1}^n \sum_{d|n} g(d) \mu(\frac{n}{d})$$
进行一些微小的变换，可以推导出
$$ \begin{aligned} && \sum_{i=1}^n \sum_{d|n} g(d) \mu(\frac{n}{d}) \\ && = \sum_{d=1}^n g(d) \sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} \mu(i) \end{aligned} $$
发现 $g$ 的前缀和十分好求，这样就可以通过数论分块+杜教筛 $\mu$ 来做了。
#include <unordered_map>
#include <algorithm>
#include <iostream>
#include <cstring>
#include <climits>
#include <cstdlib>
#include <cstdio>
#include <vector>
#include <cmath>
#include <ctime>
#include <queue>
#include <stack>
#include <map>
#include <set>
#define fi first
#define se second
#define U unsigned
#define P std::pair
#define Re register
#define LL long long
#define pb push_back
#define MP std::make_pair
#define all(x) x.begin(),x.end()
#define CLR(i,a) memset(i,a,sizeof(i))
#define FOR(i,a,b) for(Re int i = a;i <= b;++i)
#define ROF(i,a,b) for(Re int i = a;i >= b;--i)
#define DEBUG(x) std::cerr << #x << '=' << x << std::endl
const int MAXN = 2000000+5;
const int ha = 1e9 + 7;
bool p[MAXN];
int prime[MAXN],mu[MAXN],cnt;
inline int qpow(int a,int n=ha-2){
int res = 1;
while(n){
if(n & 1) res = 1ll*res*a%ha;
a = 1ll*a*a%ha;
n >>= 1;
}
return res;
}
inline void pre(int n=MAXN-5){
mu[1] = 1;p[1] = 1;
FOR(i,2,n){
if(!p[i]) prime[++cnt] = i,mu[i] = -1;
for(int j = 1;j <= cnt && i*prime[j] <= n;++j){
p[i*prime[j]] = true;
if(i%prime[j]) mu[i*prime[j]] = -mu[i];
else{
mu[i*prime[j]] = 0;
break;
}
}
}
FOR(i,1,n) mu[i] = (mu[i-1]+mu[i]+ha)%ha;
}
std::unordered_map<int,int> S;
inline int calc(int x){
if(x <= MAXN-10) return mu[x];
if(S[x]) return S[x];
int ans = 1;
for(int l = 2,r;l > 0 && l <= x;l = r+1){
r = x/(x/l);
ans = (ans-1ll*calc(x/l)*(r-l+1)%ha+ha)%ha;
}
return S[x] = ans;
}
int inv6 = qpow(6),inv2 = qpow(2);
inline int get(LL x){
int res = 1ll*x*(x+1)%ha*(2*x+1)%ha;
res = 1ll*res*inv6%ha;
res = (1ll*res+2ll*x)%ha;
int t = 3ll*x%ha*(x+1)%ha*inv2%ha;
res = (res-t+ha)%ha;
return res;
}
inline int get(int l,int r){
return (get(r)-get(l-1)+ha)%ha;
}
inline void Solve(){
int ans = 0,x;scanf("%d",&x);
for(int l = 1,r;l <= x;l = r+1){
r = x/(x/l);
int t = 1ll*calc(x/l)*get(l,r)%ha;
// DEBUG(l);DEBUG(r);DEBUG(t);
ans = ((ans+t)%ha+ha)%ha;
}
printf("%d\n",ans);
}
int main(){
pre();
int T;scanf("%d",&T);
while(T--) Solve();
return 0;
}
2 comments
%%% 杜教筛巨佬 RainAir
orz，RainAirtql