网上有关莫比乌斯反演的博客很多,但是其中只有少数一部分适合初学者学习。通过那很少的一部分和数学一本通上” 详细” 的证明,和我一个下午的颓废认真研读,终于,我还是没懂什么是莫比乌斯反演,以及为什么鲁迅会那么痴迷于枣树。

引入

什么是反演

现在我们有一个函数 $g(x)=\sum_{i=1}^x{i}$。我们做一件很蠢的事 (但是就是反演的本质),如果有一个函数 $f(x)=x$,那么 $g(x)=\sum_{i=1}^x{f(i)}$。无聊的我们暂时忘记了 f(x) 怎么求 (就是看着 f(x)=x,知道了 x 却不知道 f(x)),现在好心的鲁迅告诉我们其中一棵是枣树 g(x) 是多少,那么我们突然领悟到一个事实: 如果我们知道了 g(x),那么我们可以很快也知道 f(x)。因为对于 g(x) 我们可以这样求:$g(x)=x(1+f(x))/2$,所以有了这个我们就可以反推回 f(x),即 $f(x)=(2g(x)-x)/x$。而这就是很垃圾的反演。

刚才例子的反演变换是 $sum=x(1+x)/2$→$x=\frac{1}{2}(-1+\sqrt{1+8sum})$,虽然很无聊,但是它告诉我们:当我们知道了两个函数之间的正向关系,我们就可以利用这个关系反向求出想要的结果。

如何更好地理解 (准备学习的) 莫比乌斯反演

其实在书上总是公式连天,题解中用法无穷无尽的莫比乌斯反演,也不是很难理解的。
今天我们换个角度去理解它。

先看这么一个问题: 求有多少 (i,j) 满足 gcd(i,j)=1,i<=n,j<=m (n,m<=1000000)
首先我们有 m×n 组 (i,j),然后我们要把 gcd(i,j)!=1 的删去;
如果 i 和 j 都是一个数 x 的倍数,那么这样的 (i,j) 就可以删去,如果枚举每一个 x 就可以删去所有不该保留的 (i,j);
问题来了:有的 i,j 会被删去多次,比如 (6,12) 在 2,3,6 出都被删掉了。
我们就会想,这可能跟容斥有一点关系。
所以我们仔细观察被删除的 (i,j) 的特征。
1 如果已经用 x,y 两个数删过 (i,j),那么 xy 的倍数就被多删了一次,要补回来。
2. 如果已经用 x,y,z 三个数删过 (i,j),且已经用 xy,xz,yz 补回来了,那么 xyz 的倍数就被多补了一次,因此要再删掉一次。
3. 以此类推,我们发现了类似于容斥的规律。而这就是莫比乌斯反演的精髓。

正题

我们先上公式:

如果有
$$g(n)=\sum\limits_{d|n}{f(d)}$$
那么有
$$f(n)=\sum\limits_{d|n}{(u(\frac{n}{d})g(d))}$$
这里用到了一个叫μ的函数。这个函数的名字叫 “莫比乌斯函数”,读作”miu(第四声)”(μ)。这个函数就是莫比乌斯反演的精髓,它掌控着每一个 f(x) 的符号以及容斥与否。

