HDU 1695 GCD(容斥 or 莫比乌斯反演)

这题可以用容斥做，然而效率并不高。。

于是学了下莫比乌斯反演（资料百度找）

求出mo数组后

设f(x)为gcd为x的种数

F(x)为gcd为x倍数的种数

那么显然F(x) = (b / x) * (d / x)

莫比乌斯反演之后，得到f(x) = sum(mo[i] * F(i))。

然后还要容斥减去对称重复的。对称重复的情况为min(b, d)小的中，求一遍除2，（因为存在x = y的情况只有（1,1）一种）

最后还要注意特判下k == 0的情况

代码：

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

const int N = 100005;

int mo[N], prime[N], pn;
bool vis[N];

void Moblus() {
    memset(vis, false, sizeof(vis));
    mo[1] = 1; pn = 0;
    for(int i = 2; i < N; i++) {
        if(!vis[i]) {
            prime[pn++] = i;
            mo[i] = -1;
        }
        for(int j = 0; j < pn; j++) {
            if(i * prime[j] >= N) break;
            vis[i * prime[j]] = true;
            if(i % prime[j] == 0) {
                mo[i * prime[j]] = 0;
                break;
            } else mo[i * prime[j]] = -mo[i];
        }
    }
}

typedef long long ll;
int t, a, b, c, d, k;

int main() {
    Moblus();
    int cas = 0;
    scanf("%d", &t);
    while (t--) {
        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 ans = 0;
        for (int i = 1; i <= b; i++)
            ans += (ll)mo[i] * (b / i) * (d / i);
        ll tmp = 0;
        for (int i = 1; i <= b; i++)
            tmp += (ll)mo[i] * (b / i) * (b / i);
        printf("Case %d: %I64d\n", ++cas, ans - tmp / 2);
    }
    return 0;
}