【HDU1695】GCD（莫比乌斯反演+优化）

这道题已经拿容斥做过了，题意可以看我的另一篇博客http://blog.csdn.net/brazy/article/details/77256645

最近正好在学习莫比乌斯反演，这题就再拿莫比乌斯反演来做一遍：

设f(n)表示gcd(x,y)=n的数量，F(n)表示gcd(x,y)是n的倍数的数量。

而F(n) = (x/n) * (y/n) 要加括号，因为这是取整之后的乘积

根据定义用第二种形式：f(n) = sigma(mu[d/n]*F(d)) (n|d)

这样只要枚举k的倍数一直到min(n,m)就可以了。

可是如果k=1，那么枚举一次就是O(N)；总复杂度为O(N*N);

实际上可以继续优化；solve(n,m,k)等价于solve(n/k,m/k)表示x在[1,n/k],y在[1,m/k]，gcd(x,y)==1的对数。

由于x/i,x/(i+1),x/(i+2)...x/(i+t)存在连续相同的结果，也就是这段区间[l,r]内(n/i)*(m/i)的结果是相同的；

这样i在[l,r] 范围内的(n/i)*(m/i)*mu[i];就等价于 (n/i)*(m/i)*(sum[r]-sum[l-1])； 

sum表示mu的前缀和。所以这里可以快速处理。复杂度为sqrt(N); 总时间复杂度为N*sqrt(N);

参考：https://wenku.baidu.com/view/fbec9c63ba1aa8114431d9ac.html

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <vector>
#include <cmath>
#include <algorithm>
#include <string>
#include <cstring>
#include <map>

using namespace std;

const int MAXN = 1000000;
bool check[MAXN + 10];
int prime[MAXN + 10];
int mu[MAXN + 10];
void Moblus()
{
    memset(check, false, sizeof(check));
    mu[1] = 1;
    int tot = 0;
    for(int i = 2; i <= MAXN; ++i)
    {
        if(!check[i])
        {
            prime[tot++] = i;
            mu[i] = -1;
        }
        for(int j = 0; j < tot; ++j)
        {
            if(i * prime[j] > MAXN) break;
            check[i * prime[j]] = true;
            if(i % prime[j] == 0)
            {
                mu[i * prime[j]] = 0;
                break;
            }
            else
            {
                mu[i * prime[j]] = -mu[i];
            }
        }
        mu[i] += mu[i - 1];
    }
}
long long solve(int m, int n)
{
    long long ans = 0LL;
    int last = 0;
    for(int i = 1; i <= min(m, n); i = last + 1)
    {
        last = min(n / (n/i), m / (m/i));
        ans = ans + 1LL * (mu[last] - mu[i - 1]) * (m / i) * (n / i);
    }
    return ans;
}

int main()
{
    Moblus();
    int CASE;
    scanf("%d", &CASE);
    for(int cas = 1; cas <= CASE; ++cas)
    {
        int a, b, c, d, k;
        scanf("%d%d%d%d%d", &a, &b, &c, &d, &k);
        if(k == 0 || b < k || d < k)
        {
            printf("Case %d: 0\n", cas);
            continue;
        }
        if(b > d) swap(b, d);
        b /= k;
        d /= k;
        printf("Case %d: %lld\n", cas, solve(b, d) - solve(b, b) / 2);
    }
    return 0;
}