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

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;
}