HDOJ 1695 GCD（容斥+欧拉函数&&莫比乌斯反演+分块）

于 2017-06-28 19:15:52 发布
阅读量527 收藏 1
点赞数 1
分类专栏： 数论 文章标签： 容斥 欧拉函数 莫比乌斯反演 分块
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/my_sunshine26/article/details/73850529
版权


GCD

TimeLimit: 6000/3000 MS (Java/Others)    Memory Limit:32768/32768 K (Java/Others)
Total Submission(s): 11181    Accepted Submission(s): 4239

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 numberof choices may be very large, you're only required to output the total numberof 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 numberof 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 describedabove.

 

 

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).



题意：给出a,b,c,d,k,求满足 a<=x<=b && c<=y<=d && gcd(x,y)==k 的数对 （x,y）的对数（其中a=c=1） ，且（1，2）和（2，1）算同一种。


分析：由于求两个数的gcd==k 比较难求，我们需要把它转化为gcd==1的情况，那样就可以利用素数的性质去判断。

           转化过程为：原式变为 ：求满足 1<=x<=b/k && 1<=y<=d/k && gcd(x,y)==1 的数对 （x,y）的对数


方法一：容斥+欧拉函数

              由于这里有两个变量x,y;我们可以“定一变一”，比如枚举x的值，每次固定一个x值，求1<=y<=d/k &&gcd(x,y)==1的对数

              我们可以反向考虑，所有y的个数减去（x,y）不互质的对数，那么就是一个经典的容斥问题了：枚举x的质因子,然后根据容斥计算（1，y）中与x不互质的数的个数。

              去重：举个例子,当x=2&&1<=y<=5时 ，y可取 1, 3,5

                                        当x=3&&1<=y<=5时， y可取 1,2,4,5

                                        重复的是(2，3) 和（3，2） ,易知我们只要控制x<y即可，也就是说上面的结果减去（1，y）中与x互质且小于x的数，用欧拉函数即可解决

           PS：要特判k==0(除0错误)以及k>b||k>d

代码(1076MS)

#include <cstdio>
#include <cmath>
#include <stack>
#include <cstring>
#include <algorithm>
using namespace std;
#define mst(a,b) memset((a),(b),sizeof(a))
#define f(i,a,b) for(int i=(a);i<(b);++i)
#define rush() int T;scanf("%d",&T);while(T--)
typedef long long ll;
const int maxn= 35;
const ll mod = 1e9+7;
const int INF = 0x3f3f3f3f;
const double eps = 1e-6;

int prime[maxn];
int m;

ll euler(ll num)       //既求出euler(num),又求出num的质因子
{
    m=0;
    ll res=num;
    for(ll i=2; i*i<=num; i++)
    {
        if(num%i==0)
        {
            prime[m++]=i;
            res=res/i*(i-1);
            while(num%i==0)
                num/=i;
        }
    }
    if(num>1)
    {
        res=res/num*(num-1);
        prime[m++]=num;
    }
    return res;
}

ll solve(ll num)
{
    ll ans=0,temp;
    int flag;
    for(ll i=1; i<(ll)(1<<m); i++)
    {
        temp=1;
        flag=0;
        for(int j=0; j<m; j++)
        {
            if(i&(ll)(1<<j))
            {
                flag++;
                temp*=prime[j];
            }
        }
        if(flag&1)
            ans+=num/temp;
        else ans-=num/temp;
    }
    return ans;
}

int main()
{
    ll a,b,c,d,k;
    int cas=1;
    rush()
    {
        scanf("%I64d%I64d%I64d%I64d%I64d",&a,&b,&c,&d,&k);
        if(k==0||k>b||k>d)
        {
            printf("Case %d: 0\n",cas++);
            continue;
        }
        if(b>d)              //优化
            swap(b,d); 
        ll ans=0;
        for(ll i=k; i<=b; i+=k)
        {
            ll temp=euler(i/k);
            ans+=(d/k-solve(d/k)-temp);
        }
        printf("Case %d: %I64d\n",cas++,ans+1);
    }
    return 0;
}

