杜教筛

杜教筛

积性函数

定义

1. 定义在所有正整数上的函数称为算术函数（或数论函数）。
2. 若算术函数 \(f\) 满足对于任意两个互质的正整数 \(p\) 和 \(q\) 均有 \(f(pq)=f(p)f(q)\)，则称该函数为积性函数。
3. 若算术函数 \(f\) 满足对于任意两个正整数 \(p\) 和 \(q\) 均有 \(f(pq)=f(p)f(q)\)，则称该函数为完全积性函数。

常见积性函数

这些积性函数常常用于狄利克雷卷积。

• 元函数，也称单位函数：\(\varepsilon(n)=[n=1]\)。（完全积性函数）
• 恒等函数：\(\operatorname{id}(n)=n\)。（完全积性函数）
• 常数函数：\(1(n)=1\)。（完全积性函数）
• 约数个数函数：\(d(n)=\sum_{d\mid n}1\)。
• 约数和函数：\(\sigma(n)=\sum_{d\mid n}d\)。
• 欧拉函数：\(\varphi(n)\) 表示小于等于 \(n\) 的正整数中与 \(n\) 互质的数的个数。
• 莫比乌斯函数：\(\mu(n)=\begin{cases}1&n=1\\(-1)^m&n=p_1p_2\cdots p_m~,~\text{其中 }p_i\text{ 为不同的素数}\\0&\text{其他} \end{cases}\)。

狄利克雷卷积

定义

两个数论函数的狄利克雷卷积记为 \(f*g\)，如果满足

\[h(n)=\sum_{d\mid n}f(d)g\Big(\frac nd\Big) \]

则记 \(h=f*g\)。

性质

• 交换律：\(f*g=g*f\)。
• 结合律：\((f*g)*h=f*(g*h)\)。
• 分配律：\(f*(g+h)=(f*g)+(f*h)\)。
• 上文提到的 \(\varepsilon\) 函数为狄利克雷卷积中的单位元，即对于任何数论函数 \(f\)，有 \(f*\varepsilon=\varepsilon*f=f\)。类似地，如果两个数论函数 \(f\) 和 \(g\) 满足 \(f*g=\varepsilon\)，则称 \(g\) 是 \(f\) 的逆元，且这个逆元是唯一存在的。
  逆元的作用可以体现为：若 \(f*g=h\)，且 \(f\) 的逆元为 \(q\)，则有 \(g=h*q\)。
• 两个积性函数的狄利克雷卷积仍然是积性函数。

常用狄利克雷卷积

这部分的内容掌握熟练才能熟练运用杜教筛。

• \(\varepsilon=\mu*1\)
  这说明 \(\mu\) 函数和 \(1\) 函数互为逆元。
• \(d=1*1\)，\(\sigma=\text{id}*1\)
  这两个卷积较为容易证明，且由逆元能推导出 \(\mu*d=1\)，\(\mu*\sigma=\text{id}\)。
• \(\text{id}=\varphi*1\)，\(\varphi=\mu*\mathrm{id}\)
  这两个卷积较为重要。第一个式子由 \(n=\sum_{d\mid n}\varphi(d)\) 可得，第二个式子可由逆元得到。

杜教筛

杜教筛用于在低于线性时间内求一些积性函数的前缀和。它运用到了整除分块、狄利克雷卷积、线性筛。

杜教筛公式

\[g(1)S(n)=\sum_{i=1}^nh(i)-\sum_{i=2}^ng(i)S\Big(\Big\lfloor\frac ni\Big\rfloor\Big) \]

下面给出推导。

设 \(f\) 是一个数论函数，我们需要计算 \(S(n)=\sum_{i=1}^nf(i)\)。为了能在低于线性时间内求解，我们尝试运用类似整除分块这类的技巧加速计算。

我们需要构造两个积性函数 \(h\) 和 \(g\)，满足 \(h=g*f\)。具体构造方法基本参考上文给出的常用卷积。

根据卷积的定义有 \(h(i)=\sum_{d\mid i}g(d)f(i/d)\)。对 \(h(i)\) 求和得到

\[\begin{aligned} \sum_{i=1}^nh(i)&=\sum_{i=1}^n\sum_{d\mid i}^ig(d)f\Big(\frac id\Big)\\ &=\sum_{d=1}^ng(d)\sum_{d\mid i}^nf\Big(\frac id\Big)\\ &=\sum_{d=1}^ng(d)\sum_{i=1}^{\lfloor n/d\rfloor}f(i)\\ &=\sum_{d=1}^ng(d)S\Big(\Big\lfloor\frac nd\Big\rfloor\Big) \end{aligned} \]

注意第二行的变换比较特殊，称为杜教筛变换，具体证明略。

得到最后一行后，把等式右边的第一项提出得到

