杜教筛 以及积性函数的前世今生 --算法竞赛专题解析（4）
《算法竞赛》2022.10月出版。网购：京东 天猫 当当
作者签名书（75折+快递）请加微信（手机号就是微信）：13916333036
文章目录
- 0 杜教筛简介
- 0.1 杜教筛的核心内容
- 0.2 杜教筛之起源
- 0.3 本文的结构
- 1 积性函数
- 1.1 积性函数定义
- 1.2积性函数的基本问题
- 1.3常见积性函数
- 2 欧拉函数
- 2.1 欧拉函数的定义和性质
- 2.2求欧拉函数的通解公式
- 2.3 线性筛（欧拉筛）求1～n内所有的欧拉函数
- 2.3.1 欧拉筛
- 2.3.2 求1～n内所有的欧拉函数
- 2.4 欧拉函数前缀和
- 3 整除分块（数论分块）
- 4 狄利克雷卷积
- 5 莫比乌斯函数和莫比乌斯反演
- 5.1 莫比乌斯函数的定义和性质
- 5.2 莫比乌斯函数的由来
- 5.3 莫比乌斯反演
- 6 杜教筛
- 6.1 杜教筛公式
- 6.1.1 杜教筛公式的推导
- 6.1.2 杜教筛算法和复杂度
- 6.2 欧拉函数前缀和
- 6.3 莫比乌斯函数前缀和
- 6.4 模板代码
- 6.5 题目
- 致谢
- 参考文献
本文详细讲解了与杜教筛有关的各种知识，包括积性函数、欧拉函数、整除分块、狄利克雷卷积、莫比乌斯反演、杜教筛等。以及它们的：理论知识、典型例题、关键代码、练习题。
本文的目标是：让一名 数论小白，也能最终看懂和应用杜教筛。
0 杜教筛简介
0.1 杜教筛的核心内容
杜教筛用途：在低于线性时间里，高效率求一些积性函数的前缀和。
杜教筛算法 = 整除分块 + 狄利克雷卷积 + 线性筛。
杜教筛公式：g(1)S(n)=∑i=1nh(i)−∑i=2ng(i)S(⌊ni⌋)g(1)S(n)=\sum_{i=1}^nh(i)-\sum_{i=2}^ng(i)S(\lfloor\frac{n}{i}\rfloor)g(1)S(n)=∑i=1nh(i)−∑i=2ng(i)S(⌊in⌋)
0.2 杜教筛之起源
“杜教筛”，一个亲切的算法名字，从2016年后流行于中国的OI圈。
杜教，即信息学竞赛的大佬杜瑜皓。本文作者好奇“杜教筛”这个名称的来源，在QQ上联系到他，他说，虽然算法的原理不是他的发明，但确实是他最早应用在OI竞赛中。至于为什么被称为“杜教筛”并广为传播，可能是一个迷了。
这个算法的原始出处，据skywalkert（唐靖哲，唐老师）的博客说，“这种黑科技大概起源于Project Euler这个网站，由xudyh（杜瑜皓）引入中国的OI、ACM界，目前出现了一些OI模拟题、OJ月赛题、ACM赛题是需要这种技巧在低于线性时间的复杂度下解决一类积性函数的前缀和问题。”
Project Euler网站也是一个OJ，题库里面有很多高难度的数学题。不过这个OJ并不需要提交代码，只需要提交答案即可。这个网站被推崇的一大原因是：如果提交的答案正确，就有权限看别人对这一题的解题方法，以及上传的代码。
经典的杜教筛题目，例如洛谷P4213，求数论函数的前缀和。
| 给定一个正整数n，n≤231−1n，n≤2^{31} −1n，n≤231−1，求: ans1=∑i=1nϕ(i)，ans2=∑i=1nμ(i)ans1=\sum_{i=1}^n\phi(i)，ans2=\sum_{i=1}^n\mu(i)ans1=i=1∑nϕ(i)，ans2=i=1∑nμ(i) |
|–|–|
题目中的ϕ(i)\phi(i)ϕ(i)是欧拉函数，μ(i)\mu(i)μ(i)是莫比乌斯函数。欧拉函数、莫比乌斯函数在算法竞赛中很常见，它们都是积性函数。题目求这2个积性函数的前缀和，由于给的n很大，显然不能直接用暴力的、复杂度O(n)O(n)O(n)的线性方法算，需要用杜教筛，复杂度是O(n23)O(n^{\frac23})O(n32)。
读者看到O(n23)O(n^{\frac23})O(n32)，可能感觉和O(n)O(n)O(n)相差不大，但实际上是一个很大的优化，两者相差n13n^{\frac13}n31倍。在上面的例题中，nnn=21亿，n23n^{\frac23}n32=160万，两者相差1300倍。
关于“杜教筛”的理论知识介绍，影响较大的有：
（1）2016年信息学奥林匹克中国国家队候选队员论文集，“积性函数求和的几种方法 任之洲”，其中提到杜教筛：“这个算法在国内信息学竞赛中一般称为杜教筛，可以高效地完成一些常见数论函数的计算。”
（2）2016年，skywalkert（唐老师）的博客。
0.3 本文的结构
由于杜教筛涉及的知识比较多，不容易讲清楚。所以本文将循循善诱地完成这一困难的任务，本文的顺序是：
（1）积性函数的定义和一些性质。积性函数的常见计算。
（2）以欧拉函数为例，讲解积性函数的计算。欧拉函数有什么用、如何求、它的求解和积性函数有什么关系、为什么用到杜教筛。
（3）整数分块。
（4）狄利克雷卷积。
（5）莫比乌斯函数和莫比乌斯反演。
（6）杜教筛。杜教筛公式、复杂度证明、欧拉函数和莫比乌斯函数的杜教筛实现。
1 积性函数
如果读者没有深入学过数论，第一次看到这些词语：积性函数、欧拉函数、狄利克雷卷积、莫比乌斯函数、莫比乌斯反演、杜教筛，是不是感到高深莫测、头皮发怵？
即使读者学过一些数论，但是还没有系统读过初等数论的书，那么在阅读本文之前，最好也读一本，比如《初等数论及其应用》Kenneth H.Rosen著，夏鸿刚译，机械工业出版社。这本书很易读，非常适合学计算机的人阅读：（1）概览并证明了初等数论的理论知识；（2）理论知识的各种应用，可以直接用在算法题目里；（3）大量例题和习题；（4）与计算机算法编程有很多结合。花两天时间通读很有益处。
本文所有的数学定义都引用自这本书。
1.1 积性函数定义
定义1:
（1）定义在所有正整数上的函数称为算术函数（或数论函数）。
（2）如果算术函数f对任意两个互质（互素）的正整数ppp和qqq，均有f(pq)=f(p)f(q)f(pq)=f(p)f(q)f(pq)=f(p)f(q)，称为积性函数（multiplicative functions，或译为乘性函数）2。
（3）如果对任意两个正整数ppp和qqq，均有f(pq)=f(p)f(q)f(pq)=f(p)f(q)f(pq)=f(p)f(q)，称为完全积性函数。
性质：
【定理1.1】 积性函数的和函数也是积性函数。如果fff是积性函数，那么fff的和函数F(n)=∑d∣nf(d)F(n)=\sum_{d|n}f(d)F(n)=∑d∣nf(d)也是积性函数。
其中d∣nd|nd∣n表示ddd整除nnn。
和函数在杜教筛中起到了关键的作用。
1.2积性函数的基本问题
（1）计算积性函数fff的第nnn项f(n)f(n)f(n)。
（2）计算积性函数fff在1～n1～n1～n范围内所有项：f(1)、f(2)、⋯、f(n)f(1)、f(2)、\cdots 、f(n)f(1)、f(2)、⋯、f(n)。
（3）计算积性函数fff前nnn项的和：∑i=1nf(i)\sum_{i=1}^nf(i)∑i=1nf(i)，即前缀和。这就是杜教筛的目标。
后面以欧拉函数和莫比乌斯函数为例，解答了这几个问题。
1.3常见积性函数
下面的积性函数，在狄利克雷卷积中常常用到。
I(n)I(n)I(n)：恒等函数，I(n)=1I(n)=1I(n)=1
id(n)id(n)id(n)：单位函数，id(n)=nid(n)=nid(n)=n
Ik(n)I_k(n)Ik(n)：幂函数，Ik(n)=nkI_k(n)=n^kIk(n)=nk
ϵ(n)\epsilon(n)ϵ(n)：元函数，ϵ(n)={1,n=10,n>1\epsilon(n)= \begin{cases} 1, & \text {n=1} \\ 0, & \text{n>1} \end{cases}ϵ(n)={1,0,n=1n>1
σ(n)\sigma(n)σ(n)：因子和函数，σ(n)=∑d∣nd\sigma(n)=\sum_{d|n}dσ(n)=∑d∣nd
d(n)d(n)d(n)：约数个数，d(n)=∑d∣n1d(n)=\sum_{d|n}1d(n)=∑d∣n1
2 欧拉函数
为了彻底了解积性函数的运算，本节以欧拉函数为例进行讲解，这是一个经典的积性函数。并用这个例子引导出为什么要使用杜教筛。
2.1 欧拉函数的定义和性质
定义：设n是一个正整数，欧拉函数ϕ(n)\phi(n)ϕ(n)定义为不超过nnn且与nnn互质的正整数的个数。
定义用公式表示：ϕ(n)=∑i=1n[gcd(i,n)=1]\phi(n)=\sum_{i=1}^n[gcd(i,n)=1]ϕ(n)=∑i=1n[gcd(i,n)=1]
例如：
n = 4时，1、3这2个数与4互质，ϕ(4)=2\phi(4)=2ϕ(4)=2；
n = 9时，1、2、4、5、7、8这6个数与9互质，ϕ(9)=6\phi(9)=6ϕ(9)=6。
【定理2.1】3 设p和q是互质的正整数，那么ϕ(pq)=ϕ(p)ϕ(q)\phi(pq)=\phi(p)\phi(q)ϕ(pq)=ϕ(p)ϕ(q)。
这个定理说明欧拉函数是一个积性函数。有以下推理：
若n=p1a1×p2a2×⋯×psasn=p_1^{a_1}\times p_2^{a_2}\times \cdots\times p_s^{a_s}n=p1a1×p2a2×⋯×psas，其中p1、p2、⋯、psp_1、p_2、\cdots、p_sp1、p2、⋯、ps互质，a1、a2、⋯、asa_1、a_2、\cdots、a_sa1、a2、⋯、as是它们的幂，则ϕ(n)=ϕ(p1a1)×ϕ(p2a2)×⋯×ϕ(psas))\phi(n)=\phi(p_1^{a_1}) \times \phi(p_2^{a_2}) \times \cdots\times\phi(p_s^{a_s}) )ϕ(n)=ϕ(p1a1)×ϕ(p2a2)×⋯×ϕ(psas)) 。
如果p1、p2、...、psp_1、p_2、...、p_sp1、p2、...、ps互质，那么p1a1、p2a2、⋯、psasp_1^{a_1}、p_2^{a_2}、\cdots、p_s^{a_s}p1a1、p2a2、⋯、psas也是互质的。
以n=84n = 84n=84为例：
84=22×3×784=2^2 \times 3 \times 784=22×3×7
ϕ(84)=ϕ(22×3×7)=ϕ(22)×ϕ(3)×ϕ(7))\phi(84)=\phi(2^2 \times 3 \times 7)=\phi(2^2) \times \phi(3) \times \phi(7))ϕ(84)=ϕ(22×3×7)=ϕ(22)×ϕ(3)×ϕ(7))
【定理2.2】 设nnn为正整数，那么n=∑d∣nϕ(d)n=\sum_{d|n}\phi(d)n=∑d∣nϕ(d)
其中d∣nd|nd∣n表示ddd整除nnn，上式对nnn的正因子求和。 例如n=12n=12n=12，那么ddd是1、2、3、4、6、121、2、3、4、6、121、2、3、4、6、12，有：12=ϕ(1)+ϕ(2)+ϕ(3)+ϕ(4)+ϕ(6)+ϕ(12)12=\phi(1)+\phi(2)+\phi(3)+\phi(4)+\phi(6)+\phi(12)12=ϕ(1)+ϕ(2)+ϕ(3)+ϕ(4)+ϕ(6)+ϕ(12)
这个定理说明了n与的ϕ(n)\phi(n)ϕ(n)关系：nnn的正因子（包括111和自身）的欧拉函数之和等于nnn。
有很多方法可以证明，这里给出一种直接易懂的证明。
证明：分数1n,2n,⋯,nn\frac 1n,\frac 2n,\cdots,\frac nnn1,n2,⋯,nn，共有nnn个，且互不相等。化简这些分数，得到新的nnn个分数，它们的分母和分子互质，形如ad\frac adda，d∣nd|nd∣n且aaa与ddd互质。在所有nnn个分数中，分母为ddd的分数，数量是ϕ(d)\phi(d)ϕ(d)。所有不同分母的分数，其总数是∑d∣nϕ(d)\sum_{d|n}\phi(d)∑d∣nϕ(d)，所以n=∑d∣nϕ(d)n=\sum_{d|n}\phi(d)n=∑d∣nϕ(d)，得证。
定理2.2是欧拉函数特别重要的一个性质，这一性质被用到了杜教筛中。
2.2求欧拉函数的通解公式
欧拉函数有什么用呢？利用它可以推出欧拉定理4，欧拉定理可用于：求大数的模、求解线性同余方程等，在密码学中有重要应用。所以求ϕ(n)\phi(n)ϕ(n)是一个常见的计算。
从定理2.1可以推出定理2.3，定理2.3是计算ϕ(n)\phi(n)ϕ(n)的通解公式。
【定理2.3】 设：n=p1a1×p2a2×⋯×psasn=p_1^{a_1}\times p_2^{a_2}\times \cdots\times p_s^{a_s}n=p1a1×p2a2×⋯×psas为正整数nnn的质幂因子分解，那么：
ϕ(n)=n(1−1p1)(1−1p2)⋯(1−1pk)=n∏i=1k(1−1pi)\phi(n)=n(1-\frac 1{p_1})(1-\frac 1{p_2})\cdots(1-\frac 1{p_k})=n \prod_{i=1}^k(1-\frac 1{p_i})ϕ(n)=n(1−p11)(1−p21)⋯(1−pk1)=ni=1∏k(1−pi1)
上述公式，有两个特殊情况：
（1）若nnn是质数，ϕ(n)=n−1\phi(n)=n-1ϕ(n)=n−1。这一点很容易理解，请读者思考。
（2）若n=pkn=p^kn=pk，ppp是质数，有：ϕ(n)=ϕ(pk)=pk−pk−1=pk−1(p−1)=pk−1ϕ(p)\phi(n)=\phi(p^k)=p^k-p^{k-1}=p^{k-1}(p-1)=p^{k-1}\phi(p)ϕ(n)=ϕ(pk)=pk−pk−1=pk−1(p−1)=pk−1ϕ(p)
这两种情况将用于2.3.2节编程求欧拉函数。
所以，欧拉函数ϕ(n)\phi(n)ϕ(n)的求解，归结到了分解质因子这个问题。分解质因子有一个古老的方法，试除法5：求nnn的质因子时，逐个检查从222到n\sqrt{n}n的所有质数，如果它能整除nnn，就是一个因子。试除法的复杂度是O(n)O(\sqrt{n})O(n)，效率很低。在密码学中，有高达百位以上的十进制数分解质因子，因此发明了很多高效率的方法。不过，在算法竞赛中，数据规模不大，所以一般就用试除法。
下面是求欧拉函数的代码。先用试除法分解质因子，再用公式ϕ(n)=n∏i=1k(1−1pi)\phi(n)=n \prod_{i=1}^k(1-\frac 1{p_i})ϕ(n)=n∏i=1k(1−pi1)求得ϕ(n)\phi(n)ϕ(n)。总复杂度仍然是O(n)O(\sqrt{n})O(n)。
求单个欧拉函数的代码
int euler(int n){int ans = n;for(int p = 2; p*p <= n; ++ p){ //试除法：检查从2到sqrt(n)的每个数if(n%p == 0){ //能整除，p是一个因子，而且是质因子，请思考ans = ans/p*(p-1); //求欧拉函数的通式while(n%p == 0) //去掉这个因子的幂，并使得下一个p是质因子n /= p; //减小了n}}if(n != 1) //情况(1)：n是一个质数，没有执行上面的分解ans = ans/n*(n-1);return ans;
}
2.3 线性筛（欧拉筛）求1～n内所有的欧拉函数
现在要求1～n1～n1～n范围内所有的欧拉函数，而不是只求一个欧拉函数。前面求一个欧拉函数的复杂度是O(n)O(\sqrt{n})O(n)，如果一个一个地求1～n1～n1～n范围内所有的欧拉函数，那么nnn个的总复杂度是O(nn)O(n\sqrt{n})O(nn)。效率太低。
这个过程可以优化：
（1）利用积性函数的性质ϕ(p1p2)=ϕ(p1)ϕ(p2)\phi(p_1p_2)=\phi(p_1)\phi(p_2)ϕ(p1p2)=ϕ(p1)ϕ(p2)，求得前面的ϕ(p1)\phi(p_1)ϕ(p1)和ϕ(p2)\phi(p_2)ϕ(p2)后，可以递推出后面的ϕ(p1p2)\phi(p_1p_2)ϕ(p1p2)；
（2）求1～n1～n1～n范围内的质数，用于上一步骤的递推。
编程时，可以这样操作：
（1）求得1～n1～n1～n内每个数的最小质因子。所谓最小质因子，例如84=22×3×784=2^2 \times 3 \times 784=22×3×7，848484的最小质因子是222。这步操作可以用欧拉筛，复杂度是O(n)O(n)O(n)的，不可能更快了。
（2）在第（1）步基础上，递推求1～n1～n1～n内每个数的欧拉函数。这一步的复杂度也是O(n)O(n)O(n)的。
两者相加，总复杂度是O(n)O(n)O(n)的。
2.3.1 欧拉筛
欧拉筛是一种线性筛，也就是说，它能在O(n)O(n)O(n)的线性时间里求得1～n1～n1～n内所有的质数。
欧拉筛是对埃氏筛法的改进。埃氏筛法的原理是：每个合数都是多个质数的积；那么从最小的质数2开始，用每一个质数去筛比它大的数，就能筛掉合数。埃氏筛法低效的原因是，一个合数会被它的多个质因子重复筛。请读者自己详细了解埃氏筛法。
欧拉筛法的原理是：一个合数肯定有一个最小质因子；让每个合数只被它的最小质因子筛选一次，以达到不重复筛的目的。
具体操作步骤：
（1）逐一检查2～n2～n2～n的所有数。第一个检查的是222，它是第一个质数。
（2）当检查到第iii个数时，利用已经求得的质数去筛掉对应的合数xxx，而且是用xxx的最小质因子去筛。
下面是代码。代码很短很精妙。
欧拉筛求素数
int prime[MAXN]; //保存质数
int vis[MAXN]; //记录是否被筛
int euler_sieve(int n){ //欧拉筛。返回质数的个数。int cnt = 0; //记录质数个数memset(vis,0,sizeof(vis));memset(prime,0,sizeof(prime));for(int i=2;i<=n;i++){ //检查每个数，筛去其中的合数if(!vis[i]) prime[cnt++]=i; //如果没有筛过，是质数，记录。第一个质数是2for(int j=0; j<cnt; j++){ //用已经得到的质数去筛后面的数if(i*prime[j] >n) break; //只筛小于等于n的数vis[i*prime[j]]=1; //关键1。用x的最小质因数筛去xif(i%prime[j]==0) break; //关键2。如果不是这个数的最小质因子，结束}}return cnt; //返回小于等于n的素数的个数
}
代码中难懂的是后面2个关键行。在“关键1”， 其中的prime[j]prime[j]prime[j]是最小质因子，vis[i∗prime[j]]=1;vis[i*prime[j]]=1;vis[i∗prime[j]]=1;表示用最小质因子筛去了它的倍数。在“关键2”，及时跳出，避免重复筛。
下面用图解来说明，请读者对照代码和图解自己理解。特别注意图4，它是“关键2”。
图1. 初始化
图2. i=2，用质数2筛去4（2是4的最小质因子）
图3. i=3，用质数2筛去6，用质数3筛去9
图4. i=4，用质数2筛去8，注意质数3没有用到，循环j被跳出了
图5. i=5，用质数2筛去10，用3筛去15，用5筛去25
可以发现，每个数xxx只被筛了一次，而且是被xxx的最小质因子筛去的。这说明筛法是线性的，复杂度 O(n)O(n)O(n)。
现在用欧拉筛求1～n1～n1～n内每个数的最小质因子。只要简单修改上述程序，直接用vis[MAXN]vis[MAXN]vis[MAXN]记录最小质因子就行。代码修改如下，修改了带注释的后2行。
欧拉筛求质数+最小质因子
int prime[MAXN]; //记录质数
int vis[MAXN]; //记录最小质因子
int euler_sieve(int n){ int cnt=0;memset(vis,0,sizeof(vis));memset(prime,0,sizeof(prime));for(int i=2;i<=n;i++){ if(!vis[i]){ vis[i]=i; prime[cnt++]=i;} //vis[]记录x的最小质因子for(int j=0; j<cnt; j++){ if(i*prime[j] >n) break; vis[i*prime[j]]=prime[j]; //vis[]记录x的最小质因子if(i%prime[j]==0)break; }}return cnt;
}
2.3.2 求1～n内所有的欧拉函数
下面用一个模板题给出模板代码，它就是“欧拉筛+欧拉函数公式”。
例题：洛谷P2158 仪仗队 https://www.luogu.com.cn/problem/P2158
这是欧拉函数的模板题。get_phi()是计算欧拉函数的模板，其中判断了3种情况，细节见注释。
情况（1）：若nnn是质数，ϕ(n)=n−1\phi(n)=n-1ϕ(n)=n−1。
情况（2）：若n=pkn=p^kn=pk，ppp是质数，有ϕ(n)=pk−1ϕ(p)\phi(n)=p^{k-1}\phi(p)ϕ(n)=pk−1ϕ(p)：
情况（3）：若n=p1p2n=p_1p_2n=p1p2，用ϕ(n)=ϕ(p1)ϕ(p2)\phi(n)=\phi(p_1)\phi(p_2)ϕ(n)=ϕ(p1)ϕ(p2)递推。
洛谷P2158代码
#include<bits/stdc++.h>
using namespace std;const int maxn = 50000;
int vis[maxn]; //记录是否被筛；或者用于记录最小质因子
int prime[maxn]; //记录质数
int phi[maxn]; //记录欧拉函数
int sum[maxn]; //计算欧拉函数的和void get_phi(){ //模板：求1～maxn范围内的欧拉函数phi[1]=1;int cnt=0;for(int i=2;i<maxn;i++) {if(!vis[i]) {vis[i]=i; //vis[i]=1; //二选一：前者记录最小质因子，后者记录是否被筛prime[cnt++]=i; //记录质数phi[i]=i-1; //情况（1）：i是质数，它欧拉函数值=i-1}for(int j=0;j<cnt;j++) {if(i*prime[j] > maxn) break;vis[i*prime[j]] = prime[j]; //vis[i*prime[j]]=1; //二选一：前者记录最小质因子，后者记录是否被筛if(i%prime[j]==0){ //prime[j]是最小质因子phi[i*prime[j]] = phi[i]*prime[j];//情况（2）：i是prime[j]的k次方break;}phi[i*prime[j]] = phi[i]*phi[prime[j]]; //情况（3）：i和prime[j]互质，递推出i*prime[j]}}
}int main(){get_phi(); //计算所有的欧拉函数sum[1]=1;for(int i=2;i<=maxn;i++) //打表计算欧拉函数的和sum[i]=sum[i-1]+phi[i];int n; scanf("%d",&n);if(n==1) printf("0\n");else printf("%d\n",2*sum[n-1]+1);return 0;
}
本节的例子是用欧拉筛求欧拉函数，读者可以认识到：积性函数都可以用欧拉筛来求。后面还有一个例子：用欧拉筛求解积性函数莫比乌斯函数。
2.4 欧拉函数前缀和
在欧拉函数问题的最后，回到杜教筛的模板题，求：∑i=1nϕ(i)\sum_{i=1}^n\phi(i)∑i=1nϕ(i)。
简单的方法，是用2.3.2节计算出每个欧拉函数，再求和，复杂度是O(n)O(n)O(n)。
更快的方法，就是本篇的主题“杜教筛”，复杂度是O(n23)O(n^{\frac23})O(n32)。
“杜教筛”需要整除分块、狄利克雷卷积等前置知识。
3 整除分块（数论分块）
整除分块是杜教筛算法的基本思路。
整除分块是为了解决一个整除的求和问题：∑i=1n⌊ni⌋\sum_{i=1}^n\lfloor\dfrac{n}{i}\rfloor∑i=1n⌊in⌋
其中⌊ni⌋\lfloor\dfrac{n}{i}\rfloor⌊in⌋表示nnn除以iii，并向下取整6。nnn是给定的一个整数。
如果直接暴力计算，复杂度是O(n)O(n)O(n)的，若nnn很大会超时。但是，用整除分块算法来求解，复杂度是O(n)O(\sqrt{n})O(n)的。
下面以n=8n = 8n=8为例，列出每个情况：
| i | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
|---|---|---|---|---|---|---|---|---|
| 8i\frac 8ii8 | 8 | 4 | 2 | 2 | 1 | 1 | 1 | 1 |
对第2行求和，结果是：∑i=18⌊8i⌋=20\sum_{i=1}^8\lfloor\frac{8}{i}\rfloor=20∑i=18⌊i8⌋=20
观察上面第222行8i\dfrac 8ii8的值，发现是逐步变小的，并且很多相等，这是整除操作的规律。例如：
8/1=88/1 = 88/1=8
8/2=48/2 = 48/2=4
8/3=8/4=28/3 = 8/4 = 28/3=8/4=2
8/5=8/6=8/7=8/8=18/5 = 8/6 = 8/7 = 8/8 = 18/5=8/6=8/7=8/8=1
为了对这些整除的结果进行快速求和，自然能想到，把它们分成一块块的，其中每一块的8i\dfrac 8ii8相同，把这一块一起算，就快多了。
设每一块的左右区间是[l,r][l, r][l,r]，上面的例子可以分成4块：[1,1]、[2,2]、[3,4]、[5,8][1,1]、[2,2]、[3,4]、[5,8][1,1]、[2,2]、[3,4]、[5,8]。
每一块内部的求和计算，只要用数量乘以值就行了，例如[5,8][5,8][5,8]区间的整除和：(8−5+1)∗1=4(8-5+1)*1=4(8−5+1)∗1=4。这个计算的复杂度是O(1)O(1)O(1)的。
最后还有两个问题：
（1）把⌊ni⌋\lfloor\dfrac{n}{i}\rfloor⌊in⌋按相同值分块，一共需要分成几块？或者说，⌊ni⌋\lfloor\dfrac{n}{i}\rfloor⌊in⌋有几种取值？这决定了算法的复杂度。
【答】分块少于2n2\sqrt{n}2n种，证明如下：
（a）i≤ni\leq\sqrt{n}i≤n时，ni\dfrac{n}{i}in的值有：{n1,n2,n3⋯nn}\lbrace \frac{n}{1},\frac{n}{2},\frac{n}{3}\cdots \frac{n}{\sqrt{n}} \rbrace{1n,2n,3n⋯nn}，ni≥n\frac{n}{i}\geq\sqrt{n}in≥n，共n\sqrt{n}n个，此时⌊ni⌋\lfloor\dfrac{n}{i}\rfloor⌊in⌋有n\sqrt{n}n种取值；
（b）i>ni>\sqrt{n}i>n时，有ni<n\dfrac{n}{i}<\sqrt{n}in<n，此时⌊ni⌋\lfloor\dfrac{n}{i}\rfloor⌊in⌋也有n\sqrt{n}n种取值。
两者相加，共2n2\sqrt{n}2n种。所以，整除分块的数量是O(n)O(\sqrt{n})O(n)的。
（2）如何计算？或者说，给定每个块的第一个数lll，能推出这个块的最后一个数rrr吗？
【答】可以证明：r=n/(n/l)r = n/(n/l)r=n/(n/l)。证明略7。
下面是标程。
#include<bits/stdc++.h>
using namespace std;
int main(){long long n,ans=0;cin >> n;for(long long l=1,r;l<=n;l=r+1){r = n/(n/l); //计算r，让分块右移ans += (r-l+1)*(n/l); //求和cout << l <<""<< r <<": "<< n/r << endl; //打印分块}cout << ans; //打印和
}
整除分块习题：洛谷P1829、洛谷P2261、洛谷P3935。
4 狄利克雷卷积
狄利克雷卷积 Dirichlet Convolution，是计算求和问题的有用工具。狄利克雷卷积是杜教筛用到的核心技术之一。
定义：
设f，gf，gf，g是算术函数，记fff和ggg的狄利克雷卷积是f∗gf*gf∗g，定义为8：
(f∗g)(n)=∑d∣nf(d)g(nd)(f*g)(n)= \sum_{d|n}f(d)g(\frac{n}{d}) (f∗g)(n)=d∣n∑f(d)g(dn)
它是一个求和公式，其中d∣nd|nd∣n表示ddd整除nnn，卷积是对正因子求和。
这个卷积是什么含义呢？现在给出一个特殊的例子。定义恒等函数I(n)=nI(n)=nI(n)=n、常数函数1(n)=11(n)=11(n)=1，它们的卷积是：
(I∗1)(n)=∑d∣nI(d)1(nd)=∑d∣nd⋅1=∑d∣nd=σ(n)(I*1)(n)=\sum_{d|n}I(d)1(\frac{n}{d})=\sum_{d|n}{d} \cdot1=\sum_{d|n}d=\sigma(n)(I∗1)(n)=d∣n∑I(d)1(dn)=d∣n∑d⋅1=d∣n∑d=σ(n)
把它记为：I∗1=σI*1=\sigmaI∗1=σ
σ(n)\sigma(n)σ(n)是“因子和函数”的符号。例如n=12，那么σ(n)=1+2+3+4+6+12=28\sigma(n)=1+2+3+4+6+12=28σ(n)=1+2+3+4+6+12=28。σ(n)\sigma(n)σ(n)是积性函数。
狄利克雷卷积的性质：
（1）交换律：f∗g=g∗ff*g=g*ff∗g=g∗f
（2）结合律：(f∗g)∗h=f∗(g∗h)(f*g)*h=f*(g*h)(f∗g)∗h=f∗(g∗h)
（3）分配律：f∗(g+h)=(f∗g)+(f∗h)f*(g+h)=(f*g)+(f*h)f∗(g+h)=(f∗g)+(f∗h)
（4）元函数：ϵ(n)={1,n=10,n>1\epsilon(n)= \begin{cases} 1, & \text {n=1} \\ 0, & \text{n>1} \end{cases}ϵ(n)={1,0,n=1n>1。对任何fff，有e∗f=f∗e=fe*f=f*e=fe∗f=f∗e=f
（5）两个积性函数的狄利克雷卷积仍然是积性函数。
有一些积性函数算术函数在狄利克雷卷积中常常用到，1.3节已经给出，这里再次列出：
I(n)I(n)I(n)：恒等函数，I(n)=1I(n)=1I(n)=1
id(n)id(n)id(n)：单位函数，id(n)=nid(n)=nid(n)=n
Ik(n)I_k(n)Ik(n)：幂函数，Ik(n)=nkI_k(n)=n^kIk(n)=nk
ϵ(n)\epsilon(n)ϵ(n)：元函数，ϵ(n)={1,n=10,otherwise\epsilon(n)= \begin{cases} 1, & \text {n=1} \\ 0, & \text{otherwise} \end{cases}ϵ(n)={1,0,n=1otherwise
σ(n)\sigma(n)σ(n)：因子和函数，σ(n)=∑d∣nd\sigma(n)=\sum_{d|n}dσ(n)=∑d∣nd
d(n)d(n)d(n)：约数个数，d(n)=∑d∣n1d(n)=\sum_{d|n}1d(n)=∑d∣n1
5 莫比乌斯函数和莫比乌斯反演
如果读者困惑莫比乌斯函数有什么用，可以先看下一节“5.2 莫比乌斯函数的由来”。
5.1 莫比乌斯函数的定义和性质
莫比乌斯函数以数学家Möbius的名字命名9。
莫比乌斯函数是狄利克雷卷积的重要工具，它也是数论和组合数学的重要的积性函数。
莫比乌斯函数μ(n)\mu(n)μ(n)定义为10：
μ(n)={1,如果n=1(−1)r,如果n=p1p2...pr，其中pi是不同的质数0,其他情况\mu(n)= \begin{cases} 1, & \text {如果n=1}\\ (-1)^r, & \text{如果$n = p_1 p_2...p_r，其中p_i是不同的质数$} \\ 0, & \text{其他情况} \end{cases}μ(n)=⎩⎨⎧1,(−1)r,0,如果n=1如果n=p1p2...pr，其中pi是不同的质数其他情况
一些例子：μ(1)=1,μ(2)=−1,μ(3)=−1,μ(4)=μ(22)=0,μ(5)=−1,μ(6)=μ(2×3)=1,μ(7)=−1,μ(8)=μ(23)=0;μ(330)=μ(2×3×5×11)=(−1)4=1,μ(660)=μ(22×3×5×11)=0。\mu(1)=1,\mu(2)=-1,\mu(3)=-1,\mu(4)=\mu(2^2)=0,\mu(5)=-1,\mu(6)=\mu(2\times3)=1,\mu(7)=-1,\mu(8)=\mu(2^3)=0; \mu(330)=\mu(2\times3\times5\times11)=(-1)^4=1,\mu(660)=\mu(2^2\times3\times5\times11)=0。μ(1)=1,μ(2)=−1,μ(3)=−1,μ(4)=μ(22)=0,μ(5)=−1,μ(6)=μ(2×3)=1,μ(7)=−1,μ(8)=μ(23)=0;μ(330)=μ(2×3×5×11)=(−1)4=1,μ(660)=μ(22×3×5×11)=0。
莫比乌斯函数μ(n)\mu(n)μ(n)是积性函数。它有一个和欧拉函数的n=∑d∣nϕ(d)n=\sum_{d|n}\phi(d)n=∑d∣nϕ(d)类似的定理：
【定理5.1】莫比乌斯函数的和函数在整数nnn处的值，满足：
∑d∣nμ(d)={1,n=10,n>1\sum_{d|n}\mu(d)= \begin{cases} 1, & \text {n=1} \\ 0, & \text{n>1} \end{cases}d∣n∑μ(d)={1,0,n=1n>1
证明：
（1）n=1n=1n=1时，显然有F(1)=∑d∣nμ(1)=1F(1)=\sum_{d|n}\mu(1)=1F(1)=∑d∣nμ(1)=1。
（2）n>1n>1n>1时。根据积性函数的定义，有F(n)=F(P1a1)F(P2a2)⋯F(Ptat)F(n)=F(P_1^{a_1})F(P_2^{a_2})\cdots F(P_t^{a_t})F(n)=F(P1a1)F(P2a2)⋯F(Ptat)，其中n=P1a1P2a2⋯Ptatn=P_1^{a_1}P_2^{a_2} \cdots P_t^{a_t}n=P1a1P2a2⋯Ptat是质因子分解。如果能证明F(Pk)=0F(P^k)=0F(Pk)=0即有F(n)=0F(n)=0F(n)=0。因为当i≥2i≥2i≥2时，μ(pi)=0\mu(p^i)=0μ(pi)=0，有：
F(pk)=∑d∣pkμ(d)=μ(1)+μ(p)+μ(p2)+⋯+μ(pk)=1+(−1)+0+⋯+0=0F(p^k)=\sum_{d|p^k}\mu(d)=\mu(1)+\mu(p)+\mu(p^2)+\cdots+\mu(p^k)=1+(-1)+0+\cdots+0=0F(pk)=d∣pk∑μ(d)=μ(1)+μ(p)+μ(p2)+⋯+μ(pk)=1+(−1)+0+⋯+0=0
下面是莫比乌斯函数线性筛代码，求1～n1～n1～n内的莫比乌斯函数，和欧拉函数线性筛的代码几乎一样，此处不做解释。
bool vis[maxn];
int prime[maxn];
int Mob[maxn];
void Mobius_sieve(){int cnt = 0;vis[1] = 1;Mob[1] = 1;for(int i = 2; i <= maxn; i++){if(!vis[i])prime[cnt++] = i, Mob[i] = - 1;for(int j = 0; j < cnt && 1LL * prime[j] * i <= maxn; j++){vis[prime[j] * i] = 1;Mob[i * prime[j]] = (i % prime[j] ? -Mob[i]: 0);if(i % prime[j] == 0)break;}}
}
5.2 莫比乌斯函数的由来
看到莫比乌斯函数的定义，读者是不是感到奇怪？它是怎么来的？有什么用？
本节内容，引用《初等数论及其应用》199页，它给出了莫比乌斯函数的来龙去脉。
设fff为算术函数，fff的和函数FFF的值为F(n)=∑d∣nf(d)F(n)=\sum_{d|n}f(d)F(n)=∑d∣nf(d)，显然，FFF由fff的值决定。这种关系可以反过来吗？也就是说，是否存在一种用FFF求出fff的值的简便方法？本节给出这样的公式，看什么样的公式可行。
若fff是算术函数，FFF是它的和函数F(n)=∑d∣nf(d)F(n)=\sum_{d|n}f(d)F(n)=∑d∣nf(d)，按定义展开F(n)，n=1,2,...,8F(n)，n=1,2,...,8F(n)，n=1,2,...,8，有：
F(1)=f(1)F(1) = f(1)F(1)=f(1)
F(2)=f(1)+f(2)F(2) = f(1) + f(2)F(2)=f(1)+f(2)
F(3)=f(1)+f(3)F(3) = f(1) + f(3)F(3)=f(1)+f(3)
F(4)=f(1)+f(2)+f(4)F(4) = f(1) + f(2) + f(4)F(4)=f(1)+f(2)+f(4)
F(5)=f(1)+f(5)F(5) = f(1) + f(5)F(5)=f(1)+f(5)
F(6)=f(1)+f(2)+f(3)+f(6)F(6) = f(1) + f(2) + f(3) + f(6)F(6)=f(1)+f(2)+f(3)+f(6)
F(7)=f(1)+f(7)F(7) = f(1) + f(7)F(7)=f(1)+f(7)
F(8)=f(1)+f(2)+f(4)+f(8)F(8) = f(1) + f(2) + f(4) + f(8)F(8)=f(1)+f(2)+f(4)+f(8)
根据上面方程，可以解得：
f(1)=F(1)f(1) = F(1)f(1)=F(1)
f(2)=F(2)−F(1)f(2) = F(2) - F(1)f(2)=F(2)−F(1)
f(3)=F(3)−F(1)f(3) = F(3) - F(1)f(3)=F(3)−F(1)
f(4)=F(4)−F(2)f(4) = F(4) - F(2)f(4)=F(4)−F(2)
f(5)=F(5)−F(1)f(5) = F(5) - F(1)f(5)=F(5)−F(1)
f(6)=F(6)−F(3)−F(2)+F(1)f(6) = F(6) - F(3) - F(2) + F(1)f(6)=F(6)−F(3)−F(2)+F(1)
f(7)=F(7)−F(1)f(7) = F(7) - F(1)f(7)=F(7)−F(1)
f(8)=F(8)−F(4)f(8) = F(8) - F(4)f(8)=F(8)−F(4)
注意到f(n)f(n)f(n)等于形式为±F(nd)\pm F(\frac nd)±F(dn)的一些项之和，其中d∣nd|nd∣n。从这一结果中，可能有这样一个等式，形式是：
f(n)=∑d∣nμ(d)F(nd)f(n)=\sum_{d|n}\mu(d)F(\frac nd)f(n)=d∣n∑μ(d)F(dn)
其中μ\muμ是算术函数。如果等式成立，计算得到：μ(1)=1,μ(2)=−1,μ(3)=−1,μ(4)=0,μ(5)=−1,μ(6)=1,μ(7)=−1,μ(8)=0\mu(1)=1,\mu(2)=-1,\mu(3)=-1,\mu(4)=0,\mu(5)=-1,\mu(6)=1,\mu(7)=-1,\mu(8)=0μ(1)=1,μ(2)=−1,μ(3)=−1,μ(4)=0,μ(5)=−1,μ(6)=1,μ(7)=−1,μ(8)=0。（读者已经发现和前一节莫比乌斯函数的值一样。）
又F(p)=f(1)+f(p)F(p)=f(1)+f(p)F(p)=f(1)+f(p)，得f(p)=F(p)−F(1)f(p)=F(p)-F(1)f(p)=F(p)−F(1)，其中ppp是质数。则μ(p)=−1\mu(p)=-1μ(p)=−1。
又因为F(p2)=f(1)+f(p)+f(p2)F(p^2)=f(1)+f(p)+f(p^2)F(p2)=f(1)+f(p)+f(p2)
有：f(p2)=F(p2)−(F(p)−F(1))−F(1)=F(p2)−F(p)f(p^2)=F(p^2)-(F(p)-F(1))-F(1)=F(p^2)-F(p)f(p2)=F(p2)−(F(p)−F(1))−F(1)=F(p2)−F(p)
这要求对任意质数ppp，有μ(p2)=0\mu(p^2)=0μ(p2)=0。类似的推理得出对任意质数ppp及整数k>1k>1k>1，有μ(pk)=0\mu(p^k)=0μ(pk)=0。如果猜想μ\muμ是积性函数，则μ\muμ的值就由质数幂处的值决定。
这就是莫比乌斯函数的由来。
5.3 莫比乌斯反演
【定理5.2】莫比乌斯反演(Möbius inversion)公式。若fff是算术函数，FFF为fff的和函数，对任意正整数nnn，满足F(n)=∑d∣nf(d)F(n)=\sum_{d|n}f(d)F(n)=∑d∣nf(d)，则对任意正整数nnn，有f(n)=∑d∣nμ(d)F(nd)f(n)=\sum_{d|n}\mu(d)F(\frac nd)f(n)=∑d∣nμ(d)F(dn)。
从定理5.2可以推出下面的定理5.3。
【定理5.3】 设fff是算术函数，它的和函数为F(n)=∑d∣nf(d)F(n)=\sum_{d|n}f(d)F(n)=∑d∣nf(d)，那么如果FFF是积性函数，则fff也是积性函数。
莫比乌斯反演，不要求fff是积性函数。
6 杜教筛
6.1 杜教筛公式
6.1.1 杜教筛公式的推导
杜教筛要解决的是这一类问题：设f(n)f(n)f(n)是一个数论函数，计算S(n)=∑i=1nf(i)S(n)=\sum_{i=1}^nf(i)S(n)=∑i=1nf(i)。
由于nnn很大，要求以低于O(n)O(n)O(n)的复杂度求解，所以用前面提到的线性筛是不够的。如果能用到整除分块，每一块的值相等一起算，就加快了速度，也就是构造形如∑i=1n⌊ni⌋\sum_{i=1}^n\lfloor\frac{n}{i}\rfloor∑i=1n⌊in⌋的整除分块。
杜教筛的思路，简单地说，是把f(n)构造成能利用整除分块的新的函数，这个构造用到了卷积。
记S(n)=∑i=1nf(i)S(n)=\sum_{i=1}^nf(i)S(n)=∑i=1nf(i)。根据f(n)f(n)f(n)的性质，构造一个S(n)S(n)S(n)关于S(⌊ni⌋)S(\lfloor\frac{n}{i}\rfloor)S(⌊in⌋)的递推式，下面是构造方法11。
构造2个积性函数hhh和ggg，满足卷积：h=g∗fh=g*fh=g∗f。
根据卷积的定义，有h(i)=∑d∣ig(d)f(id)h(i)= \sum_{d|i}g(d)f(\frac{i}{d})h(i)=∑d∣ig(d)f(di)，求和：
∑i=1nh(i)=∑i=1n∑d∣iig(d)f(id)\sum_{i=1}^nh(i)=\sum_{i=1}^n\sum_{d|i}^ig(d)f(\frac{i}{d})i=1∑nh(i)=i=1∑nd∣i∑ig(d)f(di)
=∑d=1ng(d)∑d∣inf(id)\qquad\qquad =\sum_{d=1}^ng(d)\sum_{d|i}^nf(\frac{i}{d})=d=1∑ng(d)d∣i∑nf(di)
=∑d=1ng(d)∑i=1⌊nd⌋f(i)\qquad\qquad =\sum_{d=1}^ng(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}f(i)=d=1∑ng(d)i=1∑⌊dn⌋f(i)
=∑d=1ng(d)s(⌊nd⌋)\qquad\qquad =\sum_{d=1}^ng(d)s({\lfloor\frac{n}{d}\rfloor})=d=1∑ng(d)s(⌊dn⌋)
注意其中第2行变换∑i=1n∑d∣iig(d)f(id)=∑d=1ng(d)∑d∣inf(id)\sum_{i=1}^n\sum_{d|i}^ig(d)f(\frac{i}{d})=\sum_{d=1}^ng(d)\sum_{d|i}^nf(\frac{i}{d})∑i=1n∑d∣iig(d)f(di)=∑d=1ng(d)∑d∣inf(di)，网上昵称“杜教筛变换” 11。
最后一步得到：∑i=1nh(i)=∑d=1ng(d)s(⌊nd⌋)\sum_{i=1}^nh(i)=\sum_{d=1}^ng(d)s({\lfloor\frac{n}{d}\rfloor})∑i=1nh(i)=∑d=1ng(d)s(⌊dn⌋)
为计算S(n)S(n)S(n)，把式子右边的第一项提出来：
∑i=1nh(i)=g(1)S(n)+∑d=2ng(d)s(⌊nd⌋)\sum_{i=1}^nh(i)=g(1)S(n)+\sum_{d=2}^ng(d)s({\lfloor\frac{n}{d}\rfloor})i=1∑nh(i)=g(1)S(n)+d=2∑ng(d)s(⌊dn⌋)
g(1)S(n)=∑i=1nh(i)−∑d=2ng(d)s(⌊nd⌋)g(1)S(n)=\sum_{i=1}^nh(i)-\sum_{d=2}^ng(d)s({\lfloor\frac{n}{d}\rfloor})g(1)S(n)=i=1∑nh(i)−d=2∑ng(d)s(⌊dn⌋)
式中的iii和ddd无关，为方便看，改写为：
g(1)S(n)=∑i=1nh(i)−∑i=2ng(i)s(⌊ni⌋)g(1)S(n)=\sum_{i=1}^nh(i)-\sum_{i=2}^ng(i)s({\lfloor\frac{n}{i}\rfloor})g(1)S(n)=i=1∑nh(i)−i=2∑ng(i)s(⌊in⌋)
这就是杜教筛公式。
6.1.2 杜教筛算法和复杂度
一个积性函数求和问题，如果分析得到了杜教筛公式，工作已经完成了一大半，剩下的是编程实现。
公式g(1)S(n)=∑i=1nh(i)−∑i=2ng(i)s(⌊ni⌋)g(1)S(n)=\sum_{i=1}^nh(i)-\sum_{i=2}^ng(i)s({\lfloor\frac{n}{i}\rfloor})g(1)S(n)=∑i=1nh(i)−∑i=2ng(i)s(⌊in⌋)是S(n)S(n)S(n)的递归形式，编程时可以递归实现；递归的时候加上记忆化，让每个S[i]S[i]S[i]只需要算一次。
其中的h(i)h(i)h(i)和g(i)g(i)g(i)，如果构造得好，就容易计算，见后面欧拉函数和莫比乌斯函数的例子。为方便分析计算复杂度，下面略去h(i)h(i)h(i)和g(i)g(i)g(i)，分析简化后的式子S(n)=∑i=2ns(⌊ni⌋)S(n)=\sum_{i=2}^ns({\lfloor\frac{n}{i}\rfloor})S(n)=∑i=2ns(⌊in⌋)的复杂度。
设S(n)S(n)S(n)的复杂度是T(n)T(n)T(n)。
递归展开第一层：S(n)=∑i=2ns(⌊ni⌋)=s(⌊n2⌋)+s(⌊n3⌋)+⋯+s(⌊nn⌋)S(n)=\sum_{i=2}^ns({\lfloor\frac{n}{i}\rfloor})=s({\lfloor\frac{n}{2}\rfloor})+s({\lfloor\frac{n}{3}\rfloor})+\cdots+s({\lfloor\frac{n}{n}\rfloor})S(n)=∑i=2ns(⌊in⌋)=s(⌊2n⌋)+s(⌊3n⌋)+⋯+s(⌊nn⌋)
根据整除分块的原理，右边共有O(n)O(\sqrt{n})O(n)种s(⌊ni⌋)s({\lfloor\frac{n}{i}\rfloor})s(⌊in⌋)。另外，把它们加起来的时候，还有O(n)O(\sqrt{n})O(n)次合并。总时间是T(n)=O(n)+O(∑i=2nT(ni))T(n)=O(\sqrt{n})+O(\sum_{i=2}^{\sqrt {n}}T(\frac ni))T(n)=O(n)+O(∑i=2nT(in))。
再展开一层：T(ni)=O(ni)+O(∑k=2niT(nik))=O(ni)T(\frac ni)=O(\sqrt{\frac ni})+O(\sum_{k=2}^{\sqrt {\frac ni}}T(\frac {\frac ni}k))=O(\sqrt{\frac ni})T(in)=O(in)+O(∑k=2inT(kin))=O(in)，中间第2项是高阶小量，把它略去。
合起来：
T(n)=O(n)+O(∑i=2nT(ni))T(n)=O(\sqrt{n})+O(\sum_{i=2}^{\sqrt {n}}T(\frac ni))T(n)=O(n)+O(i=2∑nT(in))
=O(n)+O(∑i=2nni)\qquad\qquad=O(\sqrt{n})+O(\sum_{i=2}^{\sqrt {n}}\sqrt{\frac ni})=O(n)+O(i=2∑nin)
后一项最大，只考虑它就够了：
T(n)=∑i=2nO(ni)≈O(∫0nnxdx)=O(n34)T(n)=\sum_{i=2}^{\sqrt {n}}O(\sqrt{\frac ni})\approx O(\int_0^{\sqrt{n}} {\sqrt{\frac nx}} \,{\rm d}x)=O(n^\frac 34)T(n)=i=2∑nO(in)≈O(∫0nxndx)=O(n43)
上述计算还能优化，即预处理一部分s(n)s(n)s(n)。s(n)s(n)s(n)是一个积性函数的前缀和，先用线性筛计算一部分，假设已经预处理了前kkk个正整数的s(n)s(n)s(n)，且k≥nk\geq \sqrt{n}k≥n。因为S(1)S(1)S(1)到S(k)S(k)S(k)已经求出，那么这个式子中还没有求出的是k<⌊ni⌋≤nk<{\lfloor\frac{n}{i}\rfloor}\leq{n}k<⌊in⌋≤n，也就是要算1<i≤nk1<i\leq{\frac nk}1<i≤kn的这一部分。那么复杂度变为：
T(n)=∑i=2nkO(ni)≈O(∫0nknxdx)=O(nk)T(n)=\sum_{i=2}^{\frac nk}O(\sqrt{\frac ni})\approx O(\int_0^{\frac nk} {\sqrt{\frac nx}} \,{\rm d}x)=O(\frac n{\sqrt k})T(n)=i=2∑knO(in)≈O(∫0knxndx)=O(kn)
kkk取多大合适呢？线性筛计算也是要花时间的，而且是线性的O(k)O(k)O(k)，那么线性筛计算加上杜教筛公式计算，总时间是：
T′(n)=O(k)+T(n)=O(k)+O(nk)T'(n)=O(k)+T(n)=O(k)+O(\frac n{\sqrt k})T′(n)=O(k)+T(n)=O(k)+O(kn)
差不多k=n23k=n^{\frac 23}k=n32时，它有最小值，即：O(n23)+O(nn23)=O(n23)O(n^{\frac 23})+O(\frac n{\sqrt{n^{\frac 23}}})=O(n^{\frac 23})O(n32)+O(n32n)=O(n32)。
它比O(n34)O(n^{\frac 34})O(n43)要好。
以上就是杜教筛算法的全部。
综合起来，杜教筛算法可以概况为：用狄利克雷卷积构造递推式，编程时用整除分块和线性筛优化，算法复杂度达到了O(n23)O(n^{\frac 23})O(n32)。
杜教筛算法 = 狄利克雷卷积 + 整除分块 + 线性筛
应用杜教筛时，狄利克雷卷积是基本的工具。观察卷积的定义(f∗g)(n)=∑d∣nf(d)g(nd)(f*g)(n)= \sum_{d|n}f(d)g(\frac{n}{d})(f∗g)(n)=∑d∣nf(d)g(dn)，那么，如果求解的函数有形如∑d∣nf(d))\sum_{d|n}f(d))∑d∣nf(d))这样的性质，把这个性质与卷积的定义对照，就容易找到ggg和hhh。下面用欧拉函数和莫比乌斯函数为例说明。算法的具体实现，见6.4节的模板代码。
6.2 欧拉函数前缀和
求欧拉函数前缀和：∑i=1nϕ(i)\sum_{i=1}^n\phi(i)∑i=1nϕ(i)。
求欧拉函数前缀和的杜教筛公式，需要用到欧拉函数性质：n=∑d∣nϕ(d)n=\sum_{d|n}\phi(d)n=∑d∣nϕ(d)。
（1）直接套用杜教筛公式
按卷积定义h=f∗g=∑d∣nf(d)g(nd)h=f*g= \sum_{d|n}f(d)g(\frac{n}{d})h=f∗g=∑d∣nf(d)g(dn)，把n=∑d∣nϕ(d)n=\sum_{d|n}\phi(d)n=∑d∣nϕ(d)改写为：id=ϕ∗Iid=\phi*Iid=ϕ∗I，那么有h=id，g=Ih=id，g=Ih=id，g=I。代入杜教筛公式g(1)S(n)=∑i=1nh(i)−∑i=2ng(i)s(⌊ni⌋)g(1)S(n)=\sum_{i=1}^nh(i)-\sum_{i=2}^ng(i)s({\lfloor\frac{n}{i}\rfloor})g(1)S(n)=∑i=1nh(i)−∑i=2ng(i)s(⌊in⌋)，得：
S(n)=∑i=1ni−∑i=2ns(⌊ni⌋)S(n)=\sum_{i=1}^ni-\sum_{i=2}^ns({\lfloor\frac{n}{i}\rfloor})S(n)=i=1∑ni−i=2∑ns(⌊in⌋)
=n(n+1)2−∑i=2ns(⌊ni⌋)\qquad\qquad=\frac{n(n+1)}{2}-\sum_{i=2}^ns({\lfloor\frac{n}{i}\rfloor})=2n(n+1)−i=2∑ns(⌊in⌋)
（2）自己推导
读者如果有兴趣，也可以自己按部就班地以欧拉函数为例，做一次杜教筛公式推导的练习，推导过程是12：
n=ϕ(n)+∑d∣n,d<nϕ(d)n=\phi(n)+\sum_{d|n,d<n}\phi(d)n=ϕ(n)+d∣n,d<n∑ϕ(d)
ϕ(n)=n−∑d∣n,d<nϕ(d)\phi(n)=n-\sum_{d|n,d<n}\phi(d)ϕ(n)=n−d∣n,d<n∑ϕ(d)
∑i=1nϕ(i)=∑i=1n(i−∑d∣i,d<iϕ(d))\sum_{i=1}^n\phi(i)=\sum_{i=1}^n(i-\sum_{d|i,d<i}\phi(d))i=1∑nϕ(i)=i=1∑n(i−d∣i,d<i∑ϕ(d))
∑i=1nϕ(i)=n(n+1)2−∑i=1n∑d∣i,d<iϕ(d)\sum_{i=1}^n\phi(i)=\frac{n(n+1)}{2}-\sum_{i=1}^n\sum_{d|i,d<i}\phi(d)i=1∑nϕ(i)=2n(n+1)−i=1∑nd∣i,d<i∑ϕ(d)
现在改变枚举变量：枚举id\dfrac iddi的值，即枚举iii对ddd的倍数，因为i≠di≠di=d，所以从2开始：
∑i=1nϕ(i)=n(n+1)2−∑id=2id≤n∑d=1d≤⌊nid⌋ϕ(d)\sum_{i=1}^n\phi(i)=\frac{n(n+1)}{2}-\sum_{\frac id=2}^{\frac id \leq n}\sum_{d=1}^{ d \leq {\lfloor\frac{n}{\frac id}\rfloor} }\phi(d)i=1∑nϕ(i)=2n(n+1)−di=2∑di≤nd=1∑d≤⌊din⌋ϕ(d)
设S(n)=∑i=1nϕ(i)S(n)=\sum_{i=1}^n\phi(i)S(n)=∑i=1nϕ(i)，并把id\dfrac iddi写成i′i'i′：
S(n)=n(n+1)2−∑i’=2ns(⌊ni‘⌋)S(n)=\frac{n(n+1)}{2}-\sum_{i’=2}^ns({\lfloor\frac{n}{i‘}\rfloor})S(n)=2n(n+1)−i’=2∑ns(⌊i‘n⌋)
这样就得到了欧拉函数的“杜教筛公式”，和（1）一样。
6.3 莫比乌斯函数前缀和
求莫比乌斯函数前缀和：∑i=1nμ(i)\sum_{i=1}^n\mu(i)∑i=1nμ(i)
求莫比乌斯函数前缀和的杜教筛公式，方法和上一节欧拉函数几乎一样。
需要用到莫比乌斯函数性质：∑d∣nμ(d)={1,n=10,n>1\sum_{d|n}\mu(d)= \begin{cases} 1, & \text {n=1} \\ 0, & \text{n>1} \end{cases}∑d∣nμ(d)={1,0,n=1n>1，为方便书写，改写成∑d∣nμ(d)=[n=1]\sum_{d|n}\mu(d)=[n=1]∑d∣nμ(d)=[n=1]。
（1）套用杜教筛公式
按卷积定义h=f∗g=∑d∣nf(d)g(nd)h=f*g= \sum_{d|n}f(d)g(\frac{n}{d})h=f∗g=∑d∣nf(d)g(dn)，把∑d∣nμ(d)=[n=1]\sum_{d|n}\mu(d)=[n=1]∑d∣nμ(d)=[n=1]改写为：ϵ=μ∗I\epsilon=\mu*Iϵ=μ∗I，那么有h=ϵ，g=Ih=\epsilon，g=Ih=ϵ，g=I。代入杜教筛公式g(1)S(n)=∑i=1nh(i)−∑i=2ng(i)s(⌊ni⌋)g(1)S(n)=\sum_{i=1}^nh(i)-\sum_{i=2}^ng(i)s({\lfloor\frac{n}{i}\rfloor})g(1)S(n)=∑i=1nh(i)−∑i=2ng(i)s(⌊in⌋)，得：
S(n)=∑i=1nϵ(i)−∑i=2nS(⌊ni⌋)S(n)=\sum_{i=1}^n\epsilon(i)-\sum_{i=2}^nS({\lfloor\frac{n}{i}\rfloor})S(n)=i=1∑nϵ(i)−i=2∑nS(⌊in⌋)
=1−∑i=2nS(⌊ni⌋)\qquad\qquad=1-\sum_{i=2}^nS({\lfloor\frac{n}{i}\rfloor})=1−i=2∑nS(⌊in⌋)
（2）自己推导13
[n=1]=μ(n)+∑d∣n,d<nμ(d)[n=1]=\mu(n)+\sum_{d|n,d<n}\mu(d)[n=1]=μ(n)+d∣n,d<n∑μ(d)
μ(n)=[n=1]−∑d∣n,d<nμ(d)\mu(n)=[n=1]-\sum_{d|n,d<n}\mu(d)μ(n)=[n=1]−d∣n,d<n∑μ(d)
∑i=1nμ(i)=1−∑i=1n∑d∣i,d<iμ(d)\sum_{i=1}^n\mu(i)=1-\sum_{i=1}^n\sum_{d|i,d<i}\mu(d)i=1∑nμ(i)=1−i=1∑nd∣i,d<i∑μ(d)
∑i=1nμ(i)=1−∑id=2id≤n∑d=1d≤⌊nid⌋∑i=1nμ(d)\sum_{i=1}^n\mu(i)=1- \sum_{\frac id =2}^ {\frac id \leq n} \sum_{d=1}^{d\leq {\lfloor\frac{n}{\frac id}\rfloor} } \sum_{i=1}^n\mu(d)i=1∑nμ(i)=1−di=2∑di≤nd=1∑d≤⌊din⌋i=1∑nμ(d)
令S(n)=∑i=1nμ(i)S(n)=\sum_{i=1}^n\mu(i)S(n)=∑i=1nμ(i)，并把id\dfrac iddi写成i′i'i′，得：
S(n)=1−∑i‘=2nS(⌊ni’⌋)S(n)=1-\sum_{i‘=2}^nS({\lfloor\frac{n}{i’}\rfloor})S(n)=1−i‘=2∑nS(⌊i’n⌋)
6.4 模板代码
下面给出洛谷P4213的代码14。
注释中说明了杜教筛的三个技术：整除分块、线性筛、狄利克雷卷积。其中整除分块、线性筛的代码和前面有关的代码一样。
//代码改写自：https://blog.csdn.net/KIKO_caoyue/article/details/100061406
#include <bits/stdc++.h>
using namespace std;typedef long long ll;
const int maxn = 5e6+7; //超过n^(2/3)，够大了int prime[maxn]; //记录质数
bool vis[maxn]; //记录是否被筛；
int mu[maxn]; //莫比乌斯函数值
ll phi[maxn]; //欧拉函数值
unordered_map<int,int> summu; //莫比乌斯函数前缀和
unordered_map<int,ll> sumphi; //欧拉函数前缀和
void init(){ //线性筛预计算一部分答案int cnt = 0;vis[0] = vis[1] = 1;mu[1] = phi[1] = 1;for(int i=2;i<maxn;i++){if(!vis[i]){prime[cnt++] = i;mu[i] = -1;phi[i] = i-1;}for(int j=0;j<cnt && i*prime[j]<maxn;j++){vis[i*prime[j]] = 1;if(i%prime[j]){mu[i*prime[j]] = -mu[i];phi[i*prime[j]] = phi[i]*phi[prime[j]];}else{mu[i*prime[j]] = 0;phi[i*prime[j]] = phi[i]*prime[j];break;}}}for(int i=1;i<maxn;i++){ //最后，mu[]和phi[]改为记录1~maxn的前缀和。mu[i] += mu[i-1];phi[i] += phi[i-1];}
}
int gsum(int x){ // g(i)的前缀和return x;
}
ll getsmu(int x){if(x < maxn) return mu[x]; //预计算if(summu[x]) return summu[x]; //记忆化ll ans = 1; //杜教筛公式中的 1for(ll l=2,r;l<=x;l=r+1){ //用整除分块计算杜教筛公式r = x/(x/l);ans -= (gsum(r)-gsum(l-1))*getsmu(x/l);}return summu[x] = ans/gsum(1);
}
ll getsphi(int x){if(x < maxn) return phi[x];if(sumphi[x]) return sumphi[x]; //记忆化，每个sumphi[x]只用算一次ll ans = x*((ll)x+1)/2; //杜教筛公式中的 n(n+1)/2for(ll l=2,r;l<=x;l=r+1){ //用整除分块计算杜教筛公式，这里算 sqrt(x)次r = x/(x/l);ans -= (gsum(r)-gsum(l-1))*getsphi(x/l);}return sumphi[x] = ans/gsum(1);
}
int main(){init(); //用线性筛预计算一部分int t;scanf("%d",&t);while(t--){int n;scanf("%d",&n);printf("%lld %lld\n",getsphi(n),getsmu(n));}return 0;
}
6.5 题目
这里列出网上的一些题目。
（1）求∑i=1niϕ(i)\sum_{i=1}^ni\phi(i)∑i=1niϕ(i)
题解：http://fzl123.top/archives/898
（2）求∑i=1n∑j=1nσ(ij)\sum_{i=1}^n\sum_{j=1}^n\sigma(ij)∑i=1n∑j=1nσ(ij)，其中n≤109n≤10^9n≤109 。
题解：https://www.cnblogs.com/zhugezy/p/11312301.html
（3）hdu 6706
http://acm.hdu.edu.cn/showproblem.php?pid=6706
（4）这里列出了一些可提交的题目：
https://www.cnblogs.com/TSHugh/p/8361040.html
（4）skywalkert（唐老师）的博客给出了很多习题：
https://blog.csdn.net/skywalkert/article/details/50500009
致谢
杜瑜皓：对本文草稿提出细致的修改意见。
华东理工大学队员：傅志凌（oj.ecustacm.cn管理员）、赵李洋、李震、彭博，讨论算法复杂度。
参考文献
[1] 积性函数求和的几种方法，任之洲，2016年信息学奥林匹克中国国家队候选队员论文集
[2] 唐靖哲：https://blog.csdn.net/skywalkert/article/details/50500009
[3] 吉如一：http://jiruyi910387714.is-programmer.com/posts/195270.html
[4] 傅志凌：http://fzl123.top/archives/898
本文的注脚：
《初等数论及其应用》第7章乘性函数 ↩︎
有乘性函数，也有加性函数，参考《初等数论及其应用》182页。 ↩︎
有关欧拉函数的所有证明，见《初等数论及其应用》176-180页。 ↩︎
[欧拉定理：设mmm是一个正整数，aaa是一个整数且(a,m)=1(a, m)=1(a,m)=1，那么aϕ(m)≡1(modm)a^{\phi(m)}\equiv1(mod\ m)aϕ(m)≡1(mod m)]。 ↩︎
试除法很低效，有很多快速分解质因子的方法，参考《初等数论及其应用》93页。 ↩︎
除法运算，向上取整，英文Ceiling，用符号 ⌈⌉\lceil \rceil⌈⌉表示；向下取整，英文Floor，用符号 ⌊⌋\lfloor \rfloor⌊⌋表示。 ↩︎
证明细节见：https://blog.csdn.net/qq_41021816/article/details/84842956。 ↩︎
参考：https://brilliant.org/wiki/dirichlet-convolution/。 ↩︎
https://brilliant.org/wiki/mobius-function/。 ↩︎
《初等数论及其应用》200页。 ↩︎
https://blog.csdn.net/myjs999/article/details/78906549 ↩︎ ↩︎
引用自：https://blog.csdn.net/qq_34454069/article/details/79492437 ↩︎
引用自：https://blog.csdn.net/qq_34454069/article/details/79492437 ↩︎
改写自：https://blog.csdn.net/KIKO_caoyue/article/details/100061406 ↩︎
杜教筛 以及积性函数的前世今生 --算法竞赛专题解析（4）相关推荐
- 杜教筛--51nod1239 欧拉函数之和
求$\sum_{i=1}^{n}\varphi (i)$,$n\leqslant 1e10$. 这里先把杜教筛的一般套路贴一下: 要求$S(n)=\sum_{i=1}^{n}f(i)$,而现在有一数论 ...
- ACM-ICPC 2018 南京赛区网络预赛Sum，线性筛处理积性函数
SUM 题意:f(n)是n可以拆成多少组n=a*b,a和b都是不包含平方因子的方案数目,对于a!=b,n=a*b和n=b*a算两种方案,求∑i=1nf(i) 首先我们可以知道,n=1时f(1)=1, ...
- [复习]莫比乌斯反演,杜教筛,min_25筛
[复习]莫比乌斯反演,杜教筛,min_25筛 莫比乌斯反演 做题的时候的常用形式: \[\begin{aligned}g(n)&=\sum_{n|d}f(d)\\f(n)&=\sum_ ...
- 杜教筛（上）：整除分块，积性函数，欧拉与莫比乌斯
整除分块: 当我们求∑i=1nf(⌊ni⌋)\sum_{i=1}^nf(\lfloor\frac{n}{i}\rfloor)∑i=1nf(⌊in⌋)的时候,如果1到n求一遍感觉太傻了,因为会有很多 ...
- 牛客 华华给月月出题 （积性函数+欧拉筛+快速幂）
题目描述 华华刚刚帮月月完成了作业.为了展示自己的学习水平之高超,华华还给月月出了一道类似的题: ⊕符号表示异或和,详见样例解释. 虽然月月写了个程序暴力的算出了答案,但是为了确保自己的答案没有错,希 ...
- 积性函数的性质及证明 + 线性筛
引言 在数论问题中,积性函数有着广泛的应用. 如在莫比乌斯反演问题中,函数变换之后如何快速维护前缀和往往是最重要也是最难的一步.如果维护的函数具有积性,那就可以尝试利用线性筛在O(n)O(n)O(n) ...
- 模板：杜教筛（莫比乌斯反演、数论）
所谓杜教筛,就是dms教给我们的筛 (逃) 前言 与其说算法,不如说是技巧. 可以在低于线性的时间复杂度(准确的说是 O(n23)O(n^{\frac{2}{3}})O(n32))内完成对积性函数的 ...
- 杜教筛及其时间复杂度分析
文章目录 杜教筛 方法 举例 莫比乌斯函数 欧拉函数 时间复杂度 杜教筛 杜教筛用于求一类积性函数的前缀和,时间复杂度可以做到 O(n23)O(n^{\frac{2}{3}})O(n32). 方法 ...
- [51NOD1847]奇怪的数学题（杜教筛+min_25筛+第二类斯特林数）
f(x)f(x)f(x)表示xxx的次大约数,有f(x)=xx的最小质因数f(x)=\frac{x}{x的最小质因数}f(x)=x的最小质因数x,那么 ∑i=1n∑j=1nsgcd(i,j)k=∑i ...
- matlab狄利克雷函数,数论入门1——莫比乌斯函数,欧拉函数,狄利克雷卷积,线性筛,莫比乌斯反演,杜教筛...
数论入门1 一个菜鸡对数论的一点点理解... 莫比乌斯函数 定义函数$\mu(n)$为: 当n有平方因子时,$\mu(n)=0$. 当n没有平方因子时,$\mu(n)=(-1)^{\omega(n)} ...
最新文章
- Composer使用
- hdu1007 最近点对
- Queue：poll、offer、element、peek
- 自己整理的openresty安装步骤
- ssm实训报告心得_Java开发学习心得（一）：SSM环境搭建
- 原生html冻结表头,CSS如何实现表头冻结效果
- jmu-Java-07多线程-互斥访问 (5分)
- 美国政府召开网络安全峰会，与私营行业巨头合力提振软件供应链和开源等安全...
- Leetcode 124.二叉树中的最大路径和
- 为什么计算机能读懂 1 和 0
- 大学生应该怎么学习Java？
- CDH集群更换ip,主机名
- 大数据工程师要学的编程_每个数据工程师都应了解的ml编程技巧，第2部分
- 无损音乐知识收集1（转）
- 怎样恢复内存卡的视频文件？（图文操作解析）
- 转载知乎大神设置普通路由器支持IPV6
- 下载最新的百度地图瓦片
- 建筑CAD基础设计【3】
- 北京大学百年讲堂内听果宁法师讲人生——提得起、放得下的深刻含义(摘抄)
- 每次压力大到爆，驾校教练总爱跑敬老院干这件事