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
",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
",ans);
}
return 0;
}