【[SDOI2018]旧试题】

一道很有趣的题目。
题目大意：求

[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;
}