hdu 1695 GCD

GCD

Time Limit: 6000/3000 MS (Java/Others)    Memory Limit: 32768/32768 K (Java/Others)
Total Submission(s): 8756    Accepted Submission(s): 3248

Problem Description

Given 5 integers: a, b, c, d, k, you're to find x in a...b, y in c...d that GCD(x, y) = k. GCD(x, y) means the greatest common divisor of x and y. Since the number of choices may be very large, you're only required to output the total number of different number pairs.
Please notice that, (x=5, y=7) and (x=7, y=5) are considered to be the same.

Yoiu can assume that a = c = 1 in all test cases.

 

Input

The input consists of several test cases. The first line of the input is the number of the cases. There are no more than 3,000 cases.
Each case contains five integers: a, b, c, d, k, 0 < a <= b <= 100,000, 0 < c <= d <= 100,000, 0 <= k <= 100,000, as described above.

 

Output

For each test case, print the number of choices. Use the format in the example.

 

Sample Input

  
  

   
   2
1 3 1 5 1
1 11014 1 14409 9
  
  

 

Sample Output

  
  

   
   Case 1: 9
Case 2: 736427


   
   

    
    

     
     Hint
    
    
For the first sample input, all the 9 pairs of numbers are (1, 1), (1, 2), (1, 3), (1, 4), (1, 5), (2, 3), (2, 5), (3, 4), (3, 5).
   
   
   
    
  
  

 

Source

2008 “Sunline Cup” National Invitational Contest

 

Recommend

wangye   |   We have carefully selected several similar problems for you:   1689  1690  1693  1691  1698 

