莫比乌斯反演学习笔记

嘛......开个坑吧 咕不咕就不知道了qwq......

数论分块

这个东西之前水过一篇来着：数论分块

狄利克雷卷积

定义

对于两个数论函数(f(n),g(n))，定义他们的狄利克雷卷积（Dirichlet卷积）为

[(f*g)(n)=sumlimits_{d mid n} f(d)g(frac{n}{d}) ]

注意柿子中(d,frac{n}{d})是对称的，所以

[(f*g)(n)=(g*f)(n)=sumlimits_{dmid n} f(frac{n}{d})g(d) ]

也就是说Dirichlet卷积满足交换律，同时狄利克雷卷积也满足结合律，证明的话暴力把柿子写开就好了，这里就不写了qaq

积性函数

对于一个数论函数(f(n))，如果对于两个互质的整数(a,b)，有(f(ab)=f(a)f(b))，则称(f(n))为积性函数。

如果对于所有的(a,b)都有(f(ab)=f(a)f(b))，那么(f(n))就是一个完全积性函数。

特别的对于两个积性函数(f(n),g(n))，他们的Dirichlet卷积也是积性函数。

常见的积性函数

(varepsilon(n)=[n=1])，这个函数很重要，它也是Dirichlet卷积的单位元（所有函数卷它都是原来的函数，类比乘法中的(1)）。

(1(n)equiv 1)，常数函数，它永远等于(1)。

(id_k(n)=n^k)，恒等函数，一般直接将(id_1(n))记作(id(n))。

(sigma_k(n)=sumlimits_{dmid n}d^k)，除数函数，(sigma_0(n))就是除数个数，一般写作( au(n))或(d(n))，(sigma_1(n))即约数之和，记作(sigma(n))。

(varphi(n))，Euler函数，即(1cdots n-1)内与(n)互质的数的个数（同时也只个积性函数），同时也是个积性函数。

(mu(n))，即莫比乌斯函数，可以理解为Dirichlet卷积中(1)的逆元（类比乘法逆元），也就是说((mu*1)(n)=sumlimits_{dmid n} mu(d)=varepsilon(n)=[n=1])。（定义在下面会讲）

以及他们的狄利克雷卷积的形式

[varepsilon=mu * 1 : varepsilon(n)=sumlimits_{dmid n} mu(d) ]

[d = 1*1 : d(n)=sumlimits_{dmid n} 1 ]

[sigma=id*1 : sigma(n)=sumlimits_{dmid n} d ]

---

莫比乌斯函数

终于讲到点东西了...

上面我们知道(sumlimits_{dmid n} mu(d)=[n=1])，那么(mu)的定义到底是啥呢？

先把(n)质因数分解

[n=prod_{i=1}^k p_i^{c_i} ]

那么

[mu(n)=egin{cases} (-1)^k quad ext{所有}c_i<2 \ 0 quadquad ext{存在一个}c_i ge2 end{cases} ]

鬼知道它为什么长这样......

特别的(mu(1)=1)

如果(n)有一个平方因子的话(mu(n))就是(0)，不然就是(-1)的质因数个数次方。

首先对于(n=1)，有(sumlimits_{dmid n} mu(d)=1)

那么其他情况呢？

