莫比乌斯反演

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

f(n)=∑d|nμ(d)F(nd)f(n)=\sum_{d|n}\mu(d)F({n\over d})

在上面的公式中有一个函数，它的定义如下：

（1）若d=1d=1，那么μ(d)=1\mu(d)=1

（2）若d=p1p2..pkd=p_1p_2..p_k，pip_i均为互异素数，那么μ(d)=(−1)k\mu(d)=(-1)^k

（3）其它情况下μ(d)=0\mu(d)=0

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

（1）对任意正整数有




（2）对任意正整数有



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

[code]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];
            }
        }
    }
}

有了上面的知识，现在我们来证明莫比乌斯反演定理。

证明





证明完毕！

莫比乌斯反演是一个数学模型，很多问题都可以通过套用这个模型实现简化。

下面这题是HDU 1695

分析：对于本题，因为是GCD(x,y)=kGCD(x,y)=k，所以GCD(x/k,y/k)=1GCD(x/k,y/k)=1.

也就是说，现在问题转化为：

在区间[1,b/k][1,b/k]和[1,d/k][1,d/k]中，存在多少个有序对使得(x,y)(x,y)互质，这个问题就简单啦。做这类问题的最适合的就是莫比乌斯反演。

因为是有序对，不妨设x≤yx\le y，那么我们如果枚举每一个yy，小于yy有多少个xx与yy互素，这正是欧拉函数。所以我们可以递推法求欧拉函数，将得到的答案乘以2即可，但是这里乘以2后还有漏计算了的，那么有哪些呢？

是x==yx==y且为素数的情况，再加上就行了。

[code]#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)O(n)，来自kuangbin的博客

[code]/* ***********************************************
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)=\sum_{d|n} f(d) \Rightarrow f(n)=\sum_{d|n}\mu(d)F({n\over d})

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

F(n)=∑n|df(d)⇒f(n)=∑n|dμ(dn)F(d)F(n)=\sum_{n|d} f(d) \Rightarrow f(n)=\sum_{n|d}\mu({d\over n})F({d})

对于本题，我们设

f(d)f(d)为满足gcd(x,y)=dgcd(x,y)=d且1≤x≤n1\le x \le n和1≤y≤m1\le y \le m的(x,y)(x,y)的对数

F(d)F(d)为满足d|gcd(x,y)d|gcd(x,y)且1≤x≤n1\le x \le n和1≤y≤m1\le y\le m的(x,y)(x,y)的对数

那么，很显然F(x)=nx×mxF(x)={n\over x}\times {m\over x}，反演后得到f(x)=∑x|dμ(dx)F(d)=∑x|dμ(d|x)nd×mdf(x)=\sum_{x|d}\mu({d\over x})F(d)=\sum_{x|d}\mu(d|x){n\over d}\times {m\over d}

因为题目要求是GCD(x)=1GCD(x)=1，那么我们枚举每一个数，然后得到

ans=∑1min(b,d)μ(i)∗(b/i)∗(d/i)−∑1min(b,d)μ(i)∗(b/i)∗(b/i)ans=\sum_1^{\min(b,d)}\mu({i})*(b/i)*(d/i)-\sum_1^{\min(b,d)}\mu({i})*(b/i)*(b/i)

这就是上面O(n)O(n)的做法，但是在数据范围再大一些的时候就会TLE。

那么我们可以观察到对于b/ib/i和d/id/i，在计算机内以int计算是下取整，也就说对于i的变化，b/ib/i的变化较小，因此可以通过预处理出莫比乌斯系数的前缀和，然后通过分块计算来加速，从而达到O(n√)O(\sqrt{n})的复杂度。

[code]//      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