HAOI2011 Problem b 洛谷P2522

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)的个数。

---

Hint

100%的数据满足：1≤n≤50000，1≤a≤b≤50000，1≤c≤d≤50000，1≤k≤50000。

---

Solution

我觉得太神奇了这道题，省选题都那么恶心的吗？？？

首先需要明确的是这道题不用像BZOJ2480/HDUI695那种恶心题一样要去重。(去重语句:ret+=(1-x)*x+x*(y-x))，去重的原理是，把组合按顺序、坐标写出来，对角线以上的不要就行了。

我先就用的普通的莫比乌斯反演，f(k)表示gcd==k的(x,y)的数量，F(n)表示gcd==n的倍数，那么由反演公式F(n)=sigma n|k(f(k)) -->f(n)=sigma n|k(Mo[k/n]*F(k))，其中F(k)即k的倍数两两组合的数量，那么通过反演就计算出了f(k)的值。预处理出Mobius函数，然后就按照这个思路，需要一点容斥原理的知识，因为x的区间是[a,b]，y的区间是[c,d]，并不是BZOJ2480/HDUI695那个题那个样子的[1,b],[1,d](反正那个题恶心的地方在去重好吧)，那么这个题的答案应该就是ans=workk(b,d)+workk(a-1,c-1)-workk(b,c-1)-workk(a-1,d)，需要注意的是a和c都要-1，因为是闭区间。

接下来就是这道题很魔鬼的地方了。

用普通方法做只能A30%的数据，剩下的数据全都要T，因为数据组数的规模是50000，这个询问起来肯定要超时。于是yb老师和左哥大佬都给我说要分 块。

靠，当年的蒲公英简直是心理阴影好吗？？！！

而且分块的思路也很魔鬼，处理出Mobius函数的前缀和，然后，

然后，

然后，

块的大小为TMP=min(x/(x/i),y/(y/i))，然后只需要ret+=(x/i)*(y/i)*(sum[TMP]-sum[i-1])。

你就会发现它A了(落泪了)。

以下是AC代码:

#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define maxn 50005
using namespace std;
int Mo[maxn],primes[maxn],num[maxn],sum[maxn];
int cnt,T,a,b,c,d,k,ans;
inline void Mobius(){
	Mo[1]=1;
	for(int i=2;i<=maxn-1;i++){
		if(!num[i]){
			primes[++cnt]=i;
			Mo[i]=-1;
		}
		for(int j=1;j<=cnt&&i*primes[j]<maxn;j++){
			num[i*primes[j]]=true;
			Mo[i*primes[j]]=Mo[i]*(-1);
			if(i%primes[j]==0){
				Mo[i*primes[j]]=0;
				break;
			}
		}
	}
	for(int i=1;i<=maxn-1;i++){
		sum[i]=sum[i-1]+Mo[i];
	}
}
inline int workk(int x,int y){
	x/=k,y/=k;
	if(x>y)swap(x,y);
	int ret=0,TMP=0;
	for(int i=1;i<=x;i=TMP+1){
		TMP=min(x/(x/i),y/(y/i));
		ret+=(x/i)*(y/i)*(sum[TMP]-sum[i-1]);
	}
	return ret;
}
int main(){
	Mobius();
	scanf("%d",&T);
	for(int i=1;i<=T;i++){
		scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
		ans=workk(b,d)+workk(a-1,c-1)-workk(b,c-1)-workk(a-1,d);
		printf("%d\n",ans);
	}
	return 0;
}

以下是普通做法(30%数据):

#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define maxn 50005
using namespace std;
int Mo[maxn],primes[maxn],num[maxn],sum[maxn];
int cnt,T,a,b,c,d,k,ans;
inline void Mobius(){
	Mo[1]=1;
	for(int i=2;i<=maxn-1;i++){
		if(!num[i]){
			primes[++cnt]=i;
			Mo[i]=-1;
		}
		for(int j=1;j<=cnt&&i*primes[j]<maxn;j++){
			num[i*primes[j]]=true;
			Mo[i*primes[j]]=Mo[i]*(-1);
			if(i%primes[j]==0){
				Mo[i*primes[j]]=0;
				break;
			}
		}
	}
}
inline int workk(int x,int y){
	x/=k,y/=k;
	if(x>y)swap(x,y);
	int ret=0;
	for(int i=1;i<=x;i++){
		ret+=(x/i)*(y/i)*Mo[i];
	}
	return ret;
}
int main(){
	Mobius();
	scanf("%d",&T);
	for(int i=1;i<=T;i++){
		scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
		ans=workk(b,d)-workk(a-1,d)-workk(b,c-1)+workk(a-1,c-1);
		printf("%d\n",ans);
	}
	return 0;
}