HDU 1695 GCD 容斥原理/莫比乌斯反演

题意：

给你两个区间[a,b],[c,d]，还有一个k。让你从区间[a,b]中找出x，[c,d]中找出y，问共有多少组（x，y）使得gcd（x，y）=k。

（x，y）和（y，x）算一组。

思路：

参考：http://blog.csdn.net/yang_7_46/article/details/9072533

容斥。

普通容斥：

*如果gcd（x，y）=k，则gcd（x/k，y/k）=1。那么对于两个区间来说，我们都默认进行了b = b/k，d = d/k操作。（a，c固定为1）。（这样做并不会使得组数遗漏）

因为（x，y）和（y，x）只算一组，因此在容斥过程中，我们要保证x < y，这样才不会重复计算。

剩下就是容斥的过程，可以戳这里。

其它一些细节要自己处理一下，例如b = 1或者d = 1的时候。

莫比乌斯反演：

莫比乌斯资料：http://blog.csdn.net/acdreamers/article/details/8542292

定义f(n)：gcd(x, y)为n的方案数。

定义F(n)：gcd(x, y)是n的倍数的方案数。

则我们要求的就是f(1)。默认b < d，并且都已经除以k。

F(n) = (b/n)*(d/n);

套用莫比乌斯公式即可。

code（普通容斥）：

#include <cstring>
#include <cstdlib>
#include <cstdio>
#include <iostream>
#include <algorithm>
#include <vector>
using namespace std;

const int MAXN = 1e5+5;
typedef long long LL;

int a, b, c, d, k;
bool isp[MAXN+5];
vector <int> vec[MAXN+5];

void prime()
{
    memset(isp, false, sizeof(isp));
    for(int i = 2;i <= MAXN; i++)
    {
        if(!isp[i])
        {
            vec[i].push_back(i);
            for(int j = i*2; j <= MAXN; j += i)
            {
                isp[j] = true;
                vec[j].push_back(i);
            }
        }
    }
}
LL calc(int t, int p)
{
    LL ret, v = 1;
    int cnt = 0;
    for(int i = 0;i < vec[t].size(); i++)
    {
        if((1<<i)&p)
        {
            cnt++;
            v *= vec[t][i];
        }
    }
    if(v == 0) return 0;
    ret = (LL)(d-t)/v;
    if(cnt%2 == 0) ret = -ret;
    return ret;
}

void solve()
{
    //tepan 1
    LL res = d;
    for(int i = 2;i <= b; i++)
    {
        res += d-i;
        for(int j = 1;j < (1<<vec[i].size()); j++)
            res -= calc(i, j);
    }
    printf("%I64d\n", res);
}

int main()
{
    prime();
    int T, cas = 0;
    scanf("%d", &T);
    while(T--)
    {
        scanf("%d%d%d%d%d", &a, &b, &c, &d, &k);
        printf("Case %d: ", ++cas);
        if(k == 0)
        {
            puts("0");
            continue;
        }
        b /= k, d /= k;
        if(b > d) swap(b, d);
        if(b == 0 || d == 0)
            puts("0");
        else 
            solve();
    }
    return 0;
}

 code（莫比乌斯反演）：

#include <bits/stdc++.h>
using namespace std;

const int N = 1e5+5;
typedef long long LL;

int b, d, k;
int prime[N], cnt;
int mu[N];
bool vis[N];
void Mobius() {
    mu[1] = 1;
    cnt = 0;
    for(int i = 2;i < N; i++) {
        if(!vis[i]) {
            mu[i] = -1;
            prime[cnt++] = i;
        }
        for(int j = 0;j < cnt; j++) {
            if(i*prime[j] >= N) break;
            vis[i*prime[j]] = true;
            if(i%prime[j] != 0)
                mu[i*prime[j]] = -mu[i];
            else {
                mu[i*prime[j]] = 0;
                break;
            }
        }
    }
}

int main() {
    Mobius();
    int T, cas = 0;
    scanf("%d", &T);
    while(T--) {
        scanf("%*d%d%*d%d%d", &b, &d, &k);
        printf("Case %d: ", ++cas);
        if(k == 0) {
            puts("0");
            continue;
        }
        b /= k, d /= k;
        if(b > d) swap(b, d);
        LL ans = 0;
        for(int i = 1;i <= b; i++)
            ans += 1ll*mu[i]*(b/i)*(d/i);
        LL tmp = 0;
        for(int i = 1;i <= b; i++)
            tmp += 1ll*mu[i]*(b/i)*(b/i);
        printf("%I64d\n", ans-tmp/2);
    }
    return 0;
}