莫比乌斯反演初步

莫比乌斯反演是数论中很重要的内容,可以优雅解决很多数论的问题。之前浅学了一下,这两天寒假专题训练又拾起来刷了些题,想总结一下。由于内容较多本文证明部分会跳过一些,更着重于莫比乌斯反演在解题中的实际应用。

预备知识:

首先介绍一下什么是莫比乌斯函数

莫比乌斯函数μ(n)的定义域是N
2)μ(1)=1
3)当n存在平方因子时,μ(n)=0
4)当n是素数或奇数个不同素数之积时,μ(n)=-1
5)当n是偶数个不同素数之积时,μ(n)=1
如果大家熟悉欧拉函数的话会觉得两者非常相似,我们可以在素数的线性筛法中用O(n)的时间得到全部1-n的莫比乌斯函数。具体做法可以参考jzp大神的ppt 点击打开链接,非常推荐这个ppt,证明很完善思路也很清晰。

然后就是 莫比乌斯反演。 莫比乌斯反演有两种不同的形式,f和F为不同的函数,mu(x)表示x的莫比乌斯函数。

 

以上这两个式子就是莫比乌斯反演的基本形式了。下面从一道题目出发看莫比乌斯反演在解题的时候具体怎么用。

B - GCD

Time Limit:3000MS    Memory Limit:32768KB    64bit IO Format:%I64d & %I64u
SubmitStatus

Description

Given 5 integers: a, b, c, d, k, you're to find x in a...b, y in c...d that GCD(x, y) = k. GCD(x, y) means the greatest common divisor of x and y. Since the number of choices may be very large, you're only required to output the total number of different number pairs.
Please notice that, (x=5, y=7) and (x=7, y=5) are considered to be the same.

Yoiu can assume that a = c = 1 in all test cases.
 

Input

The input consists of several test cases. The first line of the input is the number of the cases. There are no more than 3,000 cases.
Each case contains five integers: a, b, c, d, k, 0 < a <= b <= 100,000, 0 < c <= d <= 100,000, 0 <= k <= 100,000, as described above.
 

Output

For each test case, print the number of choices. Use the format in the example.
 

Sample Input

       
       
2 1 3 1 5 1 1 11014 1 14409 9
 

Sample Output

       
       
Case 1: 9 Case 2: 736427

Hint

For the first sample input, all the 9 pairs of numbers are (1, 1), (1, 2), (1, 3), (1, 4), (1, 5), (2, 3), (2, 5), (3, 4), (3, 5).
         
                

在问题中我们需要求的是1<=x<=n,1<=y<=m中gcd(x,y)==k的组数,也就是

其中[x]表示sigma(x)。显然这里是可以用容斥搞的,但之后我们会发现莫比乌斯反演法的话更加好写,也快得多

那么这个地方怎么用莫比乌斯反演法呢

 

我们让f(n,m,k)表示1<=x<=n,1<=y<=m中gcd(x,y)==k的组数,也就是(*)式

再让F(n,m,k)表示1<=x<=n,1<=y<=m中gcd(x,y)是k的倍数的组数

 

F(n,m,k)很好求,显然F(n,m,k) =

而题目的目标是求f(n,m,k)

 

那么很显然,满足莫比乌斯反演的第二种形式

根据莫比乌斯反演法,我们得到

由于F(d)我们是知道的,所以有了这个式子f(k)就可以求出来了,只要枚举每个k的倍数作为d然后求和就可以了

到这里这道题已经可做了。

 

但这样还不够好:复杂度无法降低,莫比乌斯反演法很多题目对f(k)运算速度要求sqrt(n)甚至是更低,而现在这个的复杂度是min(n,m)/k。

 

不过幸好我们稍微转换一下思考方式就能得到一个形式更简单的算法:

 

考虑下面这句话:1<=x<=n 1<=y<=m中 gcd(x,y)==k 的对数相当于 1<=x<=n/k 1<=y<=m/k 中gcd(x,y) == 1的对数。

这句话是很显然的,因为gcd(x,y)==k每一对的x和y除k之后必然映射到一组1<=x<=n/k1<=y<=m/k中一组gcd(x,y) == 1的x和y,反之亦然。

这样的话,我们就可以把从gcd(x,y)==k的问题全部转换成在求gcd(x,y)==1的问题,只是上下界变动一下而已。

 

我们带入k=1,

我们知道f(1)表示的是gcd(x,y)==1,也就是x,y互质的对数。

结合上面的分析不难发现

那么我们就把f(k)的求解转换成了一端连续值的求和了,从编程角度比之前简单了许多。

到这一步对于解决上面一道题来说是游刃有余了。

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<string>
#include<algorithm>
#include<cmath>
using namespace std;

