HDU 1695 GCD 欧拉函数+容斥定理 || 莫比乌斯反演

GCD

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

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

/*
题意：区间x属于[1,A] ， y属于区间[1,B]
     求最大公约数是K，即gcd(x,y)=K。
     并且[1,3]和[3,1]属于同一种情况。
     
思路：HDU 4135 Co-prime 的思路在这一题有用。
     它的题意：区间[A,B],与整数N的互素的个数
     对于这一到题目：gcd(x,y)=k.
     要满足最大公约数是K，可以转化为
     [1,A],[1,B]==>[1,A/K],[1,B/K] 求互素的个数。
     好像有点难以想到。？？？？
     {{
     借鉴一下别人是说法。会更明白
     gcd(x, y) == k 说明x,y都能被k整除,但是能被k整除的未必gcd=k , 
     必须还要满足互质关系. 
     问题就转化为了求1~a/k 和 1~b/k间互质对数的问题
     }}
     这样的话，如何处理呢？
     题意要求[1,3]和[3,1]不能重复。
     对于区间[1,A/K],[1,B/K] 看成==>[1,a],[1,b] 有几种情况
     1____________a
     1____________________b
     
     
     1____________a
     1________b
     
     
     1____________a
     1____________b
     
     这三种情况。我们来个判断，总是让a<=b，用b做更大的值。就会变成
     
     1—————————a
     1—————————————————b
     在求取的过程中也是采取这样的规则。
     [?,b1];确定后一位数。表示在[1,a]中与b1互质的个数。
     那么就很好的避免了[1,3],[3,1]的情况了。
求取总和sum=sum1+sum2；
sum1=欧拉函数值[1,a]; 想想为什么？
sum2={枚举a+1--->b,与区间[1,a]互质的个数}; 
sum2就和以前的一题有关系了，要用欧拉函数+容斥定理处理。
具体的参考：http://www.cnblogs.com/tom987690183/p/3246197.html


*/

#include<stdio.h>
#include<string.h>
#include<stdlib.h>


int prime[100003],len;
bool s[100003];
int opl[100003];
int Que[100003];
int f[1000],flen;

void make_prime() //素数打表
{
   int i,j;
   len=0;
   for(i=2;i<=100000;i++)//刚开始写错，i*i<=100000;⊙﹏⊙b汗
   if(s[i]==false)
   {
       prime[++len]=i;
       for(j=i*2;j<=100000;j=j+i)
       s[j]=true;
   }
}

void make_Euler() //欧拉函数打表。
{
   int i,j;
   make_prime();
   for(i=2;i<=100000;i++)
   opl[i]=i;
   opl[1]=1;
   for(i=1;i<=len;i++)
   for(j=prime[i];j<=100000;j=j+prime[i])
   opl[j]=opl[j]/prime[i]*(prime[i]-1);
}

void make_dEuler(int n) //单点欧拉的素因子。
{
   int i;
   flen=0;
   for(i=2;i*i<=n;i++)
   {
       if(n%i==0)
       {
           while(n%i==0)
           n=n/i;
           f[++flen]=i;
       }
   }
   if(n!=1)
   f[++flen]=n;
}

int Capacity(int m)
{
   int i,j,t=0,sum=0,k;
   Que[t++]=-1;
   for(i=1;i<=flen;i++)
   {
       k=t;
       for(j=0;j<k;j++)
       Que[t++]=-1*Que[j]*f[i];
   }
   for(i=1;i<t;i++)
   sum=sum+m/Que[i];
   return sum;
}


void sc()//输出函数，测试用的。
{
   int i;
   for(i=1;i<=10;i++)
   printf("%d ",opl[i]);
   printf("
");
}

__int64 make_ini(int b,int c,int k)
{
   int i,x,y,tmp;
   __int64 sum=0;
   x=b/k;y=c/k;//加特判的用处。不能除0
   if(x>y)
   {
       tmp=x;
       x=y;
       y=tmp;
   }
   for(i=1;i<=x;i++)
   sum=sum+opl[i];//第一步
   for(i=x+1;i<=y;i++)//第二步，枚举
   {
       make_dEuler(i);
       sum=sum+(x-Capacity(x));
   }
   //sc();
   return sum;

}

int main()
{
   int T,a,b,c,d,k,i;
   __int64 sum;
   make_Euler();
   while(scanf("%d",&T)>0)
   {
       for(i=1;i<=T;i++)
       {
           sum=0;
           scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
           if(k==0) //特判，否则会Runtime Error (INTEGER_DIVIDE_BY_ZERO)
           {
               sum=0;
           }
           else sum=make_ini(b,d,k);
           printf("Case %d: %I64d
",i,sum);

       }
   }
   return 0;
}

下面再介绍一种方法。莫比乌斯反演

GCD(a,b) = d;

可以转化为

GCD(a/d,b/d) = 1;

设f(d)为(a,b) = d的种类数

F(d)为(a,b) = d 的倍数的种类数。

例如

F(2) = (a/2)*(b/2);

F(3) = (a/3)*(b/3);

由

mu[i]可以打表求出。

关于一个优化在于a/i = a/(i+k) && b/i = b/(i+k);

此时我们能节省时间来求。详见代码部分

 1 #include<iostream>
 2 #include<stdio.h>
 3 #include<cstring>
 4 #include<cstdlib>
 5 using namespace std;
 6 typedef __int64 LL;
 7 
 8 const int maxn = 1e5+3;
 9 bool s[maxn];
10 int prime[maxn],len = 0;
11 int mu[maxn];
12 int sum1[maxn];
13 void  init()
14 {
15     memset(s,true,sizeof(s));
16     mu[1] = 1;
17     for(int i=2;i<maxn;i++)
18     {
19         if(s[i] == true)
20         {
21             prime[++len]  = i;
22             mu[i] = -1;
23         }
24         for(int j=1;j<=len && (long long)prime[j]*i<maxn;j++)
25         {
26             s[i*prime[j]] = false;
27             if(i%prime[j]!=0)
28                 mu[i*prime[j]] = -mu[i];
29             else
30             {
31                 mu[i*prime[j]] = 0;
32                 break;
33             }
34         }
35     }
36     for(int i=1;i<maxn;i++)
37         sum1[i] = sum1[i-1]+mu[i];
38 }
39 LL solve(int a,int b)
40 {
41     LL sum = 0;
42     for(int i=1,la = 0;i<=a;i++,i = la+1)
43     {
44         la = min(a/(a/i),b/(b/i)); //优化部分
45         sum = sum + ((LL)(a/i))*(b/i)*(sum1[la]-sum1[i-1]);
46     }
47     return sum;
48 }
49 int main()
50 {
51     int T,l,a,b,d;
52     init();
53     scanf("%d",&T);
54     for(int t=1;t<=T;t++)
55     {
56         scanf("%d%d%d%d%d",&l,&a,&l,&b,&d);
57         LL sum = 0 ;
58         if(d==0) ;
59         else{
60             if(a>b) swap(a,b);
61             sum = solve(a/d,b/d);
62             sum = sum - solve(a/d,a/d)/2;
63         }
64         printf("Case %d: %I64d
",t,sum);
65     }
66     return 0;
67 }