令(n'=prod_{i=1}^k p_i)，显然有

[sumlimits_{d mid n}mu(d)=sumlimits_{dmid n'} mu(d) ]

考虑(sumlimits_{d mid n'} mu(d))，枚举质因子的个数，有

[sumlimits_{d mid n'} mu(d) = sumlimits_{i=0}^k (-1)^iinom{k}{i}=(1-1)^k ]

显然是(0)。

故(sumlimits_{dmid n} mu(d) =[n=1])

我们也得到了反演结论：

[varepsilon(n)=[n=1]=sumlimits_{dmid n} mu(d) ]

线性筛(mu)

(mu)是一个积性函数，所以我们考虑在线性筛的过程中算出(mu)（不知道的话可以去学一学线性筛qwq）。有如下代码

void getmu() {
	mu[1]=1;
	rep(i,2,n+1) {
		if(!vis[i]) {mu[i]=-1;p[pn++]=i;}
		for(int j=0;j<pn&&i*p[j]<=n;j++) {
			vis[i*p[j]]=1;
			if (i%p[j]) mu[i*p[j]]=mu[i]*mu[p[j]]; // 这里也可以写成mu[i*p[j]]=-mu[i]
			else {mu[i*p[j]]=0;break;}
		}
	}
}

莫比乌斯反演

对于两个数论函数(f,g)

有

[f(n)=sumlimits_{dmid n} g(d) Leftrightarrow g(n)=sumlimits_{dmid n} mu(d)f(frac{n}{d}) ]

证明的话可以考虑Dirichlet卷积

[f=g*1 ]

两边同时乘(mu)

[f*mu=g*(1*mu)=g*varepsilon=g ]

所以

[f=g*1 Leftrightarrow g=f*mu ]

其实反演的柿子并不会怎么用到2333

---

例题

还是来看题吧

LuoguP3455 [POI2007]ZAP-Queries

给(n,m,k)，让你求

[sumlimits_{i=1}^nsumlimits_{j=1}^m [gcd(i,j)=k] ]

多组数据。

不难发现这个柿子等价于

[sumlimits_{i=1}^{leftlfloorfrac{n}{k} ight floor}sumlimits_{j=1}^{leftlfloorfrac{m}{k} ight floor}[gcd(i,j)=1] ]

把这些两两互质的同时乘上(k)就是(gcd=k)的数。

那么简化问题，令

[S(n,m)=sumlimits_{i=1}^nsumlimits_{j=1}^m[gcd(i,j)=1] ]

答案就是(S(leftlfloorfrac{n}{k} ight floor,leftlfloorfrac{m}{k} ight floor))。

发现后面中括号就是一个单位函数，于是

[S(n,m)=sumlimits_{i=1}^nsumlimits_{j=1}^m sumlimits_{dmid gcd(i,j)} mu(d) ]

考虑把(d)提到外面

[sumlimits_{d} mu(d) sumlimits_{i=1}^{n}sumlimits_{j=1}^{m} [dmid i and dmid j] ]

后面两个(Sigma)就是(leftlfloorfrac{n}{d} ight floorleftlfloorfrac{m}{d} ight floor)。

所以(S(n,m)=sumlimits_{d} mu(d)leftlfloorfrac{n}{d} ight floorleftlfloorfrac{m}{d} ight floor)，(d)最大显然是到(min(n,m))

线性筛出莫比乌斯函数预处理出前缀和，数论分块就好了。

复杂度(O(Tsqrt{n}))（(T)为数据组数）

#include <bits/stdc++.h>
using namespace std;
#define rep(i,a,b) for (int i=(a);i<(b);++i)
#define per(i,a,b) for (int i=(a)-1;i>=(b);--i)
#define mp make_pair
#define pb push_back
#define fi first
#define se second
#define all(x) (x).begin(),(x).end()
#define SZ(x) ((int)(x).size())
typedef double db;
typedef long long ll;
typedef pair<int,int> PII;
typedef vector<int> VI;

const int N=5e4+10;
int mu[N],sum[N],vis[N];
int p[N],pn;

void init(int n) {
	mu[1]=1;
	rep(i,2,n+1) {
		if(!vis[i]) {
			mu[i]=-1;
			p[pn++]=i;
		}
		for(int j=0;j<pn&&i*p[j]<=n;j++) {
			vis[i*p[j]]=1;
			if (i%p[j]) mu[i*p[j]]=mu[i]*mu[p[j]];
			else {mu[i*p[j]]=0;break;}
		}
	}
	rep(i,1,n+1) sum[i]=sum[i-1]+mu[i];
}
#define S(l,r) (sum[r]-sum[l-1])

int solve(int n,int m,int d) {
	n/=d,m/=d; int tn=min(n,m);
	ll ans=0;
	for(int l=1,r=0;l<=tn;l=r+1) {
		r=min(n/(n/l),m/(m/l));
		ans+=1ll*S(l,r)*(n/l)*(m/l);
	}
	return ans;
}

int main() {
	init(50000);
	int _,n,m,d; for(scanf("%d",&_);_;_--) {
		scanf("%d%d%d",&n,&m,&d);
		printf("%d
",solve(n,m,d));
	}
	return 0;
}

LuoguP2522 [HAOI2011]Problem b

这个就是上面那题吧......类似二维前缀和容斥一下就没了，代码就不放了qwq实际上是为了提示双倍经验

LuoguP1829 [国家集训队]Crash的数字表格 / JZPTAB

给(n,m)，让你求

[sumlimits_{i=1}^n sumlimits_{j=1}^m operatorname{lcm}(i,j) ]

首先(operatorname{lcm}(i,j)=frac{ij}{gcd(i,j)})

那么就是

[sumlimits_{i=1}^nsumlimits_{j=1}^m frac{ij}{gcd(i,j)} ]

枚举(gcd)

[sumlimits_{i=1}^nsumlimits_{j=1}^m sumlimits_{dmid i,dmid j} [gcd(frac{i}{d},frac{j}{d})=1]frac{ij}{d} ]

把(d)提到外面，然后套路的变一下

[sumlimits_{d} sumlimits_{i=1}^{leftlfloorfrac{n}{d} ight floor}sumlimits_{j=1}^{leftlfloorfrac{m}{d} ight floor} [gcd(i,j)=1] frac{id imes jd}{d} ]

[= sumlimits_{d}d sumlimits_{i=1}^{leftlfloorfrac{n}{d} ight floor}sumlimits_{j=1}^{leftlfloorfrac{m}{d} ight floor} [gcd(i,j)=1] ij ]

下面令(S(n,m)=sumlimits_{i=1}^nsumlimits_{j=1}^m [gcd(i,j)=1] ij)

则(Ans=sumlimits_{d} d imes S(leftlfloorfrac{n}{d} ight floor,leftlfloorfrac{m}{d} ight floor))

我们想快一点求出(S)，然后对(d)数论分块。

推柿子

[egin{aligned} S(n,m)&=sumlimits_{i=1}^nsumlimits_{j=1}^msumlimits_{dmid gcd(i,j)} mu(d)ij \ &=sumlimits_{d} mu(d)sumlimits_{i=1}^{leftlfloorfrac{n}{d} ight floor}idsumlimits_{j=1}^{leftlfloorfrac{m}{d} ight floor}jd \ &=sumlimits_{d} mu(d)d^2sumlimits_{i=1}^{leftlfloorfrac{n}{d} ight floor}isumlimits_{j=1}^{leftlfloorfrac{m}{d} ight floor}j end{aligned} ]

后面的两个(Sigma)就直接等差数列求和吧，然后预处理出(mu(d)d^2)的前缀和，数论分块就好了。

我们能在根号的时间内求出(S)，求(Ans)的话也用数论分块就好了。

这样就是外层数论分块套内层数论分块。复杂度不太清楚(qwq)，好像上界是(O(n))的，总之跑的过就是了。

#include <bits/stdc++.h>
using namespace std;
#define rep(i,a,b) for (int i=(a);i<(b);++i)
#define per(i,a,b) for (int i=(a)-1;i>=(b);--i)
#define mp make_pair
#define pb push_back
#define fi first
#define se second
#define all(x) (x).begin(),(x).end()
#define SZ(x) ((int)(x).size())
typedef double db;
typedef long long ll;
typedef pair<int,int> PII;
typedef vector<int> VI;

const int maxn=1e7,N=maxn+10,P=20101009;
inline int add(int x,int y) {return (x+=y)>=P?x-P:x;}
inline int sub(int x,int y) {return (x-=y)<0?x+P:x;}
inline int normal(int x) {return x<0?x+P:x;}
inline int sqr(int x) {return 1ll*x*x%P;}
inline int fpow(int x,int y) {
	int ret=1; for(;y;y>>=1,x=1ll*x*x%P)
		if(y&1) ret=1ll*ret*x%P;
	return ret;
}
const int i2=(P+1)/2,i6=fpow(6,P-2);

int mu[N],sum[N],vis[N],p[N],pn;

void init(int n) {
	mu[1]=1;
	rep(i,2,n+1) {
		if(!vis[i]) {p[pn++]=i;mu[i]=-1;}
		for(int j=0;j<pn&&p[j]*i<=n;j++) {
			vis[i*p[j]]=1;
			if(i%p[j]==0) {mu[i*p[j]]=0;break;}
			else mu[i*p[j]]=mu[i]*mu[p[j]];
		}
	}
	rep(i,1,n+1) mu[i]=normal(mu[i]);
	rep(i,1,n+1) sum[i]=add(sum[i-1],1ll*mu[i]*sqr(i)%P);
}

inline int s1(int l,int r) {return 1ll*(r-l+1)*(l+r)%P*i2%P;}
inline int ss(int l,int r) {return sub(sum[r],sum[l-1]);}

int S(int n,int m) {
	int tn=min(n,m),ans=0;
	for(int l=1,r=0;l<=tn;l=r+1) {
		r=min(n/(n/l),m/(m/l));
		ans=add(ans,1ll*ss(l,r)*s1(1,n/l)%P*s1(1,m/l)%P);
	}
	return ans;
}

int solve(int n,int m) {
	int ans=0,tn=min(n,m);
	for(int l=1,r=0;l<=tn;l=r+1) {
		r=min(n/(n/l),m/(m/l));
		ans=add(ans,1ll*s1(l,r)*S(n/l,m/l)%P);
	}
	return ans;
}

int main() {
	init(maxn);
	int n,m;scanf("%d%d",&n,&m);
	printf("%d
",solve(n,m));
	return 0;
}

LuoguP3704 [SDOI2017]数字表格

让你求

[prodlimits_{i=1}^nprodlimits_{j=1}^m f_{gcd(i,j)} ]

多组数据

其中

[f_i=egin{cases}0 quadquadquadquadquad (n=0)\ 1 quadquadquadquadquad (n=1) \ f_{i-1}+f_{i-2} quad (n>1)end{cases} ]

第一次以为是Fibonacci然后调了半天样例怎么都过不去qwq......

枚举柿子中的(gcd)

[prod_{k} f_k^{sum_{i=1}^nsum_{j=1}^m [gcd(i,j)=k]} ]

发现上面的柿子我们在第一题中就推过了2333

令

[S(n,m)=sumlimits_{i=1}^nsumlimits_{j=1}^m [gcd(i,j)=1]=sumlimits_{d} mu(d)leftlfloorfrac{n}{d} ight floorleftlfloorfrac{m}{d} ight floor ]

原来的柿子就是

[prod_{k} f_k^{S(leftlfloorfrac{n}{k} ight floor,leftlfloorfrac{n}{k} ight floor)} ]

这个柿子我们已经可以数论分块套数论分块做到(O(n))的复杂度了，但是它还多测qaq，所以我们再推一下柿子

注意到(leftlfloorfrac{leftlfloorfrac{n}{a} ight floor}{b} ight floor=leftlfloorfrac{n}{ab} ight floor)（感性理解一下先每(a)个数中选一个出来再在选出的数中每(b)个中选一个，等价于在原来(n)个数中每(ab)个里面选一个出来）

[egin{aligned}Ans=prodlimits_{k} f_k^{sum_{d}mu(d)leftlfloorfrac{n}{dk} ight floorleftlfloorfrac{m}{dk} ight floor}end{aligned} ]

枚举(i=dk)，有

[Ans=prod_{i=1}^n left(prod_{j mid i}f_j^{mu(i/j)} ight)^{leftlfloorfrac{n}{i} ight floorleftlfloorfrac{m}{i} ight floor} ]

指数上面的该数论分块还是数论分块，然后发现括号内的可以在调和级数复杂度内预处理出前缀的乘积，要取一段区间的乘积的话乘上逆元就好了。复杂度(O(n ln n +Tsqrt{n}))。

#include <bits/stdc++.h>
using namespace std;
#define rep(i,a,b) for (int i=(a);i<(b);++i)
#define per(i,a,b) for (int i=(a)-1;i>=(b);--i)
#define mp make_pair
#define pb push_back
#define fi first
#define se second
#define all(x) (x).begin(),(x).end()
#define SZ(x) ((int)(x).size())
typedef double db;
typedef long long ll;
typedef pair<int,int> PII;
typedef vector<int> VI;

const int N=1e6+10,maxn=1e6,P=1e9+7;
inline int add(int x,int y) {return (x+=y)>=P?x-P:x;}
inline int sub(int x,int y) {return (x-=y)<0?x+P:x;}
inline int fpow(int x,int y) {
	int ret=1; for(;y;y>>=1,x=1ll*x*x%P)
		if(y&1) ret=1ll*ret*x%P;
	return ret;
}
inline int normal(int x) {return x<0?x+P:x;}

int vis[N],p[N],pn,mu[N];
void getmu(int n) {
	mu[1]=1;
	rep(i,2,n+1) {
		if(!vis[i]) {mu[i]=-1;p[pn++]=i;}
		for(int j=0;j<pn&&i*p[j]<=n;j++) {
			vis[i*p[j]]=1;
			if(i%p[j]==0) {mu[i*p[j]]=0;break;}
			else mu[i*p[j]]=-mu[i];
		}
	}
}

int fib[N],f[N],ifib[N];
void init(int n) {
	getmu(n);
	fib[0]=0,fib[1]=1;
	rep(i,2,n+1) fib[i]=add(fib[i-1],fib[i-2]);	
	rep(i,0,n+1) ifib[i]=fpow(fib[i],P-2),f[i]=1;
	rep(i,1,n+1) if(mu[i]!=0) {
		for(int j=i;j<=n;j+=i)
			f[j]=1ll*f[j]*(mu[i]==-1?ifib[j/i]:fib[j/i])%P;
	}
	rep(i,1,n+1) f[i]=1ll*f[i]*f[i-1]%P; 
}

int solve(int n,int m) {
	int tn=min(n,m),ans=1;
	for(int l=1,r=0;l<=tn;l=r+1) {
		r=min(n/(n/l),m/(m/l));
		ans=1ll*ans*fpow(1ll*f[r]*fpow(f[l-1],P-2)%P,1ll*(n/l)*(m/l)%(P-1))%P;
	}
	return ans;
}

int main() {
#ifdef LOCAL
	freopen("a.in","r",stdin);
#endif
	init(maxn);
	int n,m,_; for(scanf("%d",&_);_;_--) {
		scanf("%d%d",&n,&m);
		printf("%d
",solve(n,m));
	}
	return 0;
}

LuoguP2257 YY的GCD

蒟蒻的题解