莫比乌斯反演

谁践踏了优雅 2022-08-03 00:15 231阅读 0赞

定理：F(n)和f(n)是定义在非负整数集合上的两个函数，并且满足条件F(n)=∑d|nf(d)，那么我们得到结论

f(n)=∑d|nμ(d)F(nd)

在上面的公式中有一个函数，它的定义如下：  
（1）若d=1，那么μ(d)=1  
（2）若d=p1p2..pk，pi均为互异素数，那么μ(d)=(−1)k  
（3）其它情况下μ(d)=0

对于函数，它有如下的常见性质：

（1）对任意正整数有

![这里写图片描述][20140416161656796]

（2）对任意正整数有  
![这里写图片描述][20140416162025500]

线性筛选求莫比乌斯反演函数代码。

bool check[MAXN];
    int prime[MAXN];
    int mu[MAXN];
    void Moblus()
    {
        CLR(check,false);
        mu[1]=1;
        int tot=0;
        for(int i=2; i<=MAXN; i++)
        {
            if(!check[i])
            {
                prime[tot++]=i;
                mu[i]--;
            }
            for(int j=0; j<tot; j ++)
            {
                if(i*prime[j] > MAXN) break;
                check[i*prime[j]]=true;
                if(i%prime[j]==0)
                {
                    mu[i*prime[j]]=0;
                    break;
                }
                else
                {
                    mu[i*prime[j]]=-mu[i];
                }
            }
        }
    }

> 有了上面的知识，现在我们来证明莫比乌斯反演定理。  
> 证明  
> ![这里写图片描述][20140416162820796]  
> 证明完毕！

莫比乌斯反演是一个数学模型，很多问题都可以通过套用这个模型实现简化。  
下面这题是HDU 1695  
分析：对于本题，因为是GCD(x,y)=k，所以GCD(x/k,y/k)=1.  
也就是说，现在问题转化为：  
在区间\[1,b/k\]和\[1,d/k\]中，存在多少个有序对使得(x,y)互质，这个问题就简单啦。做这类问题的最适合的就是莫比乌斯反演。  
因为是有序对，不妨设x≤y，那么我们如果枚举每一个y，小于y有多少个x与y互素，这正是欧拉函数。所以我们可以递推法求欧拉函数，将得到的答案乘以2即可，但是这里乘以2后还有漏计算了的，那么有哪些呢？  
是x==y且为素数的情况，再加上就行了。

#include<stdio.h>
    #include<algorithm>
    #include<vector>
    
    using namespace std;
    
    typedef __int64 ll;
    
    ll prime[100005];
    bool flag[100005];
    ll phi[100005];
    vector<ll>link[100005];
    
    void init()//得到素数以及欧拉函数值
    {
        ll i,j,num=0;
        phi[1]=1;
        for(i=2;i<=100000;i++)
        {
            if(!flag[i])
            {
                prime[num++]=i;
                phi[i]=i-1;
            }
            for(j=0;j<num&&prime[j]*i<=100000;j++)
            {
                flag[prime[j]*i]=true;
                if(i%prime[j]==0)
                {
                    phi[i*prime[j]]=phi[i]*prime[j];
                    break;
                }
                else
                    phi[i*prime[j]]=phi[i]*(prime[j]-1);
            }
        }
        for(j=1;j<=100000;j++)//得到所有数包含的素因子
        {
            ll tmp=j;
            for(i=0;prime[i]*prime[i]<=tmp;i++)
            {
                if(tmp%prime[i]==0)
                {
                    link[j].push_back(prime[i]);
                    tmp=tmp/prime[i];
                    while(tmp%prime[i]==0)
                        tmp=tmp/prime[i];
                }
                if(tmp==1)
                    break;
            }
            if(tmp>1)
                link[j].push_back(tmp);
        }
    }
    
    ll dfs(ll x,ll b,ll now)//容斥原理
    {
        ll i,res=0;
        for(i=x;i<link[now].size();i++)
            res=res+b/link[now][i]-dfs(i+1,b/link[now][i],now);
        return res;
    }
    
    int main()
    {
        init();
        ll i,a,b,t,T,ans,c,d,k;
        while(scanf("%I64d",&T)!=EOF)
        {
            for(t=1;t<=T;t++)
            {
                ans=0;
                scanf("%I64d%I64d%I64d%I64d%I64d",&a,&b,&c,&d,&k);
                if(k==0||k>b||k>d)
                {
                    printf("Case %I64d: 0\n",t);
                    continue;
                }
                if(b>d)
                    swap(b,d);
                b=b/k,d=d/k;
                for(i=1;i<=b;i++)
                    ans=ans+phi[i];
                for(i=b+1;i<=d;i++)
                    ans=ans+b-dfs(0,b,i);
                printf("Case %I64d: %I64d\n",t,ans);
            }
        }
        return 0;
    }

