【[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(\text{gcd}(p,q))][\epsilon(\text{gcd}(p,h))][\epsilon(\text{gcd}(q,h))] \]

所以原式=

\[\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^C\sum_{p|i}\sum_{q|j}\sum_{h|k}[\epsilon(\text{gcd}(p,q))][\epsilon(\text{gcd}(p,h))][\epsilon(\text{gcd}(q,h))] \]

交换枚举顺序，枚举因数

\[\sum_{p=1}^A\sum_{q=1}^B\sum_{h=1}^C \lfloor \frac{A}{p} \rfloor \lfloor \frac{B}{q} \rfloor \lfloor \frac{C}{h} \rfloor [\epsilon(\text{gcd}(p,q))][\epsilon(\text{gcd}(p,h))][\epsilon(\text{gcd}(q,h))]\]

把后面那团\(\text{gcd}\)用\(\mu\)函数换掉

\[\sum_{p=1}^A\sum_{q=1}^B\sum_{h=1}^C \sum_{x|\text{gcd}(p,q)}\mu(x) \sum_{y|\text{gcd}(p,h)} \mu(y) \sum_{z|\text{gcd}(q,h)}\mu(z) \]

再换一下枚举顺序

\[\sum_{x=1}^{\text{min(A,B)}}\mu(x) \sum_{y=1}^{\text{min}(A,C)}\mu(y) \sum_{z=1}^{\text{min}(B,C)}\mu(z) \sum_{x|p \ \& \& \ y|p} \lfloor \frac{A}{p} \rfloor \sum_{x|q \ \& \& \ z|q} \lfloor \frac{B}{q} \rfloor \sum_{y | h \ \& \& z|h} \lfloor \frac{C}{h} \rfloor \]

\[\sum_{x=1}^{\text{min(A,B)}}\mu(x) \sum_{y=1}^{\text{min}(A,C)}\mu(y) \sum_{z=1}^{\text{min}(B,C)}\mu(z) \sum_{\text{lcm}(x,y)|p} \lfloor \frac{A}{p} \rfloor \sum_{\text{lcm}(x,z)|q} \lfloor \frac{B}{q} \rfloor \sum_{\text{lcm}(y,z)|h} \lfloor \frac{C}{h} \rfloor \]

然后我们发现这个式子其实不一样的地方只有两部分...
形如\(\sum \limits_{a|b}^c \lfloor \frac{c}{b} \rfloor\)的那三项可以\(\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} \rfloor\)中的\(a\)大于\(c\)，则这个式子的贡献一定为\(0\)，也就是说当\(\text{lcm}(x,y) \geq A\)或\(\text{lcm}(x,z) \geq B\)或\(\text{lcm}(y,z) \geq C\)时，这个式子的贡献一定为\(0\)。
于是我们就可以枚举最大公约数\(g\)，使得\(u=i \times g,v = j \times g, \mu(u) \neq 0,\mu(v) \neq 0\)并保证\(i\)与\(j\)互质且\(\text{lcm}(u,v)=i \times j \times g \leq \text{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\n", ans % mod);
}

int main() {
  sieve(100000);
  scanf("%d", &T);
  while (T--) init(), solve();
  return 0;
}