E. Number Challenge(清晰地推式子)

E. Number Challenge​

问题：${\sum }_{i=1}^a{\sum }_{j=1}^b{\sum }_{k=1}^{c}d(ijk)$，$d(n)$表示n的因子个数。

推式子：
$$\begin{gathered} \sum _{i=1}^a\sum _{j=1}^b\sum _{k=1}^{c}d(ijk) \\ \sum _{i=1}^a\sum _{j=1}^b\sum _{k=1}^c\sum _{x|i}\sum _{y|j}\sum _{z|k}[gcd(x,y)=1][gcd(y,z)=1][gcd(x,z)=1] \\ \sum _{x=1}^a\sum _{y=1}^b\sum _{z=1}^c\sum _{i=1}^a\sum _{j=1}^b\sum _{k=1}^c\sum _{x|i}\sum _{y|j}\sum _{z|k}[gcd(x,y)=1][gcd(y,z)=1][gcd(x,z)=1] \\ \sum _{x=1}^a\sum _{y=1}^b\sum _{z=1}^c\lfloor \frac{a}{x}\rfloor \lfloor \frac{b}{y}\rfloor \lfloor \frac{c}{z}\rfloor [gcd(x,y)=1][gcd(y,z)=1][gcd(x,z)=1] \\ \sum _{x=1}^a\sum _{y=1}^b\sum _{z=1}^c\lfloor \frac{a}{x}\rfloor \lfloor \frac{b}{y}\rfloor \lfloor \frac{c}{z}\rfloor \sum _{d|gcd(x,y)}\mu (d)[gcd(y,z)=1][gcd(x,z)=1] \\ \sum _{d=1}^{min(a,b)}\sum _{x=1}^a\sum _{y=1}^b\sum _{z=1}^c\lfloor \frac{a}{x}\rfloor \lfloor \frac{b}{y}\rfloor \lfloor \frac{c}{z}\rfloor \sum _{d|gcd(x,y)}\mu (d)[gcd(y,z)=1][gcd(x,z)=1] \\ \sum _{d=1}^{min(a,b)}\mu (d)\sum _{x=1,d|x}^a\sum _{y=1,d|y}^b\sum _{z=1}^c\lfloor \frac{a}{x}\rfloor \lfloor \frac{b}{y}\rfloor \lfloor \frac{c}{z}\rfloor [gcd(y,z)=1][gcd(x,z)=1] \\ \sum _{d=1}^{min(a,b)}\mu (d)\sum _{x=1}^{\lfloor \frac{a}{d}\rfloor }\sum _{y=1}^{\lfloor \frac{b}{d}\rfloor }\sum _{z=1}^c\lfloor \frac{a}{dx}\rfloor \lfloor \frac{b}{dy}\rfloor \lfloor \frac{c}{z}\rfloor [gcd(y,z)=1][gcd(x,z)=1] \\ \sum _{d=1}^{min(a,b)}\mu (d)\sum _{z=1}^c\sum _{x=1}^{\lfloor \frac{a}{d}\rfloor }[gcd(x,z)=1]\lfloor \frac{a}{dx}\rfloor \sum _{y=1}^{\lfloor \frac{b}{d}\rfloor }\lfloor \frac{b}{dy}\rfloor \lfloor \frac{c}{z}\rfloor [gcd(y,z)=1] \end{gathered}$$

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 2e3+10;
const double eps = 1e-5;
const int inf = 0x7ffffff;
const ll mod = 1073741824;
int p[N], tot, mu[N];
int a, b, c, g[N][N];
bool st[N];

//求出莫比乌斯函数mu[].
void init() {
    st[1] = mu[1] = 1;
    for(int i=2; i<N; i++) {
        if(!st[i]) p[tot++] = i, mu[i] = -1;
        for(int j=0; j<tot&&i*p[j]<N; j++) {
            st[i*p[j]] = 1;
            if(i % p[j] == 0) break;
            mu[i*p[j]] = -mu[i];
        }
    }
}

//最大公因子
int gcd(int a, int b) {
    if(!b) return a;
    else return gcd(b, a%b);
}

//打表
void get_gcd() {
    for(int i=1; i<=c; i++) {
        for(int j=1; j<=c; j++) {
            g[i][j] = gcd(i, j);
        }
    }
}

//换a,b,c的顺序.
void sort_3() {
    int A = a, B = b, C = c;
    a = min({A, B, C});
    c = max({A, B, C});
    b = A+B+C-a-c; 
}

//\sum_{i=1}^n[gcd(i, m)=1]n/x
ll f(ll n, ll m) {
    ll ans = 0;
    for(int i=1; i<=n; i++) {
        if(g[i][m] == 1) ans += n/i; 
    }
    ans = (ans + mod) % mod;
    return ans;
}


int main() {
#ifndef ONLINE_JUDGE 
    freopen("in.txt", "r", stdin);
    freopen("out.txt", "w", stdout);
#endif
    init();
    cin >> a >> b >> c;
    sort_3();
    get_gcd();
    ll ans = 0;

    //核心。
    for(int i=1; i<=a; i++) {
        for(int j=1; j<=b; j++) {
            if(g[i][j] == 1) {
                ans = (ans + mu[j]*(a/i)%mod * f(b/j, i)%mod * f(c/j, i)%mod + mod) % mod;
            }
        }
    }
    cout << ans << endl;
    return 0;
}