BZOJ 2301 Problem b（莫比乌斯反演+容斥原理）

2301: [HAOI2011]Problem b

Time Limit: 50 Sec  Memory Limit: 256 MB
Submit: 7205  Solved: 3425
[Submit][Status][Discuss]

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

 

 

要求两个区间范围内的数字的gcd为k的对个数。这种两个数字范围的题目按照套路都是要用容斥原理拆分成四个子问题。每个子问题形式为cal(a,b,k)，即两个数字的取值范围分别是[1,a]和[1,b]，然后去求解个数。

我们考虑把询问写成gcd(a,b,k)的形式，那么显然有gcd(a,b,k)==gcd(a/l,b/k,1)。我们令f(i)=gcd(a,b,i)，即f(i)为gcd(a,b)==i这样的(a,b)个数，然后令F(i)为i | gcd(a,b)这样的(a,b)的对数。显然有 。根据定义，我们可以写出以下式子：

                                                               

即我们要求的答案求和的结果是F(i)。我们发现F(i)的通项公式已经给出，但是这个f(i)却不是那么的好求。于是，我们利用莫比乌斯反演定理，利用F(i)反演求出f(i)。莫比乌斯反演了解一下：莫比乌斯反演的两种形式及其证明

                 

推出f(i)表达式之后，我们就可以对gcd(a,b,k)求解了。我们首先把问题变成gcd(a/k,b/k,1)，这样对应的就是f(1)。

                                                   

如此直接枚举d求解的话，一次询问复杂度为O(N)，总共q次询问就是O(Nq)，然而这个复杂度显然是不行的。接着，我们注意到对于一个数字n，其对应的 的取值最多只有，于是我们可以利用这个性质进行分段优化。举个例子说，对于n=100，d=34，此时 的值等于2，而当d循环到50的时候，取值还是2，这意味着这一段我就可以一起统计。预处理莫比乌斯函数的前缀和，然后对于同一段的用这一段的莫比乌斯函数之和乘以F(d)。这样，总的复杂度就可以优化到。具体见代码：

#include<bits/stdc++.h>
#define IO ios::sync_with_stdio(0);cin.tie(0);cout.tie(0)
#define LL long long
#define N 50010

using namespace std;

int s[N],mu[N],p[N];
bool isp[N];

void init()
{
    int sz=0;
    s[1]=mu[1]=1; 
    memset(isp,1,sizeof(isp));
    for(int i=2;i<N;i++)
    {
        if(isp[i]) p[++sz]=i,mu[i]=-1;
        for(int j=1;j<=sz&&(LL)i*p[j]<N;j++)
        {
            isp[i*p[j]]=0;
            if(i%p[j]==0)
            {
                mu[p[j]*i]=0; break;
            } else mu[p[j]*i]=-mu[i];
        }
        s[i]=s[i-1]+mu[i];
    }
}

int cal(int n,int m)
{
    int res=0,last;
    for(int i=1;i<=n&&i<=m;i=last+1)
    {
        last=min(n/(n/i),m/(m/i));
        res+=(s[last]-s[i-1])*(n/i)*(m/i);
    }
    return res;
}

int main()
{
    int a,b,c,d,k,n;
    IO; init(); cin>>n;
    while(n--)
    {
        cin>>a>>b>>c>>d>>k;
        a--; c--; b/=k; d/=k; a/=k; c/=k;
        cout<<cal(b,d)-cal(a,d)-cal(b,c)+cal(a,c)<<endl;
    }
    return 0;
}