Problem b（BZOJ-2301/HAOI-2011）

Problem 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

思路：

与 GCD（HDU-1695） 相似，但 x、y 不是从 1 开始，因此需要使用容斥原理来进行处理

题目实质是要求：

利用容斥定理来进行转换，设

那么有：

因此，现在只需要考虑计算  即可

设 f(d) 为满足 gcd(i,j)=d 的个数，即：

设 g(n) 为满足 gcd(i,j)=d 的倍数的个数，即：

可以看出，f(d) 与 g(n) 符合莫比乌斯反演的形式：

那么，对于  进行化简

将 f(k) 代入，有：

根据莫比乌斯反演，有：

设枚举项  为 t

那么有 

为了优化时间，可以将  的 k 提出去，即有：

此时的复杂度为 O(n)，由于是多组查询，发现式子中有整除，利用整除分块求 mu 的前缀和，进行优化即可

Source Program

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<string>
#include<cstring>
#include<cmath>
#include<ctime>
#include<algorithm>
#include<utility>
#include<stack>
#include<queue>
#include<vector>
#include<set>
#include<map>
#include<bitset>
#define EPS 1e-9
#define PI acos(-1.0)
#define INF 0x3f3f3f3f
#define LL long long
#define Pair pair<int,int>
const int MOD = 1E9+7;
const int N = 100000+5;
const int dx[] = {1,-1,0,0,-1,-1,1,1};
const int dy[] = {0,0,-1,1,-1,1,-1,1};
using namespace std;
int mu[N];
int prime[N];
bool bprime[N];
int cnt;
LL sum[N];
void getMu(int n){//线性筛求莫比乌斯函数
    cnt=0;
    mu[1]=1;//根据定义，μ(1)=1
    memset(bprime,false,sizeof(bprime));

    for(int i=2;i<=n;i++){//求2~n的莫比乌斯函数
        if(!bprime[i]){
            prime[++cnt]=i;//存储质数
            mu[i]=-1;//i为质数时，μ(1)=-1
        }
        for(int j=1;j<=cnt&&i*prime[j]<=n;j++){//枚举i之前的素数个数
            bprime[i*prime[j]]=true;//不是质数
            if(i%prime[j])//i不是prime[j]的整数倍时，i*prime[j]就不会包含相同质因子
                mu[i*prime[j]]=-mu[i];//mu[k]=mu[i]*mu[prime[j]]，因为prime[j]是质数，mu值为-1
            else{
                mu[i*prime[j]]=0;
                break;//留到后面再筛
            }
        }
    }

    for(int i=1;i<=n;i++)
        sum[i]=sum[i-1]+mu[i];
}
LL cal(LL n,LL m,LL k){//计算a在[1,n],b在[1,m]中GCD(a,b)=1的个数
    n/=k;
    m/=k;
    LL minn=min(n,m);
    LL res=0;
    for(int left=1,right;left<=minn;left=right+1){
        right=min(n/(n/left),m/(m/left));
        res+=(1LL)*(n/left)*(m/left)*(sum[right]-sum[left-1]);
    }
    return res;
}

int main(){
    getMu(100000);

    int t;
    scanf("%d",&t);
    int Case=1;
    while(t--){
        int a,b,c,d,k;
        scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);

        LL res=cal(b,d,k)-cal(a-1,d,k)-cal(c-1,b,k)+cal(a-1,c-1,k);
        printf("%lld\n",res);
    }
    return 0;
}