一道很有趣的题目。
题目大意：求
[sum_{i=1}^A sum_{j=1}^B sum_{k=1}^C sigma_0(ijk)
]
因为
[sigma_0(ijk)=sum_{p|i}sum_{q|j}sum_{h|k}[epsilon( ext{gcd}(p,q))][epsilon( ext{gcd}(p,h))][epsilon( ext{gcd}(q,h))]
]
所以原式=
[sum_{i=1}^Asum_{j=1}^Bsum_{k=1}^Csum_{p|i}sum_{q|j}sum_{h|k}[epsilon( ext{gcd}(p,q))][epsilon( ext{gcd}(p,h))][epsilon( ext{gcd}(q,h))]
]
交换枚举顺序，枚举因数
[sum_{p=1}^Asum_{q=1}^Bsum_{h=1}^C
lfloor frac{A}{p}
floor lfloor frac{B}{q}
floor lfloor frac{C}{h}
floor [epsilon( ext{gcd}(p,q))][epsilon( ext{gcd}(p,h))][epsilon( ext{gcd}(q,h))]]
把后面那团( ext{gcd})用(mu)函数换掉
[sum_{p=1}^Asum_{q=1}^Bsum_{h=1}^C sum_{x| ext{gcd}(p,q)}mu(x) sum_{y| ext{gcd}(p,h)} mu(y) sum_{z| ext{gcd}(q,h)}mu(z)
]
再换一下枚举顺序
[sum_{x=1}^{ ext{min(A,B)}}mu(x) sum_{y=1}^{ ext{min}(A,C)}mu(y) sum_{z=1}^{ ext{min}(B,C)}mu(z) sum_{x|p & & y|p} lfloor frac{A}{p}
floor sum_{x|q & & z|q} lfloor frac{B}{q}
floor sum_{y | h & & z|h} lfloor frac{C}{h}
floor
]
[sum_{x=1}^{ ext{min(A,B)}}mu(x) sum_{y=1}^{ ext{min}(A,C)}mu(y) sum_{z=1}^{ ext{min}(B,C)}mu(z) sum_{ ext{lcm}(x,y)|p} lfloor frac{A}{p}
floor sum_{ ext{lcm}(x,z)|q} lfloor frac{B}{q}
floor sum_{ ext{lcm}(y,z)|h} lfloor frac{C}{h}
floor
]
然后我们发现这个式子其实不一样的地方只有两部分...
形如(sum limits_{a|b}^c lfloor frac{c}{b}
floor)的那三项可以(mathcal{O}(n log n))预处理出来。
接下来才是这道题的精髓。
在(mathcal{O}(n^3))枚举三元组((x,y,z))的同时，我们可以发现有很多时候他们的情况是(0)。
1：(mu(x)=0)或(mu(y)=0)或(mu(z)=0)。
2：若(sum limits_{a|b}^c lfloor frac{c}{b}
floor)中的(a)大于(c)，则这个式子的贡献一定为(0)，也就是说当( ext{lcm}(x,y) geq A)或( ext{lcm}(x,z) geq B)或( ext{lcm}(y,z) geq C)时，这个式子的贡献一定为(0)。
于是我们就可以枚举最大公约数(g)，使得(u=i imes g,v = j imes g, mu(u)
eq 0,mu(v)
eq 0)并保证(i)与(j)互质且( ext{lcm}(u,v)=i imes j imes g leq ext{max}(A,B,C))，连边((u,v))，最后统计三元环的时候顺便计数就好了。
因为(mathsf{s{color{red}hadowice1984}})大佬告诉我们边数最多只有(760741)条，所以(m sqrt m)是可过的。
代码：
#include <bits/stdc++.h>
typedef long long ll;
const int N = 1e5 + 200;
const ll mod = 1e9 + 7;
typedef std::pair<int, int> pii;
int n, m, i, j, k, prim_tot, T, A, B, C, mx, mn, cnte;
ll _1[N], _2[N], _3[N];
int vis[N], deg[N];
int mu[N], prim[N];
bool _vis[N];
ll ans;
inline void sieve(int n) {
mu[1] = 1;
for (int i = 2; i <= n; i++) {
if (!_vis[i]) prim[++prim_tot] = i, mu[i] = -1;
for (int j = 1; j <= prim_tot && i * prim[j] <= n; j++) {
_vis[i * prim[j]] = 1;
if (i % prim[j] == 0) break;
else mu[i * prim[j]] = -mu[i];
}
}
}
struct edge {
int u, v, w;
edge() { u = v = w = 0; }
edge(int _u, int _v, int _w) {
u = _u, v = _v, w = _w;
}
} edges[N << 4];
std::vector<pii> G[N];
inline void init() {
scanf("%d %d %d", &A, &B, &C), ans = 0, cnte = 0;
mx = std::max(std::max(A, B), C);
mn = std::min(std::min(A, B), C);
for (int i = 1; i <= mx; i += 8) {
_1[i] = _1[i + 1] = _1[i + 2] = _1[i + 3] = _1[i + 4] = _1[i + 5] = _1[i + 6] = _1[i + 7] = 0;
_2[i] = _2[i + 1] = _2[i + 2] = _2[i + 3] = _2[i + 4] = _2[i + 5] = _2[i + 6] = _2[i + 7] = 0;
_3[i] = _3[i + 1] = _3[i + 2] = _3[i + 3] = _3[i + 4] = _3[i + 5] = _3[i + 6] = _3[i + 7] = 0;
deg[i] = deg[i + 1] = deg[i + 2] = deg[i + 3] = deg[i + 4] = deg[i + 5] = deg[i + 6] = deg[i + 7] = 0;
}
for (int i = 1; i <= mx; i++) G[i].clear();
for (int i = 1; i <= A; i++) {
for (int j = i; j <= A; j += i)
_1[i] += A / j;
}
for (int i = 1; i <= B; i++) {
for (int j = i; j <= B; j += i)
_2[i] += B / j;
}
for (int i = 1; i <= C; i++) {
for (int j = i; j <= C; j += i)
_3[i] += C / j;
}
}
inline void solve() {
for (int i = 1; i <= mn; i++) {
if (mu[i]) ans += 1ll * mu[i] * mu[i] * mu[i] * _1[i] * _2[i] * _3[i];
}
for (int g = 1; g <= mx; g++) {
for (int i = 1; i * g <= mx; i++) {
if (!mu[i * g]) continue;
for (int j = i + 1; 1ll * i * j * g <= mx; j++) {
if (mu[j * g] && std::__gcd(i, j) == 1) {
int u = i * g, v = j * g, lcm = i * j * g;
int p = mu[u] * mu[u] * mu[v], q = mu[u] * mu[v] * mu[v];
ans += p * (_1[u] * _2[lcm] * _3[lcm] + _1[lcm] * _2[u] * _3[lcm] + _1[lcm] * _2[lcm] * _3[u]);
ans += q * (_1[v] * _2[lcm] * _3[lcm] + _1[lcm] * _2[v] * _3[lcm] + _1[lcm] * _2[lcm] * _3[v]);
deg[u]++, deg[v]++;
edges[++cnte] = edge(u, v, lcm);
}
}
}
}
for (int i = 1; i <= cnte; i++) {
int u = edges[i].u, v = edges[i].v, w = edges[i].w;
if (deg[u] < deg[v] || (deg[u] == deg[v] && u > v)) std::swap(u, v);
G[u].push_back(std::make_pair(v, w));
}
for (int k = 1; k <= mx; k++) {
if (!mu[k]) continue;
int SZ = G[k].size();
for (int i = 0; i < SZ; i++) vis[G[k][i].first] = G[k][i].second;
for (int i = 0; i < SZ; i++) {
int v = G[k][i].first, w = G[k][i].second;
if (!mu[v]) continue;
for (int j = 0, _SZ = G[v].size(); j < _SZ; j++) {
int vv = G[v][j].first;
if (!vis[vv]) continue;
int wyz = G[v][j].second, wxz = vis[vv], mmu = mu[k] * mu[v] * mu[vv];
if (!mmu) continue;
ans += mmu * (_1[w] * _2[wyz] * _3[wxz] + _1[w] * _2[wxz] * _3[wyz]
+ _1[wyz] * _2[w] * _3[wxz] + _1[wyz] * _2[wxz] * _3[w]
+ _1[wxz] * _2[w] * _3[wyz] + _1[wxz] * _2[wyz] * _3[w]);
}
}
for (int i = 0; i < SZ; i++) vis[G[k][i].first] = 0;
}
printf("%lld
", ans % mod);
}
int main() {
sieve(100000);
scanf("%d", &T);
while (T--) init(), solve();
return 0;
}