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;
}
