莫比乌斯反演与杜教筛

莫比乌斯反演与杜教筛

积性函数

定义

对于一个数论函数 \(f\)，若满足 \(\forall (a,b)=1\) 都有 \(f(ab)=f(a)f(b)\)，那么称 \(f\) 为积性函数。

例子

常见的积性函数有很多，比如以下：

• \(id(x)=x\)。
• \(\epsilon(x)=[x=1]\)。
• \(d(x)=\sum_{i=1}[i|x]\)。

狄利克雷卷积

对于两个数论函数 \(f,g\) 定义他们的狄利克雷卷积

\[f*g(x)=\sum_{ij=x} f_ig_j \]

这个玩意有个漂亮的性质，我们后面再谈。

【例题】

求 \(1*\varphi\)。

解答：原式相当于 \(\sum_{i|x} \varphi(i)\)，下面介绍两种推法：

法1：展开法 考虑 \(x\) 的标准形式 \(x=p_1^{\alpha _1}p_2^{\alpha_2}\dots p_k^{\alpha_k}\)，那么观察 \(\varphi(x)\) 的计算式，发现它 \(p_i^{\alpha_i}\) 的系数都在一下的集合中挑选：

\[\{1,p-1,(p-1)p,\dots,(p-1)p^{\alpha-1}\} \]

不难发现上述集合所有元素的和恰好为 \(p^{\alpha}\)，全部乘起来结果恰为 \(n\)。

法2：一一对应法 考虑所有 \(\leq x\) 的正整数 \(n\)。那么 \(\dfrac{n}{\gcd(n,x)}\) 必然和 \(\dfrac{x}{\gcd(n,x)}\) 互质，这就形成了一个一一对应，所以 \(1*\varphi=id\)。

性质

【性质1】 只需要知道所有 \(x\) 为质数的方幂时的值就可以确定整个函数，换句话说假设 \(x\) 的标准分解为 \(p_1^{\alpha _1}p_2^{\alpha_2}\dots p_k^{\alpha_k}\) 就会有：

\[f(x)=\prod f(p_i^{\alpha_i}) \]

【性质2】 积性函数的狄利克雷卷积为积性函数。

证明：考虑所有的 \(x\) 为质数方幂 \(p^a\) 的值，\(f*g(p^a)=\sum_{i+j=a} f(p^i)g(p^j)\)。不难发现其他的 \(x\) 也可以可以写成标准分解的形式，证毕。

莫比乌斯函数 \(\mu\)

我们需要寻找一个函数 \(\mu\) 使得 \(1*\mu=\epsilon\)。

• 因为 \(1\) 为积性函数，所以 \(\mu\) 为积性函数。
• \(1*\mu(1)=1\)，所以 \(\mu(1)=1\)。
• 假设 \(p\) 为质数，那么 \(1*\mu(p)=\mu(1)+\mu(p)=0\)，所以 \(\mu(p)=-1\)。
• 假设 \(p\) 为质数，\(k\geq 2\)。那么 \(1*\mu(p)=\mu(1)+\mu(p)+\dots +\mu(p^k)\)，不难得到 \(\mu(p^k)=0\)。
• 有了前面的条件可以知道：

\[\mu(x)=\begin{cases}0\ (∃i>1,i^2|x)\\(-1)^k\ (i\text{ have } k \text{ different primeson.})\end{cases} \]

给一个线性筛的 \(\mu\) 求法：

mu[1]=1;
for(int i=2;i<MAXN;i++){
	if(!vis[i]) mu[i]=-1,p[++cnt]=i;
	for(int j=1;j<=cnt&&p[j]*i<MAXN;j++){
		vis[i*p[j]]=1;
		if(i%p[j]==0) break;
		mu[i*p[j]]=-mu[i];
	}
}

莫比乌斯反演

简介

假设我们需要函数 \(f\) 不好求，但是 \(1*f=g\) 是好求的函数，我们就可以有 \(\mu*1*f=g*\mu\) 所以 \(\epsilon*f=g*\mu\)，也就是 \(f(1)=g*\mu\)。\(f(2,3,\dots)\) 是同理的，直接暴力复杂度是 \(O(n\ln n)\) 的。