\[\begin{matrix} \displaystyle\sum_{i=1}^nh(i)=g(1)S(n)+\sum_{d=2}^ng(d)S\Big(\Big\lfloor\frac nd\Big\rfloor\Big)\\ \displaystyle g(1)S(n)=\sum_{i=1}^nh(i)-\sum_{i=2}^ng(d)S\Big(\Big\lfloor\frac nd\Big\rfloor\Big) \end{matrix} \]

这就是杜教筛公式。

复杂度

经过复杂的推导，复杂度是 \(O(n^{2/3})\) 的水平。

应用

一般我们构造的 \(g\) 和 \(h\) 都应是容易计算前缀和的函数。具体实现上采用了记忆化搜索。

比如在模板题中，需要计算 \(\sum\varphi(i)\) 和 \(\sum\mu(i)\)。前者运用 \(\text{id}=\varphi*1\) 即可，后者运用 \(\varepsilon=\mu*1\) 即可。

#include<bits/stdc++.h>
using namespace std;

using ll=long long;
constexpr int MAXN=2e6+5;
bool isp[MAXN];
int mu[MAXN];
ll phi[MAXN];
vector<int>pri;
unordered_map<int,ll>smu,sphi;

void init(int n){
	mu[1]=phi[1]=1;
	for(int i=2;i<=n;i++){
		if(!isp[i]) pri.emplace_back(i),mu[i]=-1,phi[i]=i-1;
		for(auto j:pri){
			if(i>n/j) break;
			isp[i*j]=1;
			if(i%j==0){
				mu[i*j]=0;
				phi[i*j]=phi[i]*j;
				break;
			}
			mu[i*j]=-mu[i];
			phi[i*j]=phi[i]*phi[j];
		}
	}
	for(int i=1;i<=n;i++) mu[i]+=mu[i-1],phi[i]+=phi[i-1];
}
ll gphi(ll x){
	if(x<MAXN) return phi[x];
	if(sphi.count(x)) return sphi[x];
	ll ans=1ll*x*(x+1)>>1;
	for(ll l=2,r;l<=x;l=r+1){
		r=x/(x/l);
		ans-=(r-l+1)*gphi(x/l);
	}
	return sphi[x]=ans;
}
ll gmu(ll x){
	if(x<MAXN) return mu[x];
	if(smu.count(x)) return smu[x];
	ll ans=1;
	for(ll l=2,r;l<=x;l=r+1){
		r=x/(x/l);
		ans-=(r-l+1)*gmu(x/l);
	}
	return smu[x]=ans;
}

int main(){
	init(2e6);
	int T,n;
	scanf("%d",&T);
	while(T--){
		scanf("%d",&n);
		printf("%lld %lld\n",gphi(n),gmu(n));
	}
	return 0;
}

其实不用两个杜教筛就能求解。

注意由 \(\varphi=\mu*\mathrm{id}\) 可得

\[\begin{aligned}\sum_{i=1}^n\varphi(i)&=\sum_{i=1}^n\sum_{d\mid i}d\mu\Big(\frac id\Big)\\&=\sum_{d=1}^nd\sum_{i=1}^{\lfloor n/d\rfloor}\mu(i)\end{aligned} \]

所以只需对 \(\mu\) 使用杜教筛，然后 \(\varphi\) 用整除分块直接求解即可。

#include<bits/stdc++.h>
using namespace std;

using ll=long long;
constexpr int MAXN=2e6+5;
bool isp[MAXN];
int mu[MAXN];
ll phi[MAXN];
vector<int>pri;
unordered_map<int,ll>smu;

void init(int n){
	mu[1]=1;
	for(int i=2;i<=n;i++){
		if(!isp[i]) pri.emplace_back(i),mu[i]=-1;
		for(auto j:pri){
			if(i>n/j) break;
			isp[i*j]=1;
			if(i%j==0) break;
			mu[i*j]=-mu[i];
		}
	}
	for(int i=1;i<=n;i++) mu[i]+=mu[i-1];
}
ll gmu(ll x){
	if(x<MAXN) return mu[x];
	if(smu.count(x)) return smu[x];
	ll ans=1;
	for(ll l=2,r;l<=x;l=r+1){
		r=x/(x/l);
		ans-=(r-l+1)*gmu(x/l);
	}
	return smu[x]=ans;
}
ll gphi(ll x){
	ll ans=0;
	for(ll l=1,r;l<=x;l=r+1){
		r=x/(x/l);
		ans+=(l+r)*(r-l+1)*gmu(x/l)>>1;
	}
	return ans;
}

int main(){
	init(2e6);
	int T,n;
	scanf("%d",&T);
	while(T--){
		scanf("%d",&n);
		printf("%lld %lld\n",gphi(n),gmu(n));
	}
	return 0;
}