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;
}
posted @ 2019-01-18 17:22  虚拟北方virtual_north。  阅读(101)  评论(0编辑  收藏  举报