这当然是可以过的，但是如果数据范围再大一些的话，我们试着用莫比乌斯反演来解决。  
这是最朴素的莫比乌斯反演写法，时间复杂度为O(n)，来自[kuangbin的博客][kuangbin]

/* *********************************************** Author :kuangbin Created Time :2013/8/21 19:32:35 File Name :F:\2013ACM练习\专题学习\数学\莫比乌斯反演\HDU1695GCD.cpp ************************************************ */
    
    #include <stdio.h>
    #include <string.h>
    #include <iostream>
    #include <algorithm>
    #include <vector>
    #include <queue>
    #include <set>
    #include <map>
    #include <string>
    #include <math.h>
    #include <stdlib.h>
    #include <time.h>
    using namespace std;
    const int MAXN = 1000000;
    bool check[MAXN+10];
    int prime[MAXN+10];
    int mu[MAXN+10];
    void Moblus()
    {
        memset(check,false,sizeof(check));
        mu[1] = 1;
        int tot = 0;
        for(int i = 2; i <= MAXN; i++)
        {
            if( !check[i] )
            {
                prime[tot++] = i;
                mu[i] = -1;
            }
            for(int j = 0; j < tot; j++)
            {
                if(i * prime[j] > MAXN) break;
                check[i * prime[j]] = true;
                if( i % prime[j] == 0)
                {
                    mu[i * prime[j]] = 0;
                    break;
                }
                else
                {
                    mu[i * prime[j]] = -mu[i];
                }
            }
        }
    }
    int main()
    {
        //freopen("in.txt","r",stdin);
        //freopen("out.txt","w",stdout);
        int T;
        int a,b,c,d,k;
        Moblus();
        scanf("%d",&T);
        int iCase = 0;
        while(T--)
        {
            iCase++;
            scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
            if(k == 0)
            {
                printf("Case %d: 0\n",iCase);
                continue;
            }
            b /= k;
            d /= k;
            if(b > d)swap(b,d);
            long long ans1 = 0;
            for(int i = 1; i <= b;i++)
                ans1 += (long long)mu[i]*(b/i)*(d/i);
            long long ans2 = 0;
            for(int i = 1;i <= b;i++)
                ans2 += (long long)mu[i]*(b/i)*(b/i);
            ans1 -= ans2/2;
            printf("Case %d: %I64d\n",iCase,ans1);
        }
        return 0;
    }

我们知道莫比乌斯反演的一般描述为：

F(n)=∑d|nf(d)⇒f(n)=∑d|nμ(d)F(nd)

其实它还有另一种描述，本题也是用到这种。那就是：

F(n)=∑n|df(d)⇒f(n)=∑n|dμ(dn)F(d)

对于本题，我们设  
f(d)为满足gcd(x,y)=d且1≤x≤n和1≤y≤m的(x,y)的对数  
F(d)为满足d|gcd(x,y)且1≤x≤n和1≤y≤m的(x,y)的对数  
那么，很显然F(x)=nx×mx，反演后得到f(x)=∑x|dμ(dx)F(d)=∑x|dμ(d|x)nd×md  
因为题目要求是GCD(x)=1，那么我们枚举每一个数，然后得到

ans=∑1min(b,d)μ(i)∗(b/i)∗(d/i)−∑1min(b,d)μ(i)∗(b/i)∗(b/i)

这就是上面  O(n) 的做法，但是在数据范围再大一些的时候就会TLE。  
那么我们可以观察到对于  b/i 和  d/i ，在计算机内以int计算是下取整，也就说对于i的变化，  b/i 的变化较小，因此可以通过预处理出莫比乌斯系数的前缀和，然后通过分块计算来加速，从而达到  O(n√) 的复杂度。

