专题练习---（数论）莫比乌斯反演

　　　　　　　　　　　　GCD

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

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

题意：  给出a,b,c,d,k, 使得  a<=x<=b , c<=y<=d  &&gcd(x,y)=k的数对（x,y）的对数。

　　限制：  a=c=1 ,  0<b,c<=100000 ;(n1,n2)和（n2 ,n1）算为同一种情况、

---------------------------------分割线-----------------------------------------------

个人理解：  莫比乌斯反演就是一种类似的容斥原理，首先来看看莫比乌斯函数。

   莫比乌斯函数：

　                    mu(d) =　　1　　　　　　d==1

　　　　　　　　mu(d) =　　(-1)^k         d=p1p2p3....pk (pk为素数)

　　　　　　　　mu(d) =　　0                其他 　　　　　　　

 莫比乌斯反演：

         求和函数：    

   其中  d|n 表示  n%d==0   1<=d<=n

         对于上面的公式，举例：

         F(1) = f(1)                         F(6) = f(1) +f(2) +f(3) + f(6)

         F(2) = f(1) +f(2)                 F(7) = f(1) + f(7)

　　　 F(3) = f(1) + f(3)               F(8) = f(1) +f(2)+f(4)+f(8)

　　　 F(4) = f(1) +f(2) +f(4)        F(9) = f(1)+f(3)+f(9)

　　　 F(5) = f(1) + f(5) 　　　　　F(10) = f(1)+f(2)+f(5)+f(10)

       实际中，往往是，我们求解F(N)比较容易，而求解f(x) 比较困难。但是我们亦反推求解f(x)

        f(1) = F(1)　　　　　　　　　f(6) = F(6) -F(3) -F(2) +F(1)

　　　f(2) = F(2) -F(1)                 f(7) = F(7) -F(1)

        f(3) = F(3) - F(1)　　　　　　f(8) = F(8) - F(4)

        f(4) = F(4) -F(2)　　　　　　 f(9) = F(9) - F(3)

        f(5) = F(5) -F(1) 　　　　　　f(10) = F(10) - F(5) -F(2) +F(1)

　　对于上述，我们不难总结出，已知F(x)， 求解 f(x)的公式：

       d|n 意思是   n%d==0

　　而且我们可以很欣喜的发现，U(d)其实就是我们上上面求解的mu(d)（即莫比乌斯函数）

　　其实对于上面的公式，我们如果仔细器观察，也不难发现，这个和容斥原理还是有那么点类似的。

   对于容斥原理

   　　    P(AUBUC) =P(A)+P(B)+P(C)-P(AUB)-P(BUC)-P(AUC)+P(ANBNC) ;

　　而如  f(6) = F(6) -F(3) -F(2) +F(1)  ，这样来看，还是是挺相似的。

----------------------------------------分割线-----------------------------------------------

  那么对于这道题，我们可以这样分析：

1  设f(k)为gcd(x,y)=k的数对(x,y)的对数，我们要求的是f(1)
2 设F(k)为gcd(x,y)为k的倍数的数对(x,y)的对数，可以想到F(k)=floor(b/k)*floor(d/k)，
   F(k)=E[b/k]*[d/k]*mu[k];

3 由莫比乌斯反演得：
4 令lim=min(b/k,d/k)
5 f(1)=mu[1]*F(1) + mu[2]*F[2] + ... + mu[lim]*F(lim)
6 因为(n1,n2)和(n2,n1)算为同一种情况，所以最后结果还要减掉重复的情况。

 代码:

 1 #include<cstdio>
 2 #include<cstdlib>
 3 #include<cstring>
 4 #include<iostream>
 5 #include<algorithm>
 6 #define _llint long long
 7 using namespace std;
 8 const int maxn=100010;
 9 bool jud[maxn];
10 int mu[maxn] , prime[maxn];
11 
12 void Mobius(){
13   int cnt=0;
14   memset(jud , 0 , sizeof jud);
15   mu[1]=1;   //d=1
16   for(int i=2 ; i<maxn ;i++){
17 
18       if(!jud[i]){
19           mu[i]=-mu[1];  //1*p
20           prime[cnt++]=i;
21       }
22     for(int j=0 ; j<cnt&&i*prime[j]<maxn ; j++ ){
23 
24         jud[i*prime[j]]=1;
25         if(i%prime[j]==0){
26             mu[i*prime[j]]=0 ; //²»ÊÇ»¥ÒìµÄËØÊý
27             break;
28         }else{
29              mu[i*prime[j]]=  -1*mu[i];
30         }
31     }
32   }
33   return ;
34 }
35 
36 _llint solve(int n, int m){
37 
38    _llint res=0;
39    for(int i=1; i<=n ;i++){
40      res+=(_llint)(m/i)*(n/i)*mu[i];
41    }
42    return res;
43 }
44 
45 int main()
46 {
47   int T,a,b,c,d,k;
48   Mobius();
49   scanf("%d",&T);
50   for(int i=1 ;i<=T ;i++){
51     scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
52     if(0==k)
53     {
54         printf("Case %d: 0\n",i);
55         continue;
56     }
57      b/=k  , d/=k;
58      if(b>d)    swap( b , d );
59     _llint ans1=solve(b,d);
60     _llint ans2=solve(b,b);
61      printf("Case %d: %lld\n",i,ans1-ans2/2);
62   }
63  return 0;
64 }