官方对μ(x) 的定义是这样的:(摘自 http://www.cnblogs.com/chenyang920/p/4811995.html)

(1).μ是一个关于整数的函数
(2).μ[x] = 1 当且仅当 x 能够分解成偶数个不同质数的乘积(其中 1 不能被分解,所以 1 的分解出的质数个数是 0,所以μ[1]=1)
(3).μ[x] = -1 当且仅当 x 能够分解成奇数个不同质数的乘积
(4).μ[x] = 0 除开 (2),(3) 的其他情况

然后我们要做的就是用这个函数去 A 题,当然这里还是证明一下。
首先,我们要明白一点:由于 g(x) 可以表示成很多 f(x’) 的和,而 x’ (1).μ(x)=1,x=1: 因为 g(1)=f(1), 所以μ(1)=1;

(2).μ(x)=-1,x=p(质数): 因为 g(p)=f(1)+f(p),f(p)=μ(1)g(p)+μ(p)g(1)。已知μ(1)==1,联立两式得:μ(p)=-1;

(3).μ(x)=1,x=p1*p2: 因为 g(p1p2)=f(1)+f(p1)+f(p2)+f(p1p2),f(p1p2)=μ(1)g(p1p2)+μ(p1)g(p2)+μ(p2)g(p1)+μ(p1p2)g(1)。已知μ(1)=1,μ(p1)=μ(p2)=-1,则μ(p1p2)=+1;

(4).μ(x)=(-1)n,x=p1p2…*pn: 同前两个证明,我们可以退出当 x 为 n 个不同的质数的乘积时,μ(x)=(-1)n;

(5).μ(x)=0,x=prm(r>=2): 算了我不证了,总之此时μ(x)=0; 举个例子:g(12)=f(1)+f(2)+f(3)+f(4)+f(6)+f(12) 得出 f(12)=μ(1)g(12)+μ(2)g(6)+μ(3)g(4)+μ(4)g(3)+μ(6)g(2)+μ(12)g(1)=g(12)-g(6)-g(4)+g(3)+g(2)

这些结论显然是正确的。

所以我们证明了莫比乌斯函数的正确性,它的确是连接冥界和人间的通道,它的确可以从 g(x) 反推出 f(x)。

这时候我们就可以用莫比乌斯反演来解决前文提到的题目: 求有多少 (i,j) 满足 gcd(i,j)=1,i<=n,j<=m (n,m<=1000000)
对于所有 (i,j) 我们用很多数去筛
利用莫比乌斯函数,就有
$$
ans=\sum_{x=1}^{min(n,m)}{( u(x) [\frac{n}{x}][\frac{m}{x}] )}
$$
为什么呢?具体原因在例题中解释。

那么我们又可以发现莫比乌斯函数的一个重要性质:它是一个积性函数!μ(xy)=μ(x)μ(y) (x,y 互质) 利用这一点我们就可以 O(n) 地像线性筛一样筛出 n 以内的所有μ(x) 了!

char vis[MX];
ll prm[MX],miu[MX],pnum;

void init()
{
    miu[1]=1;
    pnum=0;
    for(ll i=2;i<=MX;i++)
    {
        if(!vis[i])prm[++pnum]=i,miu[i]=-1;
        for(ll j=1;j<=pnum;j++)
        {
            if(i*prm[j]>=MX)break;
            vis[i*prm[j]]=1;
            if(i%prm[j]==0)
            {
                miu[i*prm[j]]=0;
                break;
            }
            else miu[i*prm[j]]=-miu[i];
        }
    }
}

应用

我们先来看一道题目:

HDU1695 GCD

给定 a,b,求 GCD(x,y)==k(x<=a,y<=b) 的 (x,y) 个数 (x,y) 与 (y,x) 只算一个。

我们可以知道,这道题其实再求 GCD(x,y)==1(x<=a/k,y<=b/k) 的个数。这个怎么处理呢?

我们设 f(d) 为 gcd(x,y)==d 的数对个数,g(d) 为 gcd(x,y)%d==0 的数对个数。

由于它们构成了同余的关系,所以我们知道这样可以利用反演将 f(x) 转为求 g(x)。为什么要这样转化呢?因为 g(x) 实在是太好求了。由于 gcd(x,y)%d==0,我们就可以像问题一开始消除 k 的方法一样消除 d,$\sum{[gcd(x,y)\% d==0]}(x<=a,y<=b)$等价于 $\frac{a}{d}\times \frac{b}{d}$

那么显而易见:

$$g(n)=\sum\limits_{n|d}{f(d)}(d<=min(a,b))$$
$$f(n)=\sum\limits_{n|d}{u(\frac{d}{n})g(d)}(d<=min(a,b))$$

这也是莫比乌斯反演的另一种形式。这和之前的 d|n 不同,成了 n|d,但是处理方法依旧。进而通过我们刚才发现的 g(x) 的简便求法,我们可以在 O(n) 时间内求出答案。

#include <iostream>
#include <cstdio>
#include <cmath>

#define MX 1000010
#define min(x,y) (x>y?y:x)

using namespace std;

typedef long long ll;

char vis[MX];
ll prm[MX],miu[MX],pnum;

void init()
{
    miu[1]=1,pnum=0;
    for(ll i=2;i<=MX;i++)
    {
        if(!vis[i])prm[++pnum]=i,miu[i]=-1;
        for(ll j=1;j<=pnum;j++)
        {
            if(i*prm[j]>=MX)break;
            vis[i*prm[j]]=1;
            if(i%prm[j]==0)
            {
                miu[i*prm[j]]=0;
                break;
            }
            else miu[i*prm[j]]=-miu[i];
        }
    }
}

int main()
{
    ll t,a,b,k,mn;
    init();
    scanf("%lld",&t);
    for(ll w=1;w<=t;w++)
    {
        scanf("%*d%lld%*d%lld%lld",&a,&b,&k);
        if(k==0)
        {
            cout<<"Case "<<w<<": "<<0<<endl;
            continue;
        }
        a/=k,b/=k;
        mn=min(a,b);
        ll ans=0,t2=0;
        for(ll i=1;i<=mn;i++)ans+=(ll)miu[i]*((a/i)*(b/i)),t2+=(ll)miu[i]*(mn/i)*(mn/i);
        cout<<"Case "<<w<<": "<<ans-t2/2<<endl;
    }
    return 0;
}

分类: 文章

4 条评论

konnyakuxzy · 2017年6月26日 8:04 下午

机房前面有两个蒟蒻,一个是 XZY,另一个也是 XZY

litble · 2017年6月26日 7:11 下午

其实我还是搞不懂您为什么这么痴迷与鲁迅和枣树的关系?

konnyakuxzy · 2017年6月24日 7:37 上午

哇塞,太强啦,感谢感谢 OrzOrzOrz

【算法】 如何感性理解莫比乌斯函数 – K-XZY · 2018年12月25日 8:29 下午

[…] 如果有要学习莫比乌斯反演的同学请看这里:【算法】莫比乌斯反演 -boshi […]

发表评论

电子邮件地址不会被公开。 必填项已用*标注