1 简介
莫比乌斯反演,俗称懵逼钨丝繁衍,顾名思义就是一看就叫人懵逼的算法。我在颓废了一天的莫反之后,成功地被莫反弄懵逼了。它的根本思想就是容斥原理,而莫比乌斯函数则用于辅助进行容斥的计算。
莫比乌斯函数的定义像一个容斥符号的定义:
对于
n
=
p
1
p
2
p
3
…
p
k
n=p_1p_2p_3…p_k
n=p1p2p3…pk,有
μ
(
n
)
=
(
−
1
)
k
\mu(n)=(-1)^k
μ(n)=(−1)k,否则
μ
(
n
)
=
0
\mu(n)=0
μ(n)=0
如果对于一个函数
f
(
x
)
f(x)
f(x),若有另一个函数
g
(
n
)
=
∑
d
∣
n
f
(
d
)
g(n)=\sum_{d|n}f(d)
g(n)=d∣n∑f(d)
那么就有
f
(
n
)
=
∑
d
∣
n
μ
(
d
)
g
(
n
d
)
f(n)=\sum_{d|n}\mu(d)g(\frac n d)
f(n)=d∣n∑μ(d)g(dn)
另一种反演形式:
g
(
n
)
=
∑
n
∣
d
f
(
d
)
⇒
f
(
n
)
=
∑
n
∣
d
μ
(
d
n
)
g
(
d
)
g(n)=\sum_{n|d}f(d)\Rightarrow f(n)=\sum_{n|d} \mu(\frac d n)g(d)
g(n)=n∣d∑f(d)⇒f(n)=n∣d∑μ(nd)g(d)
这个g函数的设置看似很无厘头,但是在某些时候函数g会比函数f好算得多,这可以加速计算。
2 相关证明(卷积法)
引自:知乎
2.1
定义f,g为两个函数。则其卷积运算为:
(
f
∗
g
)
(
n
)
=
∑
i
j
=
n
f
(
i
)
g
(
j
)
(f*g)(n)=\sum_{ij=n}f(i)g(j)
(f∗g)(n)=ij=n∑f(i)g(j)
2.2
我们定义卷积意义下的单位元
ι
\iota
ι,满足:
ι
∗
f
=
f
\iota*f=f
ι∗f=f
即是:
(
ι
∗
f
)
(
n
)
=
∑
i
j
=
n
ι
(
i
)
f
(
i
)
=
f
(
n
)
(\iota*f)(n)=\sum_{ij=n}\iota(i)f(i)=f(n)
(ι∗f)(n)=ij=n∑ι(i)f(i)=f(n)
由卷积的运算,易知:
ι
=
{
1
,
n
=
1
0
,
n
≠
1
\iota= \left\{ \begin{aligned} 1,n=1\\ 0,n≠1\\ \end{aligned} \right.
ι={1,n=10,n̸=1
2.3
我们定义乘法单位元为在普通乘法
(
f
g
)
(
n
)
=
f
(
n
)
g
(
n
)
(fg)(n)=f(n)g(n)
(fg)(n)=f(n)g(n)意义下的单位元,记作
u
u
u。即:
∀
i
∈
I
u
(
i
)
=
1
\forall_{i\in I}u(i)=1
∀i∈Iu(i)=1
2.4
我们将
u
u
u在卷积意义下的逆元称为莫比乌斯函数。即它们满足:
u
∗
μ
=
ι
u*\mu=\iota
u∗μ=ι
对于
μ
\mu
μ满足性质
∑
d
∣
n
μ
(
d
)
=
0
,
n
>
1
\sum_{d|n}\mu(d)=0,n>1
∑d∣nμ(d)=0,n>1
证明如下:
∑
d
∣
n
μ
(
d
)
=
μ
(
1
)
+
μ
(
p
1
)
+
…
+
μ
(
p
k
)
+
μ
(
p
1
p
2
)
+
μ
(
p
1
p
2
…
p
k
)
=
(
k
0
)
+
(
k
1
)
(
−
1
)
+
(
k
2
)
(
−
1
)
2
+
…
+
(
k
k
)
(
−
1
)
k
=
(
1
−
1
)
k
=
0
\begin{aligned} \sum_{d|n}\mu(d)&=\mu(1)+\mu(p_1)+\ldots +\mu(p_k)+\mu(p_1p_2)+\mu(p_1p_2\ldots p_k)\\ &=\left(\begin{aligned}k\\0 \end{aligned}\right)+ \left(\begin{aligned}k\\1 \end{aligned}\right)(-1)+ \left(\begin{aligned}k\\2 \end{aligned}\right)(-1)^2+\ldots+\left(\begin{aligned}k\\k \end{aligned}\right)(-1)^k\\ &=(1-1)^k=0 \end{aligned}
d∣n∑μ(d)=μ(1)+μ(p1)+…+μ(pk)+μ(p1p2)+μ(p1p2…pk)=(k0)+(k1)(−1)+(k2)(−1)2+…+(kk)(−1)k=(1−1)k=0
2.5
那么莫比乌斯反演用数学语言描述就是:
当且仅当
f
(
n
)
=
∑
d
∣
n
g
(
d
)
f(n)=\sum_{d|n}g(d)
f(n)=d∣n∑g(d)
满足
g
(
n
)
=
∑
d
∣
n
μ
(
n
d
)
f
(
d
)
g(n)=\sum_{d|n}\mu(\frac n d)f(d)
g(n)=d∣n∑μ(dn)f(d)
换而言之: f = g ∗ u ⇔ g = f ∗ μ f=g*u\Leftrightarrow g=f*\mu f=g∗u⇔g=f∗μ
将左边的式子两边同乘 μ \mu μ即可证明之
3 性质
n ∑ d ∣ n μ ( d ) d = ϕ ( n ) n\sum_{d|n} \frac {\mu(d)}d=\phi(n) n∑d∣ndμ(d)=ϕ(n)
证明
对于
ϕ
(
x
)
=
x
(
1
−
1
p
1
)
(
1
−
1
p
2
)
…
(
1
−
1
p
k
)
\phi(x)=x(1-\frac 1 {p_1})(1-\frac 1 {p_2})…(1-\frac 1 {p_k})
ϕ(x)=x(1−p11)(1−p21)…(1−pk1),我们把除x之外的所有数都乘起来,得
ϕ
(
x
)
=
x
(
1
−
∑
i
=
1
k
1
p
i
+
∑
i
,
j
=
1
,
i
̸
=
j
k
1
p
i
p
j
…
)
\phi(x)=x(1-\sum_{i=1}^k\frac 1 {p_i}+\sum_{i,j=1,i\not=j}^k\frac 1 {p_ip_j}…)
ϕ(x)=x(1−i=1∑kpi1+i,j=1,i̸=j∑kpipj1…)
而我们可以发现其式子中的
p
i
,
p
i
p
j
,
p
i
p
j
p
r
…
p_i,p_ip_j,p_ip_jp_r…
pi,pipj,pipjpr…就是x的所有约数,而其相应的符号则与上面的莫比乌斯函数之定义相符,也就是说
ϕ
(
x
)
=
∑
d
∣
x
x
d
μ
(
d
)
=
∑
d
∣
x
d
μ
(
x
d
)
\phi(x)=\sum_{d|x} \frac x d \mu(d)=\sum_{d|x} d\mu(\frac x d)
ϕ(x)=∑d∣xdxμ(d)=∑d∣xdμ(dx)
4 线性筛求莫比乌斯函数
void init()
{
miu[1]=1;
for(int i=2;i<=50000;i++)
{
if(!is[i]){pri[++tot]=i;miu[i]=-1;}
for(int j=1;j<=tot;j++)
{
int k=i*pri[j];
if(k>50000) break;
is[k]=1;
if(i%pri[j]==0){miu[k]=0;break;}
else miu[k]=-miu[i];
}
}
}
实际应用时,还可以做底数优化,对其进行前缀和处理,可将复杂度优化至 O ( n ) O(\sqrt n) O(n)。详见例题。
5 例题
[HAOI2011] Problem b(bzoj2301)
题意:给定五个整数a,b,c,d,e,求满足 g c d ( x , y ) = e , x ∈ [ a , b ] , y ∈ [ c , d ] gcd(x,y)=e,x\in[a,b],y\in[c,d] gcd(x,y)=e,x∈[a,b],y∈[c,d]的数对(x,y)对数。
题解:对于 a = c = 1 a=c=1 a=c=1的特殊情况,就相当于求满足 g c d ( x , y ) = 1 , x ∈ [ 1 , b / e ] , y ∈ [ 1 , d / e ] gcd(x,y)=1,x\in[1,b/e],y\in[1,d/e] gcd(x,y)=1,x∈[1,b/e],y∈[1,d/e]的数对对数。
我们不妨设 f ( i ) f(i) f(i)表示条件为 g c d ( x , y ) = i gcd(x,y)=i gcd(x,y)=i时的答案。则f(1)即为所求。构造 g ( i ) g(i) g(i)表示满足 g c d ( x , y ) = k ∗ i , k ∈ Z ∗ gcd(x,y)=k*i,k\in Z^* gcd(x,y)=k∗i,k∈Z∗的数对对数。那么可以知道 g ( i ) = b i ∗ d i g(i)=\frac b i*\frac d i g(i)=ib∗id,就可以愉快地求f(i)的值了。
然而时间复杂度是 O ( n ) O(n) O(n),TLE妥妥的。考虑怎么优化,注意到在某一段区间内的 n i \frac n i in取值相等(因为向下取整了嘛),那么就一起计算。然后脑补一下,你就会发现 n / ( n / i ) n/(n/i) n/(n/i)就是取值相等的段末位置。时间复杂度 O ( n ) O(\sqrt n) O(n)。
最后再用容斥乱搞一通即可。
代码:
#include <iostream>
#include <cstdio>
using namespace std;
typedef long long ll;
const int maxn=50010;
int z,a,b,c,d,e,tot,pri[maxn],is[maxn],miu[maxn],sum[maxn];
ll ans;
inline int min(int x,int y){return x<y?x:y;}
void init()
{
miu[1]=1;
for(int i=2;i<=50000;i++)
{
if(!is[i]){pri[++tot]=i;miu[i]=-1;}
for(int j=1;j<=tot;j++)
{
int k=i*pri[j];
if(k>50000) break;
is[k]=1;
if(i%pri[j]==0){miu[k]=0;break;}
else miu[k]=-miu[i];
}
}
sum[1]=miu[1];
for(int i=2;i<=50000;i++)
sum[i]=sum[i-1]+miu[i];
}
ll check(int x,int y)
{
ll res=0;
int t=min(x,y);
for(int i=1,j;i<=t;i=j+1)
{
j=min(x/(x/i),y/(y/i));
res+=(ll)(sum[j]-sum[i-1])*(x/i)*(y/i);
}
return res;
}
int main()
{
#ifndef ONLINE_JUDGE
freopen("in.txt","r",stdin);
#endif
init();
scanf("%d",&z);
while(z--)
{
scanf("%d%d%d%d%d",&a,&b,&c,&d,&e);
ans=check(b/e,d/e)-check((a-1)/e,d/e)-check(b/e,(c-1)/e)+check((a-1)/e,(c-1)/e);
printf("%lld\n",ans);
}
return 0;
}