Codeforces 235 E Number Challenge(莫比乌斯反演)

题目链接：
Codeforces 235 E Number Challenge
题意:

$$记d(i)表示i的约数个数，计算：\sum_{i=1}^{a}\sum_{j=1}^{b}\sum_{k=1}^{c}d(ijk)，a,b,c\in [1,2000]$$
分析：
$$Ans\ = \ \sum_{i=1}^{a}\sum_{j=1}^{b}\sum_{k=1}^{c}d(ijk)\ = \ \sum_{i=1}^{a}\lfloor{\frac{a}{i}}\rfloor\sum_{j=1}^{b}\lfloor{\frac{b}{j}}\rfloor\sum_{k=1}^{c}\lfloor{\frac{c}{k}}\rfloor_{gcd(i,j)= gcd(i,k)= gcd(j,k)=1}$$
对于后面的部分反演可得： $$Ans =\sum_{i=1}^{i=a}{\lfloor{\frac{a}{i}}}\rfloor\sum_{d}^{min(b,c)}μ(d)\sum_{d|j,(i,j)=1}\lfloor{\frac{b}{j}}\rfloor\sum_{d|k,(i,k)=1}\lfloor{\frac{c}{k}}\rfloor$$
$记j=dj′,k=dk′,则：记j=dj′,k=dk′,则：$
$$Ans=\sum_{i=1}^{i=a}{\lfloor{\frac{a}{i}}}\rfloor\sum_{d}^{min(b,c)}μ(d)\sum_{gcd(i,dj')=1}\lfloor{\frac{b}{dj'}}\rfloor\sum_{gcd(i,dk')=1}\lfloor{\frac{c}{dk'}}\rfloor$$ ,
因为 $gcd(i,dj')=1$，我们先保证 $gcd(i,d)=1,$然后枚举 $j':1 \rightarrow \frac{b}{d}$,保证 $gcd(i,j')=1$，这样就可以使得 $gcd(i,dj')=1$，累加即可。对于 $gcd(i,dk')=1$同样处理。
枚举 $i$和$d$的时间复杂度是是 $a$和$min(b,c)$,但是枚举 $j'和k'$的时间复杂度是 $\sum_{i=1}^{i=b}{b/i}$和 $\sum_{i=1}^{i=c}{c/i}$,所以总的时间复杂度大概是 $O(a*b*log(b))$

#include <cstdio>
#include <cstring>
#include <string>
#include <cmath>
#include <algorithm>
#include <bitset>
using namespace std;
typedef long long ll;
const int MAX_N = 2010;
const ll mod = (ll)1 << 30;

int prime_cnt, prime[MAX_N], mu[MAX_N], gcd[MAX_N][MAX_N];
bitset<MAX_N> bs;

void GetMu()
{
    prime_cnt = 0;
    bs.set();
    mu[1] = 1;
    for(int i = 2; i < MAX_N; ++i) {
        if(bs[i]) {
            prime[prime_cnt++] = i;
            mu[i] = -1;
        }
        for(int j = 0; j < prime_cnt && i * prime[j] < MAX_N; ++j) {
            bs[i * prime[j]] = 0;
            if(i % prime[j]) {
                mu[ i * prime[j]] = - mu[i];
            }else {
                mu[i * prime[j]] = 0;
                break;
            }
        }   
    }
}

inline int GCD(int a, int b)
{
    if(gcd[a][b]) return gcd[a][b];
    else return b == 0 ? a : gcd[a][b] = GCD(b, a % b);
}

inline void GetGcd()
{
    for(int i = 1; i < MAX_N; i++){
        for(int j = i; j < MAX_N; j++) {
            gcd[i][j] = gcd[j][i] = GCD(i, j);
        }
    }
}

inline ll work(int n, int x)
{
    ll res = 0;
    for(int i = 1; i <= n; ++i) {
        if(gcd[i][x] == 1) {
            res = (res + (n / i)) % mod;
        }
    }
    return res;
}

inline ll solve(int a, int b, int c)
{
    ll res = 0;
    int top = min(b, c);
    for(int i = 1; i <= a; ++i) {
        for(int d = 1; d <= top; ++d){
            if(gcd[i][d] == 1) {
                ll tmp = 0;
                tmp = (ll) (a / i) * mu[d] * work(b / d, i) % mod * work(c / d, i) % mod;
                res = ((res + tmp) % mod + mod) % mod;
            }
        }
    }
    return res;
}

int main()
{
    GetMu();
    GetGcd();
    int a, b, c;
    while(~scanf("%d%d%d", &a, &b, &c)) {
        printf("%lld\n", solve(a, b, c));
    }
    return 0;
}