方法二：莫比乌斯反演 （刚入门，感觉在解决复杂问题及时间优化上具有很大优势）

             

先贴上莫比乌斯反演的公式及其性质


公式一：

公式二：

证明：


性质：



求莫比乌斯函数的模板（同时筛出了区间内的素数）

void init()
{
    mu[1]=1;
    for(int i=2;i<maxn;i++)
    {
        if(!isprime[i])
        {
            prime[num_prime++]=i;
            mu[i]=-1;
        }
        for(int j=0;j<num_prime&&i*prime[j]<maxn;j++)
        {
            isprime[i*prime[j]]=1;
            if(i%prime[j]==0)
            {
                mu[i*prime[j]]=0;
                break;
            }
            else mu[i*prime[j]]=-mu[i];
        }
    }
}



分析：设f(k)为满足gcd（x,y）=k的（x,y）对数 ，那么我们要求的便是f（1）；

           设F(k)为满足gcd（x,y）=k的倍数 的（x,y）对数，易知F（k）=floor（b/k）*floor（d/k）。

           由莫比乌斯反演法的公式一可得 f(1)=mu[1]*F(1) + mu[2]*F[2] + ... + mu[max]*F(max)  其中max=min(b/k,d/k)

           同样，这里也需要去重： 假设1<=x<=3 , 1<=y<=5, 列举两者的结果，发现重复的有 （1，1），（1，2），（1，3），（2，3）和 （1，1），（2，1），（3，1），（3，2） 即区间【1，3】内的结果 ，区间的右边取决于min（x,y），减去这一部分的结果的1/2即可


代码（15MS） 

#include <cstdio>
#include <cmath>
#include <queue>
#include <stack>
#include <cstring>
#include <algorithm>
using namespace std;
#define mst(a,b) memset((a),(b),sizeof(a))
#define f(i,a,b) for(int i=(a);i<(b);++i)
#define rush() int T;scanf("%d",&T);while(T--)
typedef long long ll;
const int maxn= 100005;
const ll mod = 1e9+7;
const int INF = 0x3f3f3f3f;
const double eps = 1e-6;

int prime[maxn]= {0};
int num_prime=0;
int mu[maxn];
bool isprime[maxn]= {1,1};


void init()
{
    mu[1]=1;
    for(int i=2;i<maxn;i++)
    {
        if(!isprime[i])
        {
            prime[num_prime++]=i;
            mu[i]=-1;
        }
        for(int j=0;j<num_prime&&i*prime[j]<maxn;j++)
        {
            isprime[i*prime[j]]=1;
            if(i%prime[j]==0)
            {
                mu[i*prime[j]]=0;
                break;
            }
            else mu[i*prime[j]]=-mu[i];
        }
    }
}

int main()
{
    int cas=1;
    int a,b,c,d,k;
    init();
    rush()
    {
        scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
        if(k==0)
        {
            printf("Case %d: 0\n",cas++);
            continue;
        }
        b/=k,d/=k;
        if(b>d)
            swap(b,d);
        ll ans1=0;
        ll ans2=0;
        for(int i=1;i<=b;i++)
        {
            ans1+=(ll)mu[i]*(b/i)*(d/i);
            ans2+=(ll)mu[i]*(b/i)*(b/i);
        }
        ans1-=ans2/2;
        printf("Case %d: %I64d\n",cas++,ans1);
    }
    return 0;
}

方法三：莫比乌斯反演+分块


上面的方法每次查询时间复杂度为O(n)，那么有没有比这更快的解法呢？实际上是有的，利用分块可以是实现O(sqrt(n))的时间复杂度
实现思路：先来看看上面的求法，其实就是枚举公约数，一个个加上去，但考虑到 b/i 和 d/i 不能整除的特点，我们很容易发现floor（b/i）和floor（d/i）在一定的区间内的值是不变的，也就是说这里可以实现分块加速，举个例子，当b = 10，d  =  11时， i从6到10，l / i和r / i的值都是1。
实现主干代码见下 

ll ans=0;
for(int i=1; i<=n; i=last+1)
{
    last=min(n/(n/i),m/(m/i));
    ans+=(ll)(sum[last]-sum[i-1])*(n/i)*(m/i);
}


