杜教筛
杜教筛被用于处理一类数论函数的前缀和问题.对于数论函数 𝑓f,杜教筛可以在低于线性时间的复杂度内计算 𝑆(𝑛) =∑𝑛𝑖=1𝑓(𝑖)S(n)=∑i=1nf(i)
.
算法思想
我们想办法构造一个 𝑆(𝑛)S(n) 关于 𝑆(⌊𝑛𝑖⌋)S(⌊ni⌋)
的递推式.
对于任意一个数论函数 𝑔g,必满足:
𝑛∑𝑖=1(𝑓∗𝑔)(𝑖)=𝑛∑𝑖=1∑𝑑∣𝑖𝑔(𝑑)𝑓(𝑖𝑑)=𝑛∑𝑖=1𝑔(𝑖)𝑆(⌊𝑛𝑖⌋)∑i=1n(f∗g)(i)=∑i=1n∑d∣ig(d)f(id)=∑i=1ng(i)S(⌊ni⌋)
其中 𝑓 ∗𝑔f∗g 为数论函数 𝑓f
和 𝑔g
的 狄利克雷卷积.
略证
𝑔(𝑑)𝑓(𝑖𝑑)g(d)f(id) 就是对所有 𝑖 ≤𝑛i≤n
的做贡献,因此变换枚举顺序,枚举 𝑑d
,𝑖𝑑id
(分别对应新的 𝑖,𝑗i,j
)
𝑛∑𝑖=1∑𝑑∣𝑖𝑔(𝑑)𝑓(𝑖𝑑)=𝑛∑𝑖=1⌊𝑛/𝑖⌋∑𝑗=1𝑔(𝑖)𝑓(𝑗)=𝑛∑𝑖=1𝑔(𝑖)⌊𝑛/𝑖⌋∑𝑗=1𝑓(𝑗)=𝑛∑𝑖=1𝑔(𝑖)𝑆(⌊𝑛𝑖⌋)∑i=1n∑d∣ig(d)f(id)=∑i=1n∑j=1⌊n/i⌋g(i)f(j)=∑i=1ng(i)∑j=1⌊n/i⌋f(j)=∑i=1ng(i)S(⌊ni⌋)
那么可以得到递推式:
𝑔(1)𝑆(𝑛)=𝑛∑𝑖=1𝑔(𝑖)𝑆(⌊𝑛𝑖⌋)−𝑛∑𝑖=2𝑔(𝑖)𝑆(⌊𝑛𝑖⌋)=𝑛∑𝑖=1(𝑓∗𝑔)(𝑖)−𝑛∑𝑖=2𝑔(𝑖)𝑆(⌊𝑛𝑖⌋)g(1)S(n)=∑i=1ng(i)S(⌊ni⌋)−∑i=2ng(i)S(⌊ni⌋)=∑i=1n(f∗g)(i)−∑i=2ng(i)S(⌊ni⌋)
假如我们可以构造恰当的数论函数 𝑔g 使得:
- 可以快速计算 ∑𝑛𝑖=1(𝑓 ∗𝑔)(𝑖)∑i=1n(f∗g)(i)
;
- 可以快速计算 𝑔g
的前缀和,以用数论分块求解 ∑𝑛𝑖=2𝑔(𝑖)𝑆(⌊𝑛𝑖⌋)∑i=2ng(i)S(⌊ni⌋)
.
则我们可以在较短时间内求得 𝑔(1)𝑆(𝑛)g(1)S(n).
注意
无论数论函数 𝑓f 是否为积性函数,只要可以构造出恰当的数论函数 𝑔g
, 便都可以考虑用杜教筛求 𝑓f
的前缀和.
如考虑 𝑓(𝑛) =i𝜑(𝑛)f(n)=iφ(n), 显然 𝑓f
不是积性函数,但可取 𝑔(𝑛) =1g(n)=1
, 从而:
𝑛∑𝑘=1(𝑓∗𝑔)(𝑘)=i𝑛(𝑛+1)2∑k=1n(f∗g)(k)=in(n+1)2
计算 ∑𝑘≤𝑚(𝑓 ∗𝑔)(𝑘)∑k≤m(f∗g)(k) 和 ∑𝑘≤𝑚𝑔(𝑘)∑k≤mg(k)
的时间复杂度均为 𝑂(1)O(1)
, 故可以考虑使用杜教筛.
时间复杂度
令 𝑅(𝑛) ={⌊𝑛𝑘⌋:𝑘=2,3,…,𝑛}R(n)={⌊nk⌋:k=2,3,…,n}.利用数论分块的 性质 可知,对任意的 𝑚 ∈𝑅(𝑛)m∈R(n)
,都有 𝑅(𝑚) ⊆𝑅(𝑛)R(m)⊆R(n)
.也就是说,使用记忆化之后,只需要对所有 𝑘 ∈𝑅(𝑛)k∈R(n)
计算一次 𝑆(𝑘)S(k)
就可以得到 𝑅(𝑛)R(n)
的值.而这些点的数目 |𝑅(𝑛)| =𝑂(√𝑛)|R(n)|=O(n)
.
设计算 ∑𝑛𝑖=1(𝑓 ∗𝑔)(𝑖)∑i=1n(f∗g)(i) 和 ∑𝑛𝑖=1𝑔(𝑖)∑i=1ng(i)
的时间复杂度均为 𝑂(1)O(1)
. 设计算 𝑆(𝑛)S(n)
的时间复杂度为 𝑇(𝑛)T(n)
, 则:
𝑇(𝑛)=∑𝑘∈𝑅(𝑛)𝑇(𝑘)=Θ(√𝑛)+⌊√𝑛⌋∑𝑘=1𝑂(√𝑘)+⌊√𝑛⌋∑𝑘=2𝑂(√𝑛𝑘)=𝑂(∫√𝑛0(√𝑥+√𝑛𝑥)d𝑥)=𝑂(𝑛3/4).T(n)=∑k∈R(n)T(k)=Θ(n)+∑k=1⌊n⌋O(k)+∑k=2⌊n⌋O(nk)=O(∫0n(x+nx)dx)=O(n3/4).
若我们可以预处理出一部分 𝑆(𝑘)S(k), 其中 𝑘 =1,2,…,𝑚k=1,2,…,m
,𝑚 ≥⌊√𝑛⌋m≥⌊n⌋
.设预处理的时间复杂度为 𝑇0(𝑚)T0(m)
,则此时的 𝑇(𝑛)T(n)
为:
𝑇(𝑛)=𝑇0(𝑚)+∑𝑘∈𝑅(𝑛);𝑘>𝑚𝑇(𝑘)=𝑇0(𝑚)+⌊𝑛/𝑚⌋∑𝑘=1𝑂(√𝑛𝑘)=𝑂(𝑇0(𝑚)+∫𝑛/𝑚0√𝑛𝑥d𝑥)=𝑂(𝑇0(𝑚)+𝑛√𝑚).T(n)=T0(m)+∑k∈R(n);k>mT(k)=T0(m)+∑k=1⌊n/m⌋O(nk)=O(T0(m)+∫0n/mnxdx)=O(T0(m)+nm).
若 𝑇0(𝑚) =𝑂(𝑚)T0(m)=O(m)(如线性筛),由均值不等式可知:当 𝑚 =Θ(𝑛2/3)m=Θ(n2/3)
时,𝑇(𝑛)T(n)
取得最小值 𝑂(𝑛2/3)O(n2/3)
.
伪证一例
设计算 𝑆(𝑛)S(n) 的复杂度为 𝑇(𝑛)T(n)
, 则有:
𝑇(𝑛)=Θ(√𝑛)+𝑂⎛⎜ ⎜⎝⌊√𝑛⌋∑𝑖=2𝑇(⌊𝑛𝑖⌋)⎞⎟ ⎟⎠T(n)=Θ(n)+O(∑i=2⌊n⌋T(⌊ni⌋)) 𝑇(⌊𝑛𝑖⌋)=Θ(√𝑛𝑖)+𝑂⎛⎜ ⎜ ⎜⎝⌊√𝑛/𝑖⌋∑𝑗=2𝑇(⌊𝑛𝑖𝑗⌋)⎞⎟ ⎟ ⎟⎠=𝑂(√𝑛𝑖)T(⌊ni⌋)=Θ(ni)+O(∑j=2⌊n/i⌋T(⌊nij⌋))=O(ni)
其中,𝑂(∑⌊√𝑛/𝑖⌋𝑗=2𝑇(⌊𝑛𝑖𝑗⌋))O(∑j=2⌊n/i⌋T(⌊nij⌋)) 视作高阶无穷小,从而可以舍去.故:
𝑇(𝑛)=Θ(√𝑛)+𝑂⎛⎜ ⎜⎝⌊√𝑛⌋∑𝑖=2√𝑛𝑖⎞⎟ ⎟⎠=𝑂⎛⎜ ⎜⎝⌊√𝑛⌋∑𝑖=1√𝑛𝑖⎞⎟ ⎟⎠=𝑂(∫√𝑛0√𝑛𝑥d𝑥)=𝑂(𝑛3/4)T(n)=Θ(n)+O(∑i=2⌊n⌋ni)=O(∑i=1⌊n⌋ni)=O(∫0nnxdx)=O(n3/4) Bug
问题在于「视作高阶无穷小,从而可以舍去」这一处.我们将 𝑇(⌊𝑛𝑖⌋)T(⌊ni⌋) 代入 𝑇(𝑛)T(n)
的式子里,有:
𝑇(𝑛)=Θ(√𝑛)+𝑂⎛⎜ ⎜⎝⌊√𝑛⌋∑𝑖=2√𝑛𝑖⎞⎟ ⎟⎠+𝑂⎛⎜ ⎜ ⎜⎝⌊√𝑛⌋∑𝑖=2⌊√𝑛/𝑖⌋∑𝑗=2𝑇(⌊𝑛𝑖𝑗⌋)⎞⎟ ⎟ ⎟⎠=𝑂(√𝑛+∫√𝑛0√𝑛𝑥d𝑥)+𝑂⎛⎜ ⎜ ⎜⎝⌊√𝑛⌋∑𝑖=2⌊√𝑛/𝑖⌋∑𝑗=2𝑇(⌊𝑛𝑖𝑗⌋)⎞⎟ ⎟ ⎟⎠=𝑂(𝑛3/4)+𝑂⎛⎜ ⎜ ⎜⎝⌊√𝑛⌋∑𝑖=2⌊√𝑛/𝑖⌋∑𝑗=2𝑇(⌊𝑛𝑖𝑗⌋)⎞⎟ ⎟ ⎟⎠T(n)=Θ(n)+O(∑i=2⌊n⌋ni)+O(∑i=2⌊n⌋∑j=2⌊n/i⌋T(⌊nij⌋))=O(n+∫0nnxdx)+O(∑i=2⌊n⌋∑j=2⌊n/i⌋T(⌊nij⌋))=O(n3/4)+O(∑i=2⌊n⌋∑j=2⌊n/i⌋T(⌊nij⌋))
我们考虑 ⌊√𝑛⌋∑𝑖=2⌊√𝑛/𝑖⌋∑𝑗=2𝑇(⌊𝑛𝑖𝑗⌋)∑i=2⌊n⌋∑j=2⌊n/i⌋T(⌊nij⌋) 这部分,不难发现:
⌊√𝑛⌋∑𝑖=2⌊√𝑛/𝑖⌋∑𝑗=2𝑇(⌊𝑛𝑖𝑗⌋)=Ω⎛⎜ ⎜⎝⌊√𝑛⌋∑𝑖=2𝑇⎛⎜ ⎜⎝⎢ ⎢ ⎢⎣𝑛𝑖⋅⌊√𝑛𝑖⌋−1⎥ ⎥ ⎥⎦⎞⎟ ⎟⎠⎞⎟ ⎟⎠=Ω⎛⎜ ⎜⎝⌊√𝑛⌋∑𝑖=2𝑇(⌊√𝑛𝑖⌋)⎞⎟ ⎟⎠∑i=2⌊n⌋∑j=2⌊n/i⌋T(⌊nij⌋)=Ω(∑i=2⌊n⌋T(⌊ni⋅⌊ni⌋−1⌋))=Ω(∑i=2⌊n⌋T(⌊ni⌋))
由于没有引入记忆化,因此上式中的 𝑇(⌊√𝑛𝑖⌋)T(⌊ni⌋) 仍然是 Ω((𝑛𝑖)1/4)Ω((ni)1/4)
的,进而所谓的「高阶无穷小」部分是不可以舍去的.
实际上杜教筛的亚线性时间复杂度是由记忆化保证的.只有使用了记忆化之后才能保证不会出现那个多重求和的项.
例题
问题一
求 𝑆1(𝑛) =∑𝑛𝑖=1𝜇(𝑖)S1(n)=∑i=1nμ(i) 和 𝑆2(𝑛) =∑𝑛𝑖=1𝜑(𝑖)S2(n)=∑i=1nφ(i)
的值,1 ≤𝑛 <2311≤n<231
.
莫比乌斯函数前缀和欧拉函数前缀和
我们知道:
𝜖=[𝑛=1]=𝜇∗1=∑𝑑∣𝑛𝜇(𝑑)ϵ=[n=1]=μ∗1=∑d∣nμ(d)𝑆1(𝑛)=𝑛∑𝑖=1𝜖(𝑖)−𝑛∑𝑖=2𝑆1(⌊𝑛𝑖⌋)=1−𝑛∑𝑖=2𝑆1(⌊𝑛𝑖⌋)S1(n)=∑i=1nϵ(i)−∑i=2nS1(⌊ni⌋)=1−∑i=2nS1(⌊ni⌋)
时间复杂度的推导见 时间复杂度 一节.
对于较大的值,需要用 map/unordered_map 存下其对应的值,方便以后使用时直接使用之前计算的结果.
当然也可以用杜教筛求出 𝜑(𝑥)φ(x) 的前缀和,但是更好的方法是应用莫比乌斯反演.
莫比乌斯反演杜教筛
𝑛∑𝑖=1𝑛∑𝑗=1[gcd(𝑖,𝑗)=1]=𝑛∑𝑖=1𝑛∑𝑗=1∑𝑑∣𝑖,𝑑∣𝑗𝜇(𝑑)=𝑛∑𝑑=1𝜇(𝑑)⌊𝑛𝑑⌋2∑i=1n∑j=1n[gcd(i,j)=1]=∑i=1n∑j=1n∑d∣i,d∣jμ(d)=∑d=1nμ(d)⌊nd⌋2
由于题目所求的是 ∑𝑛𝑖=1∑𝑖𝑗=1[gcd(𝑖,𝑗) =1]∑i=1n∑j=1i[gcd(i,j)=1], 所以我们排除掉 𝑖 =1,𝑗 =1i=1,j=1
的情况,并将结果除以 22
即可.
观察到,只需求出莫比乌斯函数的前缀和,就可以快速计算出欧拉函数的前缀和了.时间复杂度 𝑂(𝑛23)O(n23).
求 𝑆(𝑛) =∑𝑛𝑖=1𝜑(𝑖)S(n)=∑i=1nφ(i).
同样的,𝜑 ∗1 =idφ∗1=id, 从而:
𝑆(𝑛)=𝑛∑𝑖=1𝑖−𝑛∑𝑖=2𝑆(⌊𝑛𝑖⌋)=12𝑛(𝑛+1)−𝑛∑𝑖=2𝑆(⌊𝑛𝑖⌋)S(n)=∑i=1ni−∑i=2nS(⌊ni⌋)=12n(n+1)−∑i=2nS(⌊ni⌋)
代码实现
---|---
### 问题二
[「LuoguP3768」简单的数学题](https://www.luogu.com.cn/problem/P3768)
大意:求
𝑛∑𝑖=1𝑛∑𝑗=1𝑖⋅𝑗⋅gcd(𝑖,𝑗)(mod𝑝)∑i=1n∑j=1ni⋅j⋅gcd(i,j)(modp)
其中 𝑛 ≤1010,5 ×108 ≤𝑝 ≤1.1 ×109n≤1010,5×108≤p≤1.1×109,𝑝p 是质数.
利用 𝜑 ∗1 =idφ∗1=id 做莫比乌斯反演化为:
𝑛∑𝑑=1𝐹2(⌊𝑛𝑑⌋)⋅𝑑2𝜑(𝑑)∑d=1nF2(⌊nd⌋)⋅d2φ(d)
其中 𝐹(𝑛) =12𝑛(𝑛 +1)F(n)=12n(n+1)
对 ∑𝑛𝑑=1𝐹(⌊𝑛𝑑⌋)2∑d=1nF(⌊nd⌋)2 做数论分块,𝑑2𝜑(𝑑)d2φ(d) 的前缀和用杜教筛处理:
𝑓(𝑛)=𝑛2𝜑(𝑛)=(id2𝜑)(𝑛)f(n)=n2φ(n)=(id2φ)(n)𝑆(𝑛)=𝑛∑𝑖=1𝑓(𝑖)=𝑛∑𝑖=1(id2𝜑)(𝑖)S(n)=∑i=1nf(i)=∑i=1n(id2φ)(i)
需要构造积性函数 𝑔g,使得 𝑓 ×𝑔f×g 和 𝑔g 能快速求和.
单纯的 𝜑φ 的前缀和可以用 𝜑 ∗1φ∗1 的杜教筛处理,但是这里的 𝑓f 多了一个 id2id2,那么我们就卷一个 id2id2 上去,让它变成常数:
𝑆(𝑛)=𝑛∑𝑖=1((id2𝜑)∗id2)(𝑖)−𝑛∑𝑖=2id2(𝑖)𝑆(⌊𝑛𝑖⌋)S(n)=∑i=1n((id2φ)∗id2)(i)−∑i=2nid2(i)S(⌊ni⌋)
化一下卷积:
((id2𝜑)∗id2)(𝑖)=∑𝑑∣𝑖(id2𝜑)(𝑑)id2(𝑖𝑑)=∑𝑑∣𝑖𝑑2𝜑(𝑑)(𝑖𝑑)2=∑𝑑∣𝑖𝑖2𝜑(𝑑)=𝑖2∑𝑑∣𝑖𝜑(𝑑)=𝑖2(𝜑∗1)(𝑖)=𝑖3((id2φ)∗id2)(i)=∑d∣i(id2φ)(d)id2(id)=∑d∣id2φ(d)(id)2=∑d∣ii2φ(d)=i2∑d∣iφ(d)=i2(φ∗1)(i)=i3
再化一下 𝑆(𝑛)S(n):
𝑆(𝑛)=𝑛∑𝑖=1((id2𝜑)∗id2)(𝑖)−𝑛∑𝑖=2id2(𝑖)𝑆(⌊𝑛𝑖⌋)=𝑛∑𝑖=1𝑖3−𝑛∑𝑖=2𝑖2𝑆(⌊𝑛𝑖⌋)=(12𝑛(𝑛+1))2−𝑛∑𝑖=2𝑖2𝑆(⌊𝑛𝑖⌋)S(n)=∑i=1n((id2φ)∗id2)(i)−∑i=2nid2(i)S(⌊ni⌋)=∑i=1ni3−∑i=2ni2S(⌊ni⌋)=(12n(n+1))2−∑i=2ni2S(⌊ni⌋)
分块求解即可.
代码实现
---|---
参考资料
- 任之洲,2016,《积性函数求和的几种方法》,2016 年信息学奥林匹克中国国家队候选队员论文
- 杜教筛的时空复杂度分析 - riteme.site
本页面最近更新: 2026/1/7 08:56:54,更新历史 发现错误?想一起完善?在 GitHub 上编辑此页! 本页面贡献者:StudyingFather, Tiphereth-A, hsfzLZH1, Ir1d, Enter-tainer, sshwy, c-forrest, Marcythm, MegaOwIer, Henry-ZHR, Xeonacid, Backl1ght, Great-designer, huayucaiji, kenlig, ksyx, Menci, Nanarikom, nanmenyangde, ouuan, purple-vine, shawlleyw, Sshwy 本页面的全部内容在CC BY-SA 4.0 和 SATA 协议之条款下提供,附加条款亦可能应用