【习题】补全代码

• 给定数论函数 \(g\) 求 \(f\) 使得 \(1*f=g\)。
• \(n\leq 10^6\)。

应用

【例题1】P3455 [POI2007]ZAP-Queries 改编

• 求 \(\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=1]\)。
• 多测，\(T,n,m\leq 5\times 10^4\)。

对 \(\gcd=1\) 问题一般使用莫反处理

考虑如下变换：

\[\begin{aligned}\sum_{i=1}^a\sum_{j=1}^b[\gcd(i,j)=1]&=\sum_{i=1}^a\sum_{j=1}^b\sum_{k|\gcd(i,j)}\mu(k)\\&=\sum_{k=1}\sum_{k|i}^a\sum_{k|j}^b\mu(k)\\&=\sum_{k=1}\ \mu(k)\times [\dfrac ak]\times[\dfrac bk]\end{aligned} \]

这里第一步 \(\epsilon=1*\mu\) 的性质。对最后的式子数论分块可以做到 \(O(T\sqrt n)\)。

ll solve(){
    ll ans=0;
	for(int l=1,r;l<=min(a,b);l=r+1){
		r=min(a/(a/l),b/(b/l));
		ans+=(a/l)*(b/l)*(mu[r]-mu[l-1]);
	}
    return ans;
}

【例题2】P2257 YY的GCD

• 求 \(\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)\in \textbf P]\)。
• 多测，\(T\leq 10^4,n,m\leq10^7\)。

我们还是像上一题一样向枚举 \(\gcd\) 得到这个奇怪式子：

\[\sum_{p\in\textbf P}\sum_{k=1}\ \mu(k)\times [\dfrac a{pk}]\times[\dfrac b{pk}] \]

这个时间复杂度是 \(T\sqrt n\ln n\) 的，不可接受，考虑简化。

我们可以尝试枚举 \(pk\)，这样子就变成了：

\[\begin{aligned} \sum_{T}\sum_{pk=T} \mu(k)\times [\dfrac aT]\times [\dfrac bT]&=\sum_T\sum_{k|T} \mu(k)\times [\dfrac aT]\times [\dfrac bT] \end{aligned} \]

不难发现我们可以把 \(\sum_{k|T} \mu(k)\) 提前预处理出来，最后直接乘起来就好了，时间复杂度是 \(O(T\sqrt n)\) 可以接受。

#include<iostream>
using namespace std;
#define ll long long
const int MAXN=1e7+5;
ll mu[MAXN],pr[MAXN],cnt=0;bool vis[MAXN];
ll f[MAXN];
void init(){
	mu[1]=1;
	for(int i=2;i<MAXN;i++){
		if(!vis[i]) pr[++cnt]=i,mu[i]=-1;
		for(int j=1;j<=cnt&&i*pr[j]<MAXN;j++){
			vis[i*pr[j]]=1;
			if(i%pr[j]==0) break;
			mu[i*pr[j]]=-mu[i];
		}
	}
	for(int i=1;i<MAXN;i++)
		for(int j=1;j<=cnt&&i*pr[j]<MAXN;j++)
			f[i*pr[j]]+=mu[i];
	for(int i=1;i<MAXN;i++)
		f[i]+=f[i-1];
}
void solve(){
	ll a,b;cin>>a>>b;
	ll ans=0;
	for(int l=1,r;l<=min(a,b);l=r+1){
		r=min(a/(a/l),b/(b/l));
		ans+=(a/l)*(b/l)*(f[r]-f[l-1]);
	}
	cout<<ans<<endl;
	return;
}
int main(){
	init();
	int _;cin>>_;
	while(_--) solve();
}

【习题】P2522 [HAOI2011]Problem b

【习题】P1829 [国家集训队]Crash的数字表格 / JZPTAB

