题目链接：LOJ 2565

时光匆匆，转眼间又是一年省选季……

这是小 Q 同学第二次参加省队选拔赛。今年，小 Q 痛定思痛，不再冒险偷取试题，而是通过练习旧试题提升个人实力。可是旧试题太多了，小 Q 没日没夜地做题，却看不到前方的光明在哪里。

一天，因做题过度而疲惫入睡的小 Q 梦到自己在考场上遇到了一道好像做过的题目，却怎么也想不起曾经自己是怎么解决它的，直到醒来还心有余悸。

小 Q 眉头一皱，感觉事情不妙，于是他找到了你，希望你能教他解决这道题目。小 Q 依稀记得题目要计算如下表达式的值

$$ \left(\sum_{i = 1}^{A} \sum_{j = 1}^{B} \sum_{k = 1}^{C} d(i j k) \right) \bmod (10^9 + 7) $$

其中 $d(i j k)$ 表示 $i\times j\times k$ 的约数个数。

数据范围：$1 \le T \le 10, 1 \le A, B, C \le 10^5, 1 \le \sum \max(A, B, C) \le 2 \cdot 10^5$。

Solution

对于 $d(i j k)$，有如下等式：

$$ d(i j k) = \sum_{u \mid i} \sum_{v \mid j} \sum_{w \mid k} [(u, v) = 1] [(v, w) = 1] [(w, u) = 1] $$

推式子：

$$ \sum_{i = 1}^{A} \sum_{j = 1}^{B} \sum_{k = 1}^{C} \sum_{u \mid i} \sum_{v \mid j} \sum_{w \mid k} [(u, v) = 1] [(v, w) = 1] [(w, u) = 1] \\ \sum_{i = 1}^{A} \sum_{j = 1}^{B} \sum_{k = 1}^{C} \sum_{u \mid i} \sum_{v \mid j} \sum_{w \mid k} \varepsilon((u, v)) \varepsilon((v, w)) \varepsilon((w, u)) \\ \sum_{u = 1}^{A} \sum_{v = 1}^{B} \sum_{w = 1}^{C} \left\lfloor\frac{A}{u}\right\rfloor \left\lfloor\frac{B}{v}\right\rfloor \left\lfloor\frac{C}{w}\right\rfloor \varepsilon((u, v)) \varepsilon((v, w)) \varepsilon((w, u)) \\ \sum_{u = 1}^{A} \sum_{v = 1}^{B} \sum_{w = 1}^{C} \left\lfloor\frac{A}{u}\right\rfloor \left\lfloor\frac{B}{v}\right\rfloor \left\lfloor\frac{C}{w}\right\rfloor \sum_{i \mid (u, v)} \mu(i) \sum_{j \mid (v, w)} \mu(j) \sum_{k \mid (w, u)} \mu(k) \\ \sum_{i = 1}^{A} \sum_{j = 1}^{B} \sum_{k = 1}^{C} \mu(i) \mu(j) \mu(k) \sum_{\operatorname{lcm}(k, i) \mid u} \left\lfloor\frac{A}{u}\right\rfloor \sum_{\operatorname{lcm}(i, j) \mid v} \left\lfloor\frac{B}{v}\right\rfloor \sum_{\operatorname{lcm}(j, k) \mid w} \left\lfloor\frac{C}{w}\right\rfloor $$

将所有的 $\sum$ 的枚举上界都变为 $n = \max(A, B, C)$，这对式子的值没有影响。

设 $f_A(x) = \sum_{x \mid i} \left\lfloor\frac{A}{i}\right\rfloor$，原式化为：

$$ \sum_{i = 1}^{n} \sum_{j = 1}^{n} \sum_{k = 1}^{n} \mu(i) \mu(j) \mu(k) f_A(\operatorname{lcm}(k, i)) f_B(\operatorname{lcm}(i, j)) f_C(\operatorname{lcm}(j, k)) $$

直接枚举有序三元组 $(i, j, k)$ 对复杂度没有优化，但是可以发现很多三元组 $(i, j, k)$ 是无用的。

• 含有任意一个数使得 $\mu(i) = 0$ 的三元组。
• 含有任意两个数使得 $\operatorname{lcm}(i, j) > n$ 的三元组。

如果把 $\mu(i) = 0$ 的数都删掉，在 $\operatorname{lcm}(i, j) \le n$ 的数之间连边，问题就几乎转化为：对于这张图中的所有三元环，计算对答案的贡献。（「几乎」的意思是，还需要额外计算 $(i, i, i)$ 和 $(i, i, j)$ 等类型的三元组。）

连边时，先枚举 $u, v$ 的 $\gcd(u, v) = d$，再枚举 $i = u / d, j = v / d$ 且满足 $\gcd(i, j) = 1$，则 $\operatorname{lcm}(u, v) = ijd \le n$，可以发现边的数量最多只有 $760741$ 条。