// whn6325689
    // Mr.Phoebe
    // http://blog.csdn.net/u013007900
    #include <algorithm>
    #include <iostream>
    #include <iomanip>
    #include <cstring>
    #include <climits>
    #include <complex>
    #include <fstream>
    #include <cassert>
    #include <cstdio>
    #include <bitset>
    #include <vector>
    #include <deque>
    #include <queue>
    #include <stack>
    #include <ctime>
    #include <set>
    #include <map>
    #include <cmath>
    #include <functional>
    #include <numeric>
    #pragma comment(linker, "/STACK:1024000000,1024000000")
    
    using namespace std;
    
    #define eps 1e-9
    #define PI acos(-1.0)
    #define INF 0x3f3f3f3f
    #define LLINF 1LL<<62
    #define speed std::ios::sync_with_stdio(false);
    
    typedef long long ll;
    typedef unsigned long long ull;
    typedef long double ld;
    typedef pair<ll, ll> pll;
    typedef complex<ld> point;
    typedef pair<int, int> pii;
    typedef pair<pii, int> piii;
    typedef vector<int> vi;
    
    #define CLR(x,y) memset(x,y,sizeof(x))
    #define CPY(x,y) memcpy(x,y,sizeof(x))
    #define clr(a,x,size) memset(a,x,sizeof(a[0])*(size))
    #define cpy(a,x,size) memcpy(a,x,sizeof(a[0])*(size))
    
    #define mp(x,y) make_pair(x,y)
    #define pb(x) push_back(x)
    #define lowbit(x) (x&(-x))
    
    #define MID(x,y) (x+((y-x)>>1))
    #define ls (idx<<1)
    #define rs (idx<<1|1)
    #define lson ls,l,mid
    #define rson rs,mid+1,r
    #define root 1,1,n
    
    template<class T>
    inline bool read(T &n)
    {
        T x = 0, tmp = 1;
        char c = getchar();
        while((c < '0' || c > '9') && c != '-' && c != EOF) c = getchar();
        if(c == EOF) return false;
        if(c == '-') c = getchar(), tmp = -1;
        while(c >= '0' && c <= '9') x *= 10, x += (c - '0'),c = getchar();
        n = x*tmp;
        return true;
    }
    template <class T>
    inline void write(T n)
    {
        if(n < 0)
        {
            putchar('-');
            n = -n;
        }
        int len = 0,data[20];
        while(n)
        {
            data[len++] = n%10;
            n /= 10;
        }
        if(!len) data[len++] = 0;
        while(len--) putchar(data[len]+48);
    }
    //-----------------------------------
    
    
    const int MAXN=100010;
    bool check[MAXN];
    int prime[MAXN];
    int mu[MAXN];
    void Moblus()
    {
        CLR(check,false);
        mu[1]=1;
        int tot=0;
        for(int i=2; i<=MAXN; i++)
        {
            if(!check[i])
            {
                prime[tot++]=i;
                mu[i]--;
            }
            for(int j=0; j<tot; j ++)
            {
                if(i*prime[j] > MAXN) break;
                check[i*prime[j]]=true;
                if(i%prime[j]==0)
                {
                    mu[i*prime[j]]=0;
                    break;
                }
                else
                {
                    mu[i*prime[j]]=-mu[i];
                }
            }
        }
    }
    int sum[MAXN];
    
    long long solve(int n,int m)
    {
        long long ans=0;
        if(n>m)swap(n,m);
        for(int i=1, la=0; i<=n; i=la+1)
        {
            la=min(n/(n/i),m/(m/i));
            ans+=(long long)(sum[la]-sum[i-1])*(n/i)*(m/i);
        }
        return ans;
    }
    int main()
    {
    
        Moblus();
        sum[0]=0;
        for(int i=1;i<=MAXN;i++)
            sum[i]=sum[i-1]+mu[i];
        int a,b,c,d,k;
        int T;
        read(T);
        while(T--)
        {
            read(a),read(b),read(c),read(d),read(k);
            ll ans=solve(b/k,d/k) - solve((a-1)/k,d/k) - solve(b/k,(c-1)/k) + solve((a-1)/k,(c-1)/k);
            printf("%lld\n",ans);
        }
        return 0;
    }

一定要用分块优化的题目[BZOJ 2301 : Problem b][BZOJ 2301 _ Problem b]

[20140416161656796]: /images/20220731/4c7566f062ee4c4f801428a2c262350f.png
[20140416162025500]: /images/20220731/603ad85487bc4fcb933d9ea16ab3478d.png
[20140416162820796]: /images/20220731/042c686e4f75481097afe99980d0fdca.png
[kuangbin]: http://www.cnblogs.com/kuangbin/archive/2013/08/21/3273440.html
[BZOJ 2301 _ Problem b]: http://www.lydsy.com/JudgeOnline/problem.php?id=2301