题目大意：sigma(i=a..b) sigma(j=c..d) [gcd(i,j)==k]  (注意a=c=1, gcd(x,y)=gcd(y,x)=k 只计算一次）

另外最后的结果需要用long long 

题解:  sigma(i=a..b) sigma(j=c..d) [gcd(i,j)==k] 

         =sigma(i=a..b)sigma(j=c..d) [gcd(i/k,j/k)==1]

       设i=i'*k ,j=j'*k   b'=b/k  d'=d/k

        =sigma(i=1..b')sigma(j=1..d')  [gcd(i,j)==1]

        =sigma(i=1..b')sigma(j=1..d')sigma(k|(i,j)) mu(k)

        =sigma(i=1..b')sigma(j=1..d')sigma(k=1..min(b',d')) [k|i] [k|j] mu[k]

        =sigma(k=1..min(b',d')) (b'/k)*(d'/k)*mu[k]

这样时间复杂度就很科学了，就可以搞了。

那么如何去重呢？ 我们发现重复的只可能在b,d 中较小范围内，即只可能是t=min(b,d)  sigma(i=1..t) sigma(j=1..t) [gcd(i,j)==k] 

那么最后的答案就是   

sigma(k=1..min(b',d')) (b'/k)*(d'/k)*mu[k] -   sigma(k=1..min(b',d'))  (min(b',d')/k)*(min(b',d')/k)*mu[k]/2   (因为这么算也是重复的，所以需要除以2）

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
int n,m,t;
int phi[100003],mu[100003],p[100003],prime[100003];
int a,b,c,d,nn,mm,k;
long long ans,ans1;
void calc()
{
  p[1]=1; phi[1]=1; mu[1]=1;
  for (int i=2;i<=100000;i++)
   {
   	if (!p[i])
   	mu[i]=-1, phi[i]=i-1, prime[++prime[0]]=i;
    for(int j=1;j<=prime[0];j++)
    {
     if (prime[j]*i>100000) break;
     p[prime[j]*i]=1;
     if (i%prime[j]==0)
     {
      mu[i*prime[j]]=0;
      phi[i*prime[j]]=phi[i]*prime[j];
      break;
     }
     else
     {
      mu[i*prime[j]]=-mu[i];
      phi[i*prime[j]]=phi[i]*(prime[j]-1);
     }
    }
   }
}
int main()
{
  scanf("%d",&t);
  calc();
  for (int i=1;i<=t;i++)
  {
   scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
   if (!k||!b||!d)
    {
      printf("Case %d: 0\n",i);
      continue;
    }
   nn=b/k; mm=d/k; ans=0; ans1=0;
   int t=min(nn,mm);
   for (int j=1;j<=t;j++)
    ans+=(long long)(nn/j)*(mm/j)*mu[j];
   for (int j=1;j<=t;j++)
    ans1+=(long long)(t/j)*(t/j)*mu[j];
   ans1/=2;
   printf("Case %d: %I64d\n",i,ans-ans1);
  }
}

这道题虽然可以用上面那种比较暴力的方式水过，但是显然有更加优越的做法，就是利用莫比乌斯反演。

【简单的莫比乌斯反演题。莫比乌斯反演可以说是对函数的一种变形，利用Dirichlet卷积进行有某种关系的函数之间的相互转化，可以简化运算。它的基本公式有两种形式，一是已知F(n)=sigma_(d|n)f(n)，则f(n)=sigma_(d|n) μ(n/d)F(d)。其中μ(i)是莫比乌斯函数并且是一个积性函数，当i为质数时值为-1，当i的最小质因子次数大于1时值为0。以上公式化成卷积的形式就是已知F=f×1，则f=μ×F。这里的1是将所有实数映射到1的一个函数，说白了就是所有自变量对应的函数值都是1。由定义式可以看出卷积是满足交换律和结合律的。以上结论的证明如下：】

【<1>.已知F=f×1，则等式两边同卷μ，得F×μ=f×1×μ。】

【<2>.而1×μ=e，e函数的性质是：当自变量为1时函数值为1，否则为0。】

【<3.>则F×μ=f×e。等式两边看作函数，右边展开可以得到式子sigma_(d|n)f(d)e(n/d)，可以看出若自变量为n，则只有d=n时e的值为1，会累加一个f(n)的答案，所以f×e=f。即f=F×μ。】

【上面的公式是枚举自变量n的因数，而枚举自变量n的倍数也可以得到一个公式，证明是类似的。这道题要求sigma_(i=1..b)sigma_(j=1..d)[(i,j)=k]，公式变形如下：】

【<1>.求gcd的式子两边同时除以k变成：sigma_(i=1..b)sigma_(j=1..d)[k|i][k|j][(i/k,j/k)=1]。】

【<2>.直接枚举k的倍数，设i’=i*k，j’=j*k，替换可得：sigma_{(i’=1..b/k)}sigma_{(j’=1..d/k)}[(i’,j’)=1]。】

【于是对于这道题，我们可以设f(n)为最大公约数为n的数对个数，F(n)为公约数为n的数对个数，即最大公约数为n的倍数的数对个数，则f(1)就是我们要求的答案。于是显然存在关系F(n)=sigma_(n|d)f(d)。利用莫比乌斯反演的第二个公式可以得到f(n)=sigma_(n|d) μ(d)F(d/n)。接下来的问题就是求解F(n)，也就是求b和d范围内可以同时被n整除的数对个数。而我们如果知道[1..b]内可以被n整除的数的个数为p，[1..d]内可以被n整除的数的个数为q，利用乘法原理显然有F(n)=p*q，而利用结论：p=floor(b/n)，q=floor(d/n)，我们可以在O(1)的时间内得到任意的F值。那么我们如果要求解f(n)，那么枚举n的倍数d就可以了。而枚举的d看起来是没有上限的，实际上在计算F(d)时，如果d大于min{b,d}，那么整除的结果就会是0，就不会再累加答案了，所以我们最大枚举到min{b,d}就可以了。】

【我们把b和d中较小的赋给b，较大的赋给d，于是可以通过枚举1的倍数即全部i=1..b求出f(1)。这里有一个问题就是由F函数的定义可以注意到我们把相同的数对例如(2,3)和(3,2)重复计算了，并且重复计算的部分只有[1..b]范围内的那些数对，因为[1..d]范围内某个数的倍数个数显然是大于[1..b]范围内的个数的，所以[b+1..d]范围内的那些数只被枚举了一次。所以我们需要计算出[1..b]范围内不重复的满足条件的数对的个数，再用刚才算出的f(1)减去它。而这就相当于d和b相等，再用同样的方法计算一遍f(1)即可。因为两个上界的范围相同，都是b，所以所有数对都被计算了两次，可以直接除以2来获得不重复的答案。时间复杂度是O(n)的，看起来好像有点不科学实际上还是跑的很快的。】（以上证明过程感谢ATP神犇的援助><)

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#define N 100003
#define ll long long
using namespace std;
int mu[N],prime[N],pd[N];
int n,m,t,a,b,c,d,k;
long long ans;
void calc()
{
	mu[1]=1;
	for (int i=2;i<=N;i++)
	 {
	 	if (!pd[i])
	 	 {
	 	 	prime[++prime[0]]=i; mu[i]=-1;
	 	 }
	 	for (int j=1;j<=prime[0];j++)
	 	 {
	 	 	if (i*prime[j]>N) break;
	 	 	pd[i*prime[j]]=1;
	 	 	if (i%prime[j]==0)
	 	 	 {
	 	 	 	mu[i*prime[j]]=0;
	 	 	 	break;
	 	 	 }
	 	 	else
	 	 	 mu[i*prime[j]]=-mu[i];
	 	 }
	 }
}
ll F(int x)
{
	return (ll)floor(b/x)*floor(d/x);
}
int main()
{
	scanf("%d",&t);
	calc();
	for (int T=1;T<=t;T++)
	 {
	 	scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
	 	ans=0; 
	 	if (!k)
	 	 {
	 	 	printf("Case %d: %d\n",T,0);
	 	 	continue;
	 	 }
	 	b/=k; d/=k;
	 	if (b>d) swap(b,d);
	 	for (int i=1;i<=b;i++)
	 	 ans+=(ll)mu[i]*F(i);
	 	ll t=0;
	    d=b;
	 	for (int i=1;i<=b;i++)
	 	 t+=(ll)mu[i]*F(i);
	 	ans-=t/2;
	 	printf("Case %d: %lld\n",T,ans);
	 }
}