完整代码（0MS）


#include <cstdio>
#include <cmath>
#include <queue>
#include <stack>
#include <cstring>
#include <algorithm>
using namespace std;
#define mst(a,b) memset((a),(b),sizeof(a))
#define f(i,a,b) for(int i=(a);i<(b);++i)
#define rush() int T;scanf("%d",&T);while(T--)
typedef long long ll;
const int maxn= 100005;
const ll mod = 1e9+7;
const int INF = 0x3f3f3f3f;
const double eps = 1e-6;

int prime[maxn]= {0};
int num_prime=0;
int mu[maxn];
bool isprime[maxn]= {1,1};
int sum[maxn];
int a,b,c,d,k;


void init()
{
    mu[1]=1;
    for(int i=2; i<maxn; i++)
    {
        if(!isprime[i])
        {
            prime[num_prime++]=i;
            mu[i]=-1;
        }
        for(int j=0; j<num_prime&&i*prime[j]<maxn; j++)
        {
            isprime[i*prime[j]]=1;
            if(i%prime[j]==0)
            {
                mu[i*prime[j]]=0;
                break;
            }
            else mu[i*prime[j]]=-mu[i];
        }
    }
    sum[0]=0;
    for(int i=1; i<maxn; i++)
    {
        sum[i]=sum[i-1]+mu[i];
    }
}

ll solve(int x,int y)
{
    int n=x/k,m=y/k;
    int last;
    if(n>m)  swap(n,m);
    ll ans=0;
    for(int i=1; i<=n; i=last+1)
    {
        last=min(n/(n/i),m/(m/i));
        ans+=(ll)(sum[last]-sum[i-1])*(n/i)*(m/i);
    }
    return ans;
}

int main()
{
    int cas=1;
    init();
    rush()
    {
        scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
        if(k==0)
        {
            printf("Case %d: 0\n",cas++);
            continue;
        }
        ll ans1=solve(b,d);
        if(b>d) swap(b,d);
        ans1-=solve(b,b)/2;
        printf("Case %d: %I64d\n",cas++,ans1);
    }
    return 0;
}






















                

        

        

    




    

      

        

            

              
                
                  Dust_Heart
                
              
            

            

                关注
              
            

        

        

          
    
            
  • 
              
                
                
                
                
                    1
                
              
              
    点赞
    
            
    • 
            
  • 
              
                
                
                
              
              
    踩
    
            
    • 
            
  • 
              
                
                
                
                
                    1
                
              
              
    
                
    
                  
                    收藏
                  
                
    
              
    
              
    
                
    
                  觉得还不错?
                  
                    一键收藏
                  
                 
                
    
              
    
            
    • 
          
  • 
            
    
              
              
            
    
              
              
              
                    0
              
            
            
    评论
    
          
    • 
          
  • 
          
    • 
          
  • 
            
              
            
              
    
            
              
                
            
    
          
    • 
        

      

      

            

                专栏目录
          










              

  



          




        



    





    



      

      

        
      

      





	

		评论	
	

	
	




  

    

      添加红包
      
    

    

      

        
        

          
          
        

        
请填写红包祝福语或标题

      

      

        
        

          
          个
        

        
红包个数最小为10个

      

      

        
        

          
          元
        

        
红包金额最低5元

      

      

        
        

          当前余额3.43元
          前往充值 >
        

      

      

        

          需支付：10.00元
        

        
        
      

    

  





  

    
    
  

  

    
    
  








  

    

      

        

          

            
            
成就一亿技术人!

          

          

        

        

          

          

            领取后你会自动成为博主和红包主的粉丝
            规则
          

        

      

      

        

          

            
              
            
            

              
hope_wisdom
 发出的红包
            

          

        

        

          

          

          

        

      

    

    

  



      
      

      
实付元

      
使用余额支付

      

        

          

            
            点击重新获取
          

        

        
扫码支付

      

      

        
          
            
          
          
        
      

      

        
        钱包余额
          0
          

            
            

              

                
抵扣说明：

                
 1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
 2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

              

            

          

      

      余额充值