【例题3】P3327 [SDOI2015]约数个数和

• 求 \(\sum_{i=1}^n\sum_{j=1}^md(ij)\)。
• 多测，\(T,n,m\leq5\times 10^4\)。

为了解决本题需要给一个有力的结论把 \(d\) 函数转化成 \(\gcd\) 函数：

\[d(xy)=\sum_{i|x}\sum_{j|y} [\gcd(i,j)=1] \]

证明留给读者，一般来说解决数论等式的证明问题有两种手法，一种是直接展开，一种是构造一一映射。

后面的推导留给读者。

【习题】AT_agc038_c [AGC038C] LCMs

【习题】P3312 [SDOI2014]数表

【习题】P3704 [SDOI2017]数字表格

杜教筛

我们需要解决数论函数快速求前缀和问题，我们需要得到某一项的前缀和，伟大的杜老师发明了 \(O(n^{0.66})\) 筛出来的方法。

简单来说，对于给定函数 \(f\)，我们需要找一个函数 \(g\) 使得 \(g,f*g\) 的前缀和都能 \(O(1)\)。那么考虑如下变换：

\[\sum_{i=1}^n f*g(i)=\sum_{i=1}^n\sum_{xy=i} f(x)g(y)=\sum_{i=1}^n (\sum_{j=1}^if(i))g([\dfrac ni]) \]

假设 \(f\) 的前缀函数叫 \(S\)，那么不难

\[\sum_{i=1}^n f*g(i)=g(1)S(n)+\sum_{i=2}^n S([\dfrac ni])g(i) \]

分治，等式左边我们知道，最右边那一坨可以分治求，解决。时间复杂度：\(O(n^{0.66})\)。具体来说，我们要设置一个阈值 \(X\) 使得 \(\leq X\) 的前缀和暴力求，\(>X\) 的用 unordered_map 记录维护即可。

【例题】P4213 【模板】杜教筛（Sum）

• 求 \(\varphi,\mu\) 的前缀和。
• 多测，\(T\leq 10,n\leq 2^{31}\)。

不难发现：

• \(\varphi*1=id\)（上面例题）
• \(\mu *1=\epsilon\)（\(\mu\) 定义）

套公式即可。

const int MAXN=1e7+5;
ll phi[MAXN],mu[MAXN];bool vis[MAXN];
int p[MAXN],cnt=0;
int tmp=0;
unordered_map<int,ll> M1,M2;
void init(){
	phi[1]=mu[1]=1;
	for(int i=2;i<MAXN;i++){
		if(!vis[i]) p[++cnt]=i,phi[i]=i-1,mu[i]=-1;
		for(int j=1;j<=cnt&&i*p[j]<MAXN;j++){
			vis[i*p[j]]=1;
			if(i%p[j]==0){
				phi[i*p[j]]=phi[i]*p[j];
				break;
			}
			phi[i*p[j]]=phi[i]*(p[j]-1);
			mu[i*p[j]]=-mu[i];
		}
	}
	for(int i=1;i<MAXN;i++)
		phi[i]+=phi[i-1],mu[i]+=mu[i-1];
}
ll pre_phi(ll x){
	if(x<MAXN) return phi[x];
	if(M1[x]) return M1[x];
	ll ans=x*(x+1)/2;
	for(ll l=2,r;l<=x;l=r+1){
		r=x/(x/l);
		ans-=(r-l+1)*pre_phi(x/l);
	}
	return M1[x]=ans;
}
ll pre_mu(ll x){
	if(x<MAXN) return mu[x];
	if(M2[x]) return M2[x];
	ll ans=1;
	for(ll l=2,r;l<=x;l=r+1){
		r=x/(x/l);
		ans-=(r-l+1)*pre_mu(x/l);
	}
	return M2[x]=ans;
}
int main(){
	init();
	int _;cin>>_;
	while(_--){
		ll n;cin>>n;
		cout<<pre_phi(n)<<" "<<pre_mu(n)<<endl;
	} 
}

【习题】P3768 简单的数学题