什么是莫比乌斯反演
关于莫比乌斯反演
莫比乌斯反演,又称懵逼钨丝繁衍,是一种看了就一脸懵逼的东西。
好吧好吧,严肃点。(不过因为本蒟蒻真的很懵逼,所以错误之处请大神指出)
莫比乌斯反演就是下面这个式子:
如果存在函数F(x)和f(x),满足
F ( n ) = ∑ d ∣ n f ( d ) F(n)=\sum_{d|n} f(d) F(n)=d∣n∑f(d)
那么就有:
f ( n ) = ∑ d ∣ n μ ( d ) F ( n d ) f(n)=\sum_{d|n} \mu(d)F(\frac{n}{d}) f(n)=d∣n∑μ(d)F(dn)
或者如果F(x)和f(x)满足:
F ( n ) = ∑ n ∣ d f ( d ) F(n)=\sum_{n|d}f(d) F(n)=n∣d∑f(d)
那么:
f ( n ) = ∑ n ∣ d μ ( d n ) F ( d ) f(n)=\sum_{n|d}\mu(\frac{d}{n})F(d) f(n)=n∣d∑μ(nd)F(d)
其中 μ ( ) \mu() μ()函数是莫比乌斯函数,定义是:
如果 d = 1 d=1 d=1 , μ ( d ) = 1 \mu(d)=1 μ(d)=1
如果 d d d为互异质数 p 1 p_1 p1, p 2 p_2 p2… p k p_k pk的乘积,则 μ ( d ) = ( − 1 ) k \mu(d)=(-1)^k μ(d)=(−1)k
否则, μ ( d ) = 0 \mu(d)=0 μ(d)=0
所以线性筛莫比乌斯函数往后看,几乎每段代码里都有。
两条性质
好了,现在补充两条 (我不会证,不过脑补一下就好了吧的) 性质:
1.如果 n > 1 n>1 n>1且 n n n为正整数,则有:
∑ d ∣ n μ ( d ) = 0 \sum_{d|n}\mu(d)=0 d∣n∑μ(d)=0
如果 n = 1 n=1 n=1则上式为1.
2.对于任意正整数 n n n均有:
∑ d ∣ n μ ( d ) d = ϕ ( n ) n \sum_{d|n}\frac{\mu(d)}{d}=\frac{\phi(n)}{n} d∣n∑dμ(d)=nϕ(n)
证明莫比乌斯反演
有了上面两条 我不会证的 性质之后,我们就可以证一下莫比乌斯反演了。
由F(x)函数的定义可以得到:
∑ d ∣ n μ ( d ) F ( n d ) = ∑ d ∣ n μ ( d ) ∑ d ′ ∣ n d f ( d ′ ) \sum_{d|n}\mu(d)F(\frac{n}{d})=\sum_{d|n}\mu(d)\sum_{d'|\frac{n}{d}}f(d') d∣n∑μ(d)F(dn)=d∣n∑μ(d)d′∣dn∑f(d′)
接下来我们让 n d = k d ′ \frac{n}{d}=kd' dn=kd′,那么 d = n k d ′ d=\frac{n}{kd'} d=kd′n,所以就有 d ∣ n d ′ d|\frac{n}{d'} d∣d′n,而每个 f ( d ′ ) f(d') f(d′)一定会和一个 μ ( d ) \mu(d) μ(d)乘一次,所以:
∑ d ∣ n μ ( d ) ∑ d ′ ∣ n d f ( d ′ ) = ∑ d ′ ∣ n f ( d ′ ) ∑ d ∣ n d ′ μ ( d ) \sum_{d|n}\mu(d)\sum_{d'|\frac{n}{d}}f(d')=\sum_{d'|n}f(d')\sum_{d|\frac{n}{d'}}\mu(d) d∣n∑μ(d)d′∣dn∑f(d′)=d′∣n∑f(d′)d∣d′n∑μ(d)
再看上面的两条性质里的性质2,明白了吧QWQ!
∑ d ∣ n μ ( d ) F ( n d ) = ∑ d ′ ∣ n f ( d ′ ) ∑ d ∣ n d ′ μ ( d ) = f ( n ) \sum_{d|n}\mu(d)F(\frac{n}{d})=\sum_{d'|n}f(d')\sum_{d|\frac{n}{d'}}\mu(d)=f(n) d∣n∑μ(d)F(dn)=d′∣n∑f(d′)d∣d′n∑μ(d)=f(n)
得证!
莫比乌斯反演的应用
例题:HDU1695
求对于区间[1,b]内的整数x和[1,d]内的y,满足gcd(x,y)=k的数对的个数。
只要让F(t)=满足gcd(x,y)%t==0的数对个数
f(t)=满足gcd(x,y)=t的数对个数,则F(t)和f(t)就存在莫比乌斯反演的关系了。显然F(t)=(b/t)*(d/t)
因为如果 g c d ( x , y ) = 1 gcd(x,y)=1 gcd(x,y)=1,则 g c d ( x ∗ k , y ∗ k ) = k gcd(x*k,y*k)=k gcd(x∗k,y∗k)=k,所以我们把b和d同时除以k,得到的f(1)再去重就是答案。令lim=min(b/k,d/k),而根据上面的莫比乌斯反演第二个公式,
f ( 1 ) = ∑ d l i m μ ( d ) ∗ F ( d ) f(1)=\sum_d^{lim}\mu(d)*F(d) f(1)=d∑limμ(d)∗F(d)
代码:
#include<iostream>
#include<cstdio>
#include<climits>
#include<algorithm>
#include<cstring>
using namespace std;
const int maxn=100005;
int T,a,b,c,d,e,tot;
long long ans1,ans2;
bool is[maxn];
int pri[maxn],miu[maxn];
void init(){//首先把莫比乌斯函数筛出来
miu[1]=1;
for(int i=2;i<=100000;i++){
if(!is[i]){pri[++tot]=i;miu[i]=-1;}
for(int j=1;j<=tot;j++){
int k=pri[j]*i;if(k>100000)break;
is[k]=1;
if(i%pri[j]==0){miu[k]=0;break;}
else miu[k]=-miu[i];
}
}
}
int main()
{
int i,j,cnt=0;
init();
scanf("%d",&T);
while(T--){
cnt++;ans1=ans2=0;
scanf("%d%d%d%d%d",&a,&b,&c,&d,&e);
if(!e){ printf("Case %d: 0\n",cnt);continue;}
b/=e;d/=e;//如果gcd(x,y)=1,那么gcd(x*e,y*e)=e;
if(b>d)swap(b,d);
for(i=1;i<=b;i++)ans1+=(long long)miu[i]*(b/i)*(d/i);
for(i=1;i<=b;i++)ans2+=(long long)miu[i]*(b/i)*(b/i);
printf("Case %d: %lld\n",cnt,ans1-ans2/2);
}
return 0;
}
不过如果你遇上的是洛谷P3455什么的,那就还要用一个类似于分块的公共乘积优化啦,具体看下面的第4条,代码也贴着:
for(i=1;i<=a;i=j+1){
j=min(a/(a/i),b/(b/i));
ans+=(long long)(sum[j]-sum[i-1])*(a/i)*(b/i);
//sum是miu的前缀和
}
例题:bzoj2301/洛谷2522 problem b
这题和上题差不多吧…加上一个容斥就行了。容斥是这样的:
ans=find(b/k,d/k)-find((a-1)/k,d/k)-find((c-1)/k,b/k)+find((a-1)/k,(c-1)/k);
而find函数就是上文代码中的这一段:
for(i=1;i<=min(x,y);i=j+1){
j=min(x/(x/i),y/(y/i));
re+=(long long)(sum[j]-sum[i-1])*(x/i)*(y/i);
}
例题:bzoj2820/洛谷P2257 YY的GCD
题目大意:给定N, M,求1<=x<=N, 1<=y<=M且gcd(x, y)为质数的(x, y)有多少对
假设
n
<
m
n<m
n<m,此处/号代表向下取整
根据两题我们知道,通过
f
(
n
)
=
∑
n
∣
d
μ
(
d
n
)
F
(
d
)
f(n)=\sum_{n|d}\mu(\frac{d}{n})F(d)
f(n)=∑n∣dμ(nd)F(d),然后设
f
(
x
)
f(x)
f(x)为
g
c
d
(
a
,
b
)
=
x
gcd(a,b)=x
gcd(a,b)=x的数对数,
F
(
x
)
F(x)
F(x)为
x
∣
g
c
d
(
a
,
b
)
x|gcd(a,b)
x∣gcd(a,b)的数对数,用莫比乌斯反演第二条可以求出式子:
f
(
x
)
=
∑
x
∣
d
μ
(
d
x
)
n
d
∗
m
d
f(x)=\sum_{x|d}\mu(\frac{d}{x})\frac{n}{d}*\frac{m}{d}
f(x)=x∣d∑μ(xd)dn∗dm
再令上一个式子中的
x
=
p
x=p
x=p,
d
=
p
i
d=pi
d=pi可得:
a
n
s
=
∑
i
s
p
r
i
m
e
(
p
)
n
∑
i
n
(
n
i
p
)
(
m
i
p
)
μ
(
i
)
ans=\sum_{isprime(p)}^n \sum_i^n ( \frac{n}{ip})(\frac{m}{ip})\mu(i)
ans=isprime(p)∑ni∑n(ipn)(ipm)μ(i)
意会一下就可以发现,假如设
T
=
i
p
T=ip
T=ip,那么就有:
a n s = ∑ T n ( n T ) ( m T ) ∑ p ∣ T 且 i s p r i m e ( p ) μ ( T p ) ans=\sum_T^n ( \frac{n}{T})(\frac{m}{T}) \sum_{p|T且isprime(p)} \mu(\frac{T}{p}) ans=T∑n(Tn)(Tm)p∣T且isprime(p)∑μ(pT)
这样我们预处理一下后面那团东西的前缀和即可,就这样。
#include<iostream>
#include<cstdio>
#include<cstdlib>
using namespace std;
#define ll long long
int read(){
int q=0;char ch=' ';
while(ch<'0'||ch>'9')ch=getchar();
while(ch>='0'&&ch<='9')q=q*10+ch-'0',ch=getchar();
return q;
}
const int maxn=10000005;
int n,m,T,tot,lim=10000000;
bool is[maxn];ll sum[maxn];
int pri[maxn>>1],mu[maxn];
void init(){
int i,j,p;ll k;
mu[1]=1;
for(i=2;i<=lim;i++){
if(!is[i]){pri[++tot]=i;mu[i]=-1;}
for(j=1;j<=tot;j++){
k=(ll)i*pri[j];if(k>lim)break;
is[k]=1;
if(i%pri[j])mu[k]=-mu[i];
else break;
}
}
for(i=1;i<=tot;i++)//预处理前缀和的办法
for(j=1;pri[i]*j<=lim;j++)sum[pri[i]*j]+=mu[j];
for(i=1;i<=lim;i++)sum[i]=sum[i-1]+sum[i];
}
int main()
{
init();T=read();
while(T--){
n=read();m=read();
if(n>m)swap(n,m);
ll ans=0;
for(int i=1,j;i<=n;i=j+1){
j=min(n/(n/i),m/(m/i));//一毛一样的优化
ans+=(ll)(n/i)*(ll)(m/i)*(ll)(sum[j]-sum[i-1]);
}
printf("%lld\n",ans);
}
return 0;
}
例题:bzoj2005/洛谷1447 能量采集
dalao们的题解经常长这样:水题,贴代码
水吗…蒟蒻泪流满面。
好吧好吧,那蒟蒻就来跟着dalao推一推。
1.容易意会得到,一颗坐标为(x,y)植物连线上的植物数量就是gcd(x,y)-1;
题目中的式子变一下,主要矛盾是求
∑
x
=
1
n
∑
y
=
1
m
g
c
d
(
x
,
y
)
\sum_{x=1}^n \sum_{y=1}^m gcd(x,y)
x=1∑ny=1∑mgcd(x,y)
2.首先有一条结论:一个数的所有因子的欧拉函数之和等于这个数本身。
然后x和y的公因子肯定是他们最大公约数的因子,所以就是求 ∑ x = 1 n ∑ y = 1 m ∑ d = 1 n [ d ∣ n 且 d ∣ m ] ϕ ( d ) \sum_{x=1}^n \sum_{y=1}^m \sum_{d=1}^n[d|n且d|m] \phi(d) x=1∑ny=1∑md=1∑n[d∣n且d∣m]ϕ(d)
3.随便交换律一下得到: ∑ d = 1 n ∑ i = 1 n [ d ∣ n ] ∑ j = 1 m [ d ∣ m ] ϕ ( d ) ∗ ( n / i ) ∗ ( m / i ) \sum_{d=1}^n \sum_{i=1}^n[d|n] \sum_{j=1}^m[d|m]\phi(d)*(n/i)*(m/i) d=1∑ni=1∑n[d∣n]j=1∑m[d∣m]ϕ(d)∗(n/i)∗(m/i)
注意:n/i和m/i是向下取整,下同
4.现在求起来就方便多了,不过呢,还有一个优化,对于一个数i来说,设它的 n / i = k 1 , m / i = k 2 n/i=k_1,m/i=k_2 n/i=k1,m/i=k2,则同样满足 n / j = k 1 , m / j = k 2 n/j=k_1,m/j=k_2 n/j=k1,m/j=k2的最大 j = m i n ( n / k 1 , m / k 2 ) ; j=min(n/k_1,m/k_2); j=min(n/k1,m/k2);
5.最后把得到的结果乘以2再减去m*n即可
#include<iostream>
#include<cstdio>
#include<climits>
#include<algorithm>
#include<cstring>
using namespace std;
#define ll long long
ll n,m,ans;int tot;
bool is[1000010];
int pri[1000010],phi[10000100];
ll sum[1000010];
void init(){
ll i,j,k;
phi[1]=1;
for(i=2;i<=n;i++){
if(!is[i]){phi[i]=i-1;pri[++tot]=i;}
for(j=1;j<=tot;j++){
k=pri[j]*i;if(k>n)break;
is[k]=1;
if(i%pri[j]==0){phi[k]=(ll)phi[i]*pri[j];break;}
else phi[k]=(ll)(pri[j]-1)*phi[i];
}
}
for(i=1;i<=n;i++)sum[i]=sum[i-1]+phi[i];
}
int main()
{
ll i,j;
scanf("%lld%lld",&n,&m);
if(n>m)swap(n,m);
init();
for(i=1;i<=n;i=j+1){
j=min(m/(m/i),n/(n/i));
ans+=(ll)(sum[j]-sum[i-1])*(m/i)*(n/i);
}
ans=(ll)ans*2-(ll)n*m;
printf("%lld",ans);
return 0;
}
其他练习
相信你经过例题的磨炼,已经对莫比乌斯反演的应用有了基本的了解,下面我们再来几道水题练练手吧。
bzoj2154/bzoj2693/洛谷P1829 Crash的数字表格:题解戳我QvQ
bzoj3994/洛谷P3327 约数个数和:题解戳我QvQ
bzoj3529/洛谷P3312 数表:题解戳我QvQ
洛谷P3768 简单的数学题:题解戳我QvQ
bzoj4816/洛谷P3704/loj2000 数字表格:虽然利用莫比乌斯反演的方法很常规,但是这种相乘的变化还是挺新鲜的。