题意：给定a,b,c,d,k
x属于[1 , c]，y属于[1 , d]，求满足gcd(x,y)=k的对数。其中<x,y>和<y,x>算相同。
解法一：不妨设c<d,x<=y。问题可以转化为x属于[1，c / k ]，y属于[1，d/k ]，x和y互质的对数。
那么假如y<=c/k,那么对数就是y从1到c/k欧拉函数的和。如果y>c/k，就只能从[ c/k+1 ， d ]枚举，然后利用容斥。详见代码：
/********************************************************* file name: hdu1695.cpp author : kereo create time: 2015年02月11日 星期三 18时08分43秒*********************************************************/#include<iostream>#include<cstdio>#include<cstring>#include<queue>#include<set>#include<map>#include<vector>#include<stack>#include<cmath>#include<string>#include<algorithm>using namespace std;typedef long long ll;const int sigma_size=26;const int N=100+50;const int MAXN=100000+50;const int inf=0x3fffffff;const double eps=1e-8;const int mod=1000000000+7;#define L(x) (x<<1)#define R(x) (x<<1|1)#define PII pair<int, int>#define mk(x,y) make_pair((x),(y))int primecnt,factcnt;int prime[MAXN],euler[MAXN],factor[N][2];void getprime(){primecnt=0;memset(prime,0,sizeof(prime));for(int i=2;i<MAXN;i++){if(!prime[i])prime[primecnt++]=i;for(int j=0;j<primecnt && prime[j]<=MAXN/i;j++){prime[prime[j]*i]=1;if(i%prime[j] == 0)break;}}}void geteuler(){memset(euler,0,sizeof(euler));euler[1]=1;for(int i=2;i<MAXN;i++){if(!euler[i]){for(int j=i;j<MAXN;j+=i){if(!euler[j])euler[j]=j;euler[j]=euler[j]/i*(i-1);}}}}int getfactor(int x){factcnt=0;for(int i=0;prime[i]<=x/prime[i];i++){factor[factcnt][1]=0;if(x%prime[i] == 0){factor[factcnt][0]=prime[i];while(x%prime[i] == 0){x/=prime[i];factor[factcnt][1]++;}factcnt++;}}if(x!=1){factor[factcnt][0]=x;factor[factcnt++][1]++;}return factcnt;}int solve(int n,int m){int ans=0;getfactor(m);for(int i=1;i<(1<<factcnt);i++){int cnt=0,res=1;for(int j=0;j<factcnt;j++){if(i&(1<<j)){cnt++;res*=factor[j][0];}}if(cnt & 1)ans+=n/res;elseans-=n/res;}return n-ans;}int main(){int T,kase=0;scanf("%d",&T);getprime();geteuler();while(T–){printf("Case %d: ",++kase);int a,b,c,d,k;scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);if(k == 0 || k>b || k>d){printf("0\n");continue;}if(b>d)swap(b,d);b/=k; d/=k;ll ans=0;for(int i=1;i<=b;i++)ans+=euler[i];for(int i=b+1;i<=d;i++)ans+=solve(b,i);printf("%I64d\n",ans);}return 0;}
解法二：莫比乌斯反演。参见：
其中"设F(a,b,k)表示有多少组x≤a,y≤b，且Gcd(a,b)≥k"的"Gcd(a,b)>=k"应该是k | Gcd(x,y)。
/********************************************************* file name: hdu1695.cpp author : kereo create time: 2015年02月12日 星期四 09时08分41秒*********************************************************/#include<iostream>#include<cstdio>#include<cstring>#include<queue>#include<set>#include<map>#include<vector>#include<stack>#include<cmath>#include<string>#include<algorithm>using namespace std;typedef long long ll;const int sigma_size=26;const int N=100+50;const int MAXN=100000+50;const int inf=0x3fffffff;const double eps=1e-8;const int mod=1000000000+7;#define L(x) (x<<1)#define R(x) (x<<1|1)#define PII pair<int, int>#define mk(x,y) make_pair((x),(y))int primecnt;int vis[MAXN],mu[MAXN],prime[MAXN],sum[MAXN];void Mobius(){primecnt=0;memset(vis,0,sizeof(vis));mu[1]=1;for(int i=2;i<MAXN;i++){if(!vis[i]){prime[primecnt++]=i;mu[i]=-1;}for(int j=0;j<primecnt && i*prime[j]<MAXN;j++){vis[i*prime[j]]=1;if(i%prime[j])mu[i*prime[j]]=-mu[i];else{mu[i*prime[j]]=0;break;}}}}ll solve(int l,int r){ll ans=0;if(l>r)swap(l,r);int last;for(int i=1;i<=l;i=last+1){last=min(l/(l/i),r/(r/i));ans+=(ll)(sum[last]-sum[i-1])*(ll)(l/i)*(ll)(r/i);}return ans;}int main(){int T,kase=0;Mobius();sum[0]=0;for(int i=1;i<MAXN;i++)sum[i]=sum[i-1]+mu[i];scanf("%d",&T);while(T–){int a,b,c,d,k;printf("Case %d: ",++kase);scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);if(k == 0 || k>b || k>d){printf("0\n");continue;}if(b>d)swap(b,d);b/=k; d/=k;printf("%I64d\n",solve(b,d)-solve(b,b)/2);}return 0;}
，你并不一定会从此拥有更美好的人生，
