杜教筛
/ / 阅读耗时 8 分钟 杜教筛可以在低于线性复杂度$(O(n^{\frac {2} {3}}))$下求积性函数前缀和。
先来看一个简单的问题:如何求莫比乌斯函数或者欧拉函数前缀和?
直接线性筛出两个函数,然后求前缀和即可,时间复杂度$O(n)$。这样做的效率确实很优秀,但是在ACM中就是会有人去卡这个复杂度,所以需要低于线性复杂度的算法,这就是杜教筛。
积性函数
对于一个数论函数$f(x)$,如果满足以下性质:
则称$f(x)$为积性函数。莫比乌斯函数和欧拉函数都是积性函数,另外约数个数$d(x)$以及约数和$\sigma(x)$也是积性函数。
积性函数与积性函数的乘积仍然是积性函数。
如果对于任意的整数$x$和$y$都有$f(xy)=f(x)f(y)$,则称$f(x)$为完全积性函数,它有以下几个:
- $\epsilon(x)$,元函数。在$x=1$是为1,否则为0,可以写成$\epsilon(x)=[x=1]$。
- $I(x)$,恒等函数。恒为1,$I(x)=1$。
- $id(x)$,单位函数。满足$id(x)=x$。
狄利克雷卷积
在之前的FFT/NTT以及FWT中介绍过很多卷积的快速求法,这里的狄利克雷卷积也是一种卷积,定义如下:
这里的$*$运算定义为$f$卷积$g$,它满足以下运算律:
- 交换律。即有$f*g=g*f$
- 结合律。即有$(f*g)*h=f*(g*h)$
- 对加法的分配律。即有$(f+g)*h=f*h+g*h$
狄利克雷卷积中存在幺元,即元函数$\epsilon$,满足$f*\epsilon=\epsilon*f=f$。
杜教筛
杜教筛巧妙利用积性函数以及狄利克雷卷积的性质将复杂的前缀和求解问题转化为简单问题,进而降低时间复杂度。
假如现在我们需要求积性函数$f(x)$的前缀和$\displaystyle\sum_{i=1}^nf(i)$。
构造两个积性函数满足$h=f*g$,那么有:
枚举因子:
记$S(x)=\displaystyle\sum_{i=1}^xf(i)$,那么有:
提出第一项:
移项:
化到这一步,可知如果$h$和$g$的前缀和易求,那么$S(x)$就可以很快地求出来,这就是杜教筛的思想。
模板
现在来用杜教筛求莫比乌斯函数的前缀和。注意到莫比乌斯函数的性质$\displaystyle\sum_{d|n}\mu(d)=[n=1]$,可知有$\mu*I=\epsilon$,那么直接代入上文推出的式子:
后者整除分块处理。
下面考虑求欧拉函数前缀和,注意到欧拉函数的一个性质$\displaystyle\sum_{d|n}\phi(d)=n$,即有$\phi*I=id$,那么代入公式,得到:
前者可以$O(1)$去求,后者整除分块。
最后来看如何求$x\phi(x)$的前缀和。可以将$g$设置成单位函数$id$,这样狄利克雷卷积就是:
所以$h(x)=x^2$,代入公式得到:
前者有公式:
然后后者整除分块即可。
那么如何代码实现呢?首先我们需要预处理$\sqrt n$范围的前缀和,这一步线性筛出,然后利用上面的公式递推,期间需要用hash表保存前缀和结果,可以用map也可以用平衡树。
附模板题代码:1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
using namespace std;
typedef unsigned long long ull;
unordered_map<int, int> mmp1;
unordered_map<int, ull> mmp2;
int n, prim[N + 5], tot = 0;
int ms[N + 5];
ull ps[N + 5];
bool vis[N + 5];
void init() {//预处理,线性筛
ms[1] = ps[1] = 1;
for (register int i = 2; i <= N; i++) {
if (!vis[i])ps[i] = i - 1, ms[i] = -1, prim[++tot] = i;
for (register int j = 1; j <= tot; j++) {
if (prim[j] * i > N)break;
vis[i * prim[j]] = true;
if (i % prim[j])ms[i * prim[j]] = -ms[i], ps[i * prim[j]] = ps[i] * (prim[j] - 1);
else {
ms[i * prim[j]] = 0, ps[i * prim[j]] = ps[i] * prim[j];
break;
}
}
}
for (int i = 2; i <= N; i++)ps[i] += ps[i - 1], ms[i] += ms[i - 1];//前缀和
}
int getS1(int x) {//莫比乌斯函数
if (x <= N)return ms[x];
if (mmp1.count(x) != 0)return mmp1[x];
register int ans = 0, l = 2, r;
while (l <= x) {
r = x / (x / l), ans += (r - l + 1) * getS1(x / l), l = r + 1;
if (l < 0)break;//防止溢出出错
}
return mmp1[x] = 1 - ans;
}
ull getS2(int x) {//欧拉函数
if (x <= N)return ps[x];
if (mmp2.count(x) != 0)return mmp2[x];
register ull ans = 0;
register int l = 2, r;
while (l <= x) {
r = x / (x / l), ans += (r - l + 1) * getS2(x / l), l = r + 1;
if (l < 0)break;
}
return mmp2[x] = x * (x + 1ll) / 2 - ans;
}
int main() {
init();
register int t;
cin >> t;
while (t--) {
cin >> n;
cout << getS2(n) << " " << getS1(n) << "\n";
}
return 0;
}