hdu 1695 GCD (欧拉函数+容斥原理+线性筛法)

题目：http://acm.hdu.edu.cn/showproblem.php?pid=1695

GCD

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

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

 

在1--b和1--d之间寻找最大公倍数是k的两个数的对数。小技巧：两个区间同时除以k则变成1--b/k和1--d/k中各挑选一个数看有多少对互质，简化了问题。因为1,3和3,1算作同一对，那么我们可以约定X<Y（同时在一开始就处理成b<d）。当Y<=b时就是欧拉函数问题。对正整数n，欧拉函数是少于或等于n的数中与n互质的数的数目【关于等于n的情况：φ(1)=1】那么Y>b/k时就得逐个素因子分解，利用容斥原理求出非互质的情况，再算出互质的数目。为了快速，可以先用线性筛法求出1---maxn的素数，然后再对每个大于b/k的数素因子分解。【线性筛法的时间复杂度：O(n)，埃拉托色尼筛法：O(nloglog(n))，下面用到的筛法求欧拉函数就是在埃拉托色尼筛法基础上实现的】

#include <iostream>
#include <cstdio>
#include <algorithm>
using namespace std;
const int maxn=1e6+10;
typedef long long LL;
LL prime[maxn],top;
bool notprime[maxn];
void getprime(){
    top=0;
    for(LL i=2;i<maxn;i++){
        if(!notprime[i])prime[top++]=i;
        for(LL j=0;j<top&&prime[j]*i<maxn;j++){
            notprime[prime[j]*i]=1;
            if(i%prime[j]==0) break;
        }
    }
}
LL factor[1005],faccnt;
void resolve(LL x){
    faccnt=0;
    for(int i=0;prime[i]*prime[i]<=x;i++){
        if(x%prime[i]==0){
            factor[faccnt++]=prime[i];
            while(x%prime[i]==0){
                x/=prime[i];
            }
        }
    }
    if(x>1)factor[faccnt++]=x;
}
int euler[100010];
void getEuler()
{
    euler[1] = 1;
    for(int i = 2;i <= 100000;i++)
        if(!euler[i])
            for(int j = i; j <= 100000;j += i)
            {
                if(!euler[j])
                    euler[j] = j;
                euler[j] = euler[j]/i*(i-1);
            }
}
int main()
{
    //freopen("cin.txt","r",stdin);
    getprime();
    getEuler();
    LL t,ca;
    LL a,b,c,d,k;
    LL kb,kd;
    cin>>t;
    for(ca=1;ca<=t;ca++){
         scanf("%lld%lld%lld%lld%lld",&a,&b,&c,&d,&k);
         if(k==0||k>b||k>d) {  //除以0会出现0x0000005c的错误
              printf("Case %lld: 0\n",ca);
              continue;
         }
         if(b>d) swap(b,d);
         b=b/k;
         d=d/k;
         LL ans=0;
         for(int i=1;i<=b;i++){
             ans+=euler[i];
         }
         for(int i=b+1;i<=d;i++){
             resolve(i);
             LL q1=0;
             for(int j=1;j<(1<<faccnt);j++){
                 LL temp=1,sum=0;
                 for(int k=0;k<faccnt;k++){
                     if((1<<k)&j){
                         sum++;
                         temp*=factor[k];
                     }
                 }
                 if(sum&1) q1=q1+b/temp;
                 else q1=q1-b/temp;
             }
             ans=ans+b-q1;
         }
         printf("Case %lld: %lld\n",ca,ans);
    }
    return 0;
}