typedef long long ll;
#define FILL(x) memset((x),0,sizeof(x))
const int maxn = 100000+20;
bool check[maxn];
int mu[maxn];
int prime[maxn];
void Mobius(int N)
{
    FILL(check);
    mu[1] = 1;
    int tot = 0;
    for(int i=2;i<=N;i++){
        if(!check[i]){
            prime[tot ++] = i;
            mu[i] = -1;
        }
        for(int j=0;j<tot;j++){
            if(i * prime[j] > N)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 a,b,c,d,k;
int cases;
ll ans;
void solve()
{
    ans = 0;
    if(k == 0){
        printf("Case %d: 0\n",cases);
        return ;
    }
    b/=k;d/=k;
    if(b>d)swap(b,d);
    for(int i=1;i<=b;i++){
        ans += (ll)mu[i] * (b/i) * (d/i);
    }
    ll ans2 = 0;
    for(int i=1;i<=b;i++){
        ans2 += (ll)mu[i] * (b/i) * (b/i);
    }
    ans -= ans2/2;
    printf("Case %d: %I64d\n",cases,ans);
}
int main()
{
    Mobius(100001);
    int T;
    scanf("%d",&T);
    cases = 0;
    while(T--){
        cases ++;
        scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
        solve();
    }
    return 0;
}


 但是正如之前所说有些题目对f(k)的求解速度要求比较严格,比如下面这道题:


G -Problem b

Time Limit:50000MS    Memory Limit:262144KB    64bit IO Format:%lld & %llu

Description

对于给出的n个询问,每次求有多少个数对(x,y),满足a≤x≤b,c≤y≤d,且gcd(x,y) = k,gcd(x,y)函数为x和y的最大公约数。

Input

第一行一个整数n,接下来n行每行五个整数,分别表示a、b、c、d、k 

Output

共n行,每行一个整数表示满足要求的数对(x,y)的个数


Sample Input

2
2 5 1 5 1
1 5 1 5 2

Sample Output

14
3

Hint

100%的数据满足:1≤n≤50000,1≤a≤b≤50000,1≤c≤d≤50000,1≤k≤50000



     这道题数据量很大,每组数据O(n)处理的话是不够的,我们需要分块加速,这也是转换成求f(1)的最大好处之一。

转换成gcd(x,y)==1之后我们发现d的取值是连续的,n/dm/dd在一定范围的时候是不会变化的,我们称这样的区间叫不变块。我们可以把每一个不变块中的mu[i]求和之后乘上n/dm/d加到答案里。

           每个不变块中的求和工作只要预处理一下mu[i]的前缀和就能轻松搞定,这样的话我们求和时就可以每次求一块的和,而不是一个个求和了,速度大大加快,达到sqrt(n)。

            而d属于[x,n/(n/x)]这个区间中时n/d的值是不变的,对m也是一样。那么假设当前块的初始点为x,这个块的末尾点就是min(n/(n/x),m/(m/x))(因为n/d和m/d都要保持不变)。然后就可以一块一块求和了。

     说上去也许稍微有点复杂,但看一下代码应该就很容易理解了

           

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<string>
#include<algorithm>
#include<cmath>
using namespace std;

typedef long long ll;
#define FILL(x) memset((x),0,sizeof(x))
const int maxn = 100000+20;
bool check[maxn];
int mu[maxn];
int prime[maxn];
int sum[maxn];
void Mobius(int N)
{
    FILL(check);
    mu[1] = 1;
    int tot = 0;
    for(int i=2;i<=N;i++){
        if(!check[i]){
            prime[tot ++] = i;
            mu[i] = -1;
        }
        for(int j=0;j<tot;j++){
            if(i * prime[j] > N)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 a,b,c,d,k;
int cases;
ll calc(int b,int d)
{
    ll ans = 0;
    b/=k;d/=k;
    if(b>d)swap(b,d);
    for(int i=1,last=i;i<=b;i=last+1){//i是不变块的左边界
        last = min(b/(b/i),d/(d/i));//last是不变块的右边界
        ans += (ll)(sum[last]-sum[i-1]) * (b/i) * (d/i);
    }
    return ans;
}
void solve()
{
    if(k == 0){
        printf("0\n");
        return ;
    }
    ll ans1 = calc(b,d);
    ll ans2 = calc(a-1,d);
    ll ans3 = calc(b,c-1);
    ll ans4 = calc(a-1,c-1);
    ll ans = ans1 - ans2 - ans3 + ans4;
    printf("%lld\n",ans);
}
int main()
{
    Mobius(100001);
    sum[0] = 0;
    for(int i=1;i<maxn;i++)sum[i]=sum[i-1]+mu[i];
    int T;
    scanf("%d",&T);
    cases = 0;
    while(T--){
        cases ++;
        scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
        solve();
    }
    return 0;
}

到这里应该能基本解决求GCD(x,y)==k的组数这类问题,也是莫比乌斯反演法最基本的应用了。更高级的用法实际上都是从这个基础出发的。更多用到的是莫比乌斯反演最初始的定义式构造函数进行推导了。

    再次强烈推荐jzp大神的ppt点击打开链接,做莫比乌斯反演首先要会线性筛求出莫比乌斯函数。从他的ppt中可以看到另外一种应用方法,紧紧围绕下面两个重要公式(上下分别为莫比乌斯函数和欧拉函数):


   有了这两个式子在做题时很多时候可以省去构造函数的过程,直接套用,在ppt中的应用也非常精彩,推荐推荐。

   内容比较初步,后面可能带来剩下几题题解,大家请随意吐槽指教。


评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值