对于找三元环的过程，直接采用「三元环计数」的方法即可。

时间复杂度：$O(m \sqrt m)$，其中 $m = 760741$。

Code

#include <bits/stdc++.h>
typedef long long int64;

const int N = 1e5 + 5, M = 760741 + 5, MOD = 1e9 + 7;
int A, B, C, n, tot, cnt, prime[N], d[N], mu[N], p[N], deg[N], vis[N], fa[N], fb[N], fc[N];
std::vector<std::pair<int, int>> E[N];
struct Edge {
    int u, v, l;
} e[M];

void init(int n) {
    std::bitset<N> flg;
    flg.set();
    mu[1] = 1;
    p[cnt = 1] = 1;
    for (int i = 2; i <= n; i++) {
        if (flg[i]) {
            prime[++tot] = i;
            mu[i] = -1;
        }
        for (int j = 1; j <= tot && i * prime[j] <= n; j++) {
            flg[i * prime[j]] = 0;
            if (i % prime[j] == 0) {
                mu[i * prime[j]] = 0;
                break;
            }
            mu[i * prime[j]] = -mu[i];
        }
        if (mu[i]) {
            p[++cnt] = i;
        }
    }
    p[cnt + 1] = INT_MAX;
}
void calc(int *f, int m) {
    std::fill(f + 1, f + n + 1, 0);
    for (int i = 1; i <= m; i++) {
        for (int j = i; j <= m; j += i) {
            f[i] += m / j;
        }
    }
}
void solve() {
    scanf("%d%d%d", &A, &B, &C);
    n = std::max({A, B, C});
    calc(fa, A);
    calc(fb, B);
    calc(fc, C);
    int cntE = 0;
    std::fill(deg + 1, deg + n + 1, 0);
    for (int i = 1; p[i] <= n; i++) {
        E[p[i]].clear();
    }
    for (int g = 1; p[g] <= n; g++) {
        for (int i = 1; p[i] * p[g] <= n; i++) {
            if (mu[p[i] * p[g]]) {
                int lim = n / p[i] / p[g];
                for (int j = i + 1; p[j] <= lim; j++) {
                    if (std::__gcd(p[i], p[j]) == 1 && mu[p[j] * p[g]]) {
                        int u = p[i] * p[g], v = p[j] * p[g], l = p[i] * p[j] * p[g];
                        e[++cntE] = Edge{u, v, l};
                        deg[u]++;
                        deg[v]++;
                    }
                }
            }
        }
    }
    for (int i = 1; i <= cntE; i++) {
        int u = e[i].u, v = e[i].v;
        if (deg[u] < deg[v] || (deg[u] == deg[v] && u > v)) {
            std::swap(u, v);
        }
        E[u].emplace_back(v, e[i].l);
    }
    int64 ans = 0;
    for (int i = 1; p[i] <= n; i++) {
        int64 cur = (int64)fa[p[i]] * fb[p[i]] * fc[p[i]];
        mu[p[i]] > 0 ? ans += cur : ans -= cur;
    }
    for (int i = 1; i <= cntE; i++) {
        int u = e[i].u, v = e[i].v, l = e[i].l;
        int64 cur = (int64)fa[u] * fb[l] * fc[l] + (int64)fa[l] * fb[u] * fc[l] + (int64)fa[l] * fb[l] * fc[u];
        mu[v] > 0 ? ans += cur : ans -= cur;
        cur = (int64)fa[v] * fb[l] * fc[l] + (int64)fa[l] * fb[v] * fc[l] + (int64)fa[l] * fb[l] * fc[v];
        mu[u] > 0 ? ans += cur : ans -= cur;
    }
    for (int i = 1; p[i] <= n; i++) {
        int u = p[i];
        for (auto e : E[u]) {
            vis[e.first] = e.second;
        }
        for (auto e1 : E[u]) {
            int v = e1.first;
            for (auto e2 : E[v]) {
                int w = e2.first;
                if (vis[w]) {
                    int x = vis[w], y = e1.second, z = e2.second;
                    int64 cur = (int64)fa[x] * ((int64)fb[y] * fc[z] + (int64)fb[z] * fc[y])
                          + (int64)fb[x] * ((int64)fc[y] * fa[z] + (int64)fc[z] * fa[y])
                          + (int64)fc[x] * ((int64)fa[y] * fb[z] + (int64)fa[z] * fb[y]);
                    mu[u] * mu[v] * mu[w] > 0 ? ans += cur : ans -= cur;
                }
            }
        }
        for (auto e : E[u]) {
            vis[e.first] = 0;
        }
    }
    printf("%lld\n", ans % MOD);
}
int main() {
    init(N - 5);
    int T;
    scanf("%d", &T);
    for (int cs = 1; cs <= T; cs++) {
        solve();
    }
    return 0;
}