P4213 【模板】杜教筛
Rigel
·
·
题解
杜教筛是一种能以 \mathcal{O}(n^{\frac{2}{3}}) 的时间复杂度快速计算数论函数前缀和的算法。
0x00 前置知识
0x01 积性函数
首先定义定义域为 \mathbb{N}_{+} 的函数为数论函数。
-
若数论函数 f 对于任意两个互质的正整数 p,q 均有 f(pq)=f(p)f(q)，则称 f 为积性函数。
-
若数论函数 f 对于任意两个正整数 p,q 均有 f(pq)=f(p)f(q)，则称 f 为完全积性函数。
-
若 f 为积性函数，则 f 的和函数 F(n)=\displaystyle\sum\limits_{d \mid n}f(d) 也是积性函数。
以下是一些常见的积性函数。
-
-
-
-
- $\sigma_0(n)$ 亦记作 $d(n)$ 或 $\tau(n)$。
- $\sigma_1(n)$ 亦记作 $\sigma(n)$。
-
-
0x02 欧拉函数
定义欧拉函数 \varphi(n) 为小于等于 n 的正整数中与 n 互质的数的个数，即：
\varphi(n)=\sum\limits_{i=1}^{n}[\gcd(i,n)=1].
欧拉函数是积性函数。
对于质数 p，有：
以下是欧拉函数的一些性质。
-
对于正整数 n，有 \displaystyle\sum\limits_{d \mid n}\varphi(d)=n。
-
欧拉定理：若 m \in \mathbb{N}_{+}，a \in \mathbb{Z}，\gcd(a,m)=1，则 a^{\varphi(m)}\equiv 1\pmod m。
-
设 n=\displaystyle\prod\limits_{i=1}^{k}p_i^{c_i} 为正整数 n 的标准分解式，则 \displaystyle\varphi(n)=n\prod\limits_{i=1}^{k}\left( 1-\frac{1}{p_i} \right)。
0x03 莫比乌斯函数
定义莫比乌斯函数 \mu(n) 为：
1, & n=1,\\
(-1)^k, & n=\displaystyle\prod\limits_{i=1}^{k}p_i,\\
0, & \text{otherwise}.
\end{cases}
其中 p_i 为不同的质数。
莫比乌斯函数是积性函数。
对于质数 p，有：
莫比乌斯函数的一个重要性质：\displaystyle\sum\limits_{d \mid n}\mu(d)=\varepsilon(n)。
0x04 狄利克雷卷积
设 f,g 为数论函数，记 f,g 的狄利克雷卷积为 f*g，定义为：
(f*g)(n)=\sum\limits_{d \mid n}f(d)g\left( \frac{n}{d} \right).
以下是狄利克雷卷积的一些性质。
- 交换律：f*g=g*f。
- 结合律：(f*g)*h=f*(g*h)。
- 分配律：(f+g)*h=f*h+g*h。
- 对于任意数论函数 f，有 f*\varepsilon=f。
- 两个积性函数的狄利克雷卷积仍然是积性函数。
0x05 数论分块
数论分块可以快速计算形如 \displaystyle\sum\limits_{i=1}^{n}f(i)g\left( \left\lfloor\frac{n}{i}\right\rfloor \right) 的式子，时间复杂度为 \mathcal{O}(\sqrt{n})。
注意前提条件是能 \mathcal{O}(1) 求出 \displaystyle\sum\limits_{i=l}^{r}f(i) 或已经预处理出 f 的前缀和。
将 \left\lfloor\dfrac{n}{i}\right\rfloor 按相同值分块，可以证明块的数量 \le 2\sqrt{n}。
若一个块的第一个数为 l，可以证明这个块的最后一个数 r=\left\lfloor\dfrac{n}{\lfloor\frac{n}{l}\rfloor}\right\rfloor。
下面给出计算 \displaystyle\sum\limits_{i=1}^{n}\left\lfloor\frac{n}{i}\right\rfloor 的代码。
int ans=0;
for(int l=1,r;l<=n;l=r+1){
r=n/(n/l);
ans+=(n/l)*(r-l+1);
}
0x10 杜教筛原理与应用
0x11 杜教筛公式
我们要计算的是 S(n)=\displaystyle\sum\limits_{i=1}^{n}f(i)，其中 f 为数论函数。
我们需要构造数论函数 g,h，满足 h=g*f。
则 h(i)=\displaystyle\sum\limits_{d \mid i}g(d)f\left( \frac{i}{d} \right)。
对 h(i) 求和，得：
\sum\limits_{i=1}^{n}h(i)
&=\sum\limits_{i=1}^{n}\sum\limits_{d \mid i}g(d)f\left( \frac{i}{d} \right)\\
&=\sum\limits_{d=1}^{n}g(d)\sum\limits_{\substack{i=1\\d \mid i}}^{n}f\left( \frac{i}{d} \right)\\
&=\sum\limits_{d=1}^{n}g(d)\sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}f(i)\\
&=\sum\limits_{d=1}^{n}g(d)S\left( \left\lfloor\frac{n}{d}\right\rfloor \right)\\
&=g(1)S(n)+\sum\limits_{d=2}^{n}g(d)S\left( \left\lfloor\frac{n}{d}\right\rfloor \right).
\end{aligned}
于是：
g(1)S(n)=\sum\limits_{i=1}^{n}h(i)-\sum\limits_{d=2}^{n}g(d)S\left( \left\lfloor\frac{n}{d}\right\rfloor \right).
即：
\boxed{g(1)S(n)=\sum\limits_{i=1}^{n}h(i)-\sum\limits_{i=2}^{n}g(i)S\left( \left\lfloor\frac{n}{i}\right\rfloor \right).}
0x12 求欧拉函数前缀和
我们知道 \displaystyle\sum\limits_{d \mid n}\varphi(d)=n，即 \mathrm{id}=I*\varphi。
将 f=\varphi，g=I，h=\mathrm{id} 代入杜教筛公式，得：
I(1)S_{\varphi}(n) &= \sum\limits_{i=1}^{n}\mathrm{id}(i)-\sum\limits_{i=2}^{n}I(i)S_{\varphi}\left( \left\lfloor\frac{n}{i}\right\rfloor \right)\\
S_{\varphi}(n) &= \frac{n(n+1)}{2}-\sum\limits_{i=2}^{n}S_{\varphi}\left( \left\lfloor\frac{n}{i}\right\rfloor \right).
\end{aligned}
0x13 求莫比乌斯函数前缀和
我们知道 \displaystyle\sum\limits_{d \mid n}\mu(d)=\varepsilon(n)，即 \varepsilon=I*\mu。
将 f=\mu，g=I，h=\varepsilon 代入杜教筛公式，得：
I(1)S_{\mu}(n) &= \sum\limits_{i=1}^{n}\varepsilon(i)-\sum\limits_{i=2}^{n}I(i)S_{\mu}\left( \left\lfloor\frac{n}{i}\right\rfloor \right)\\
S_{\mu}(n)&=1-\sum\limits_{i=2}^{n}S_{\mu}\left( \left\lfloor\frac{n}{i}\right\rfloor \right).
\end{aligned}
0x20 时间复杂度
考虑单次计算 S(n)=\displaystyle\sum\limits_{i=1}^{n}f(i) 的时间复杂度。
::::error[错误证明]
设计算 \displaystyle\sum\limits_{i=1}^{n}h(i) 与 \displaystyle\sum\limits_{i=1}^{n}g(i) 的时间复杂度均为 \mathcal{O}(1)。
设计算 S(n) 的复杂度为 T(n)，则：
T(n)=\mathcal{O}(\sqrt{n})+\mathcal{O}\left( \sum\limits_{i=2}^{\lfloor\sqrt{n}\rfloor}T\left( \left\lfloor\frac{n}{i}\right\rfloor \right) \right).
其中：
T\left( \left\lfloor\frac{n}{i}\right\rfloor \right)=\mathcal{O}\left( \sqrt{\frac{n}{i}} \right)+\mathcal{O}\left( \sum\limits_{j=2}^{\left\lfloor\sqrt{\frac{n}{i}}\right\rfloor}T\left( \left\lfloor\frac{n}{ij}\right\rfloor \right) \right).
将 \mathcal{O}\left( \displaystyle\sum\limits_{j=2}^{\left\lfloor\sqrt{\frac{n}{i}}\right\rfloor}T\left( \left\lfloor\frac{n}{ij}\right\rfloor \right) \right) 视作高阶无穷小，于是：
T\left( \left\lfloor\frac{n}{i}\right\rfloor \right)=\mathcal{O}\left( \sqrt{\frac{n}{i}} \right).
将其代入，得：
T(n)&=\mathcal{O}(\sqrt{n})+\mathcal{O}\left( \sum\limits_{i=2}^{\lfloor\sqrt{n}\rfloor}\sqrt{\frac{n}{i}} \right)\\
&=\mathcal{O}\left( \sum\limits_{i=1}^{\lfloor\sqrt{n}\rfloor}\sqrt{\frac{n}{i}} \right)\\
&=\mathcal{O}\left( \int_{1}^{\sqrt{n}}\sqrt{\frac{n}{x}}\,\mathrm{d}x \right)\\
&=\mathcal{O}(n^{\frac{3}{4}}).
\end{aligned}
:::warning[错误点]{open}
问题在于“将 \mathcal{O}\left( \displaystyle\sum\limits_{j=2}^{\left\lfloor\sqrt{\frac{n}{i}}\right\rfloor}T\left( \left\lfloor\frac{n}{ij}\right\rfloor \right) \right) 视作高阶无穷小”。
将 T\left( \left\lfloor\dfrac{n}{i}\right\rfloor \right) 代入 T(n) 的式子中，有：
T(n)&=\mathcal{O}(\sqrt{n})+\mathcal{O}\left( \sum\limits_{i=2}^{\lfloor\sqrt{n}\rfloor}\left( \sqrt{\frac{n}{i}}+\sum\limits_{j=2}^{\left\lfloor\sqrt{\frac{n}{i}}\right\rfloor}T\left( \left\lfloor\frac{n}{ij}\right\rfloor \right) \right) \right)\\
&=\mathcal{O}(\sqrt{n})+\mathcal{O}\left( \sum\limits_{i=2}^{\lfloor\sqrt{n}\rfloor}\sqrt{\frac{n}{i}} \right)+\mathcal{O}\left( \sum\limits_{i=2}^{\lfloor\sqrt{n}\rfloor}\sum\limits_{j=2}^{\left\lfloor\sqrt{\frac{n}{i}}\right\rfloor}T\left( \left\lfloor\frac{n}{ij}\right\rfloor \right) \right)\\
&=\mathcal{O}(n^{\frac{3}{4}})+\mathcal{O}\left( \sum\limits_{i=2}^{\lfloor\sqrt{n}\rfloor}\sum\limits_{j=2}^{\left\lfloor\sqrt{\frac{n}{i}}\right\rfloor}T\left( \left\lfloor\frac{n}{ij}\right\rfloor \right) \right).
\end{aligned}
考虑 \displaystyle\sum\limits_{i=2}^{\lfloor\sqrt{n}\rfloor}\sum\limits_{j=2}^{\left\lfloor\sqrt{\frac{n}{i}}\right\rfloor}T\left( \left\lfloor\frac{n}{ij}\right\rfloor \right)：
\sum\limits_{i=2}^{\lfloor\sqrt{n}\rfloor}\sum\limits_{j=2}^{\left\lfloor\sqrt{\frac{n}{i}}\right\rfloor}T\left( \left\lfloor\frac{n}{ij}\right\rfloor \right)&=\Omega\left( \sum\limits_{i=2}^{\lfloor\sqrt{n}\rfloor}T\left( \left\lfloor\frac{\frac{n}{i}}{\left\lfloor\sqrt{\frac{n}{i}}\right\rfloor}\right\rfloor \right) \right)\\
&=\Omega\left( \sum\limits_{i=2}^{\lfloor\sqrt{n}\rfloor}T\left( \left\lfloor\sqrt{\frac{n}{i}}\right\rfloor \right) \right).
\end{aligned}
注意到 T\left( \left\lfloor\sqrt{\dfrac{n}{i}}\right\rfloor \right) 在没有使用记忆化时是 \Omega\left( \left( \dfrac{n}{i} \right)^{\frac{1}{4}} \right) 的，因此不能将 \mathcal{O}\left( \displaystyle\sum\limits_{j=2}^{\left\lfloor\sqrt{\frac{n}{i}}\right\rfloor}T\left( \left\lfloor\frac{n}{ij}\right\rfloor \right) \right) 视作高阶无穷小。
:::
::::
::::success[正确证明]{open}
令 R(n)=\left\{\left\lfloor\dfrac{n}{i}\right\rfloor \mid i=2,3,\cdots,n\right\}。
由数论分块的性质，对任意 k \in R(n)，都有 R(k) \subseteq R(n)。
也就是说，使用记忆化之后，仅需对所有 k \in R(n) 计算一次 S(k)，之后就可以反复使用。而 \lvert R(n)\rvert = \mathcal{O}(\sqrt{n})。
设计算 \displaystyle\sum\limits_{i=1}^{n}h(i) 与 \displaystyle\sum\limits_{i=1}^{n}g(i) 的时间复杂度均为 \mathcal{O}(1)。
设计算 S(n) 的时间复杂度为 T(n)，则：
T(n)&=\sum\limits_{k \in R(n)}T(k)\\
&=\mathcal{O}(\sqrt{n})+\sum\limits_{k=1}^{\lfloor\sqrt{n}\rfloor}\mathcal{O}(\sqrt{k})+\sum\limits_{k=2}^{\lfloor\sqrt{n}\rfloor}\mathcal{O}\left( \sqrt{\frac{n}{k}} \right)\\
&=\mathcal{O}\left( \int_{1}^{\sqrt{n}}\left( \sqrt{x}+\sqrt{\frac{n}{x}} \right)\,\mathrm{d}x \right)\\
&=\mathcal{O}(n^{\frac{3}{4}}).
\end{aligned}
考虑预处理一部分 S(k) 以优化复杂度，其中 k \in [1,m]，m \ge \lfloor\sqrt{n}\rfloor。此时：
T(n)&=\mathcal{O}(m)+\sum\limits_{\substack{k \in R(n)\\k>m}}T(k)\\
&=\mathcal{O}(m)+\sum\limits_{k=1}^{\left\lfloor\frac{n}{m}\right\rfloor}\mathcal{O}\left( \sqrt{\frac{n}{k}} \right)\\
&=\mathcal{O}\left( m+\int_{1}^{\frac{n}{m}}\sqrt{\frac{n}{x}}\,\mathrm{d}x \right)\\
&=\mathcal{O}\left( m+\frac{n}{\sqrt{m}} \right).
\end{aligned}
当 m=\Theta(n^{\frac{2}{3}}) 时，T(n) 取到最小值：
\boxed{T(n)=\mathcal{O}\left( n^{\frac{2}{3}}+\frac{n}{\sqrt{n^{\frac{2}{3}}}} \right)=\mathcal{O}(n^{\frac{2}{3}}).}
::::
由于本题为多测，故整体时间复杂度为 \mathcal{O}\left( m+T \cdot\dfrac{n}{\sqrt{m}} \right)。（此处 T 为数据组数）
当 m =\Theta\left( (Tn)^{\frac{2}{3}} \right) 时，理论最优时间复杂度为 \mathcal{O}\left( (Tn)^{\frac{2}{3}}+T \cdot\dfrac{n}{(Tn)^{\frac{1}{3}}} \right)=\mathcal{O}\left( (Tn)^{\frac{2}{3}} \right)。
0x30 Code
首先预处理出 i\in[1,m] 的 \varphi(i),\mu(i) 及其前缀和 S_{\varphi}(n),S_{\mu}(n)。
利用数论分块递归求 S_{\varphi}(n) 以及 S_{\mu}(n)。具体实现见代码。
#include<bits/stdc++.h>
#define int long long
const int _N=15;
const int M=8e6+10;
using namespace std;
int T,m,ans1,ans2;
int _in[_N];
int phi[M],mu[M],s_phi[M],s_mu[M];
int cnt,p[M];
bool vis[M];
unordered_map<int,int>t_phi,t_mu;
inline int read(){
int ret=0,f=1;char ch=getchar();
while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
while(ch>='0'&&ch<='9')ret=ret*10+(ch&15),ch=getchar();
return ret*f;
}
void _init(){
phi[1]=mu[1]=1;
for(int i=2;i<=m;i++){
if(!vis[i]){
p[++cnt]=i;
phi[i]=i-1; // 质数 i 的 φ(i)=i-1
mu[i]=-1; // 质数 i 的 μ(i)=-1
}
for(int j=1;j<=cnt&&i*p[j]<=m;j++){
vis[i*p[j]]=1;
if(i%p[j]){ // p_j 与 i 互质
phi[i*p[j]]=phi[i]*(p[j]-1);
// φ(i*p_j)=φ(i)*φ(p_j)=φ(i)*(p_j-1)
mu[i*p[j]]=-mu[i];
// μ(i*p_j)=μ(i)*μ(p_j)=μ(i)*(-1)=-μ(i)
}else{ // p_j 是 i 的因数，设 i=k*p_j
phi[i*p[j]]=phi[i]*p[j];
// φ(i*p_j)=φ(k*((p_j)^2))=φ(k*p_j)*p_j=φ(i)*p_j
mu[i*p[j]]=0;
// μ(i*p_j)=μ(k*((p_j)^2))=0
break;
}
}
}
for(int i=1;i<=m;i++)s_phi[i]=s_phi[i-1]+phi[i];
for(int i=1;i<=m;i++)s_mu[i]=s_mu[i-1]+mu[i];
}
int calc_phi(int x){
if(x<=m)return s_phi[x];
if(t_phi.count(x))return t_phi[x];
int ret=(x*(x+1))>>1;
for(int l=2,r;l<=x;l=r+1){ // 数论分块
r=x/(x/l);
ret-=calc_phi(x/l)*(r-l+1);
}
return t_phi[x]=ret; // 记忆化
}
int calc_mu(int x){
if(x<=m)return s_mu[x];
if(t_mu.count(x))return t_mu[x];
int ret=1;
for(int l=2,r;l<=x;l=r+1){ // 数论分块
r=x/(x/l);
ret-=calc_mu(x/l)*(r-l+1);
}
return t_mu[x]=ret; // 记忆化
}
void _solve(int n){
ans1=calc_phi(n),ans2=calc_mu(n);
printf("%lld %lld\n",ans1,ans2);
}
signed main(){
T=read();
int mx=0;
for(int i=1;i<=T;i++)_in[i]=read(),mx=max(mx,_in[i]);
m=pow(T*mx,2.0/3);
_init();
for(int i=1;i<=T;i++)_solve(_in[i]);
return 0;
}
0xFF 参考文献
本文参考了一些资料 ^{\textup{\textmd{\textsf{[1-3]}}}}。
[1] 罗勇军, 郭卫斌. 算法竞赛[M]. 北京: 清华大学出版社, 2022: 445-461.
[2] OI Wiki. 数论分块[DB/OL]. (2025-09-23)[2025-10-28]. https://oi-wiki.org/math/number-theory/sqrt-decomposition/.
[3] OI Wiki. 杜教筛[DB/OL]. (2025-10-20)[2025-10-28]. https://oi-wiki.org/math/number-theory/du/.
版权声明：本文中 0x05 数论分块与 0x20 时间复杂度参考自 OI Wiki，其内容在 CC BY-SA 4.0 协议之条款下提供。本文亦采用 CC BY-SA 4.0 协议发布。