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

引入

什么是反演

现在我们有一个函数$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'<=x,所以我们在逆求f(x)之前,所有f(x’)和g(x’)和μ(x’)都是已知的。我们就可以用已知的f(x’)g(x’)μ(x’)和未知的μ(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=nm-\sum{( 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;
}


分享至ヾ(≧∇≦*)ゝ:
分类: 所有

3 条评论

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

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

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

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

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

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

发表评论

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

你是机器人吗? =。= *