upd: （与 zky 讨论）其实预处理部分可以 fft 然后修误差，这样就能做到 $\tilde O(N)$ 了。

uupd：其实是 zky 的 subset。

• 参考资料：《数论选讲》in APIO2024 by zhoukangyang.
• Notations：$N$ 为 $A,B,C$ 的上界；记前缀和 $S_F(n)=\sum_{i=1}^nF(i)$。

讲义里提到本题的 $O(N\sqrt{N})$ 做法，但没有细说。我将首先详细地陈述此做法，然后提出一个 $O(N^{5/4}+N\log^3N)$ 的做法。

• 多元积性函数

  函数 $f(a,b,c)$ 满足当 $abc\perp xyz$ 时，有 $f(ax,by,cz)=f(a,b,c)f(x,y,z)$，则称 $f$ 为（三元）积性函数。

• Observation：令 $f(i,j,k)=d(ijk)$，则 $f$ 是积性函数。问题转化为三元积性函数的矩阵和。

• 积性分解

  将 $x,y,z$ 分解为 $x=\prod_{i=1}^mp_i^{a_i}$，$y=\prod_{i=1}^mp_i^{b_i}$，$z=\prod_{i=1}^mp_i^{c_i}$，对于积性函数 $f$

  $$f(x,y,z)=\prod_{i=1}^mf(p^{a_i},p^{b_i},p^{c_i})$$

• 狄利克雷卷积

  三元积性函数的狄利克雷卷积定义为

  $$(f*g)(x,y,z)=\sum_{a|x}\sum_{b|y}\sum_{c|z}f(a,b,c)g(x/a,y/b,z/c)$$

  类似一元的情况，有如下性质成立：

  • 积性函数的卷积是积性函数。
  • 积性函数的逆是积性函数。

• 贝尔级数

  在素数 $p$ 处观测积性函数 $f(x,y,z)$，定义贝尔级数为三元形式幂级数

  $$\mathcal F_p(u_1,u_2,u_3)=\sum_{i=0}^{+\infty}\sum_{j=0}^{+\infty}\sum_{k=0}^{+\infty}f(p^i,p^j,p^k)u_1^iu_2^ju_3^k$$

  类似一元的情况，两个函数卷积，对应的贝尔级数相乘。

• 低阶拟合法

  模仿 Powerful Number 的思路，构造一个数论函数 $g$ 使得 $f,g$ 在素数幂次较低时相等，然后做除法 $h=f/g$，$h$ 的非零元素有希望较少。

  计算 $f$ 的矩阵和时

  $$\begin{aligned} S_f(A,B,C)&=\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^Cf(i,j,k)\\ &=\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^C\sum_{x|i}\sum_{y|j}\sum_{z|j}h(x,y,z)g(i/x,j/y,k/z)\\ &=\sum_{x=1}^A\sum_{y=1}^B\sum_{z=1}^Ch(x,y,z)\sum_{i=1}^{\lfloor A/x\rfloor}\sum_{j=1}^{\lfloor B/y\rfloor}\sum_{j=1}^{\lfloor C/z\rfloor}g(i,j)\\ &=\sum_{x=1}^A\sum_{y=1}^B\sum_{z=1}^Ch(x,y,z)S_g(\lfloor A/x\rfloor,\lfloor B/y\rfloor,\lfloor C/z\rfloor)\\ \end{aligned}$$

  如果 $h$ 稀疏，且 $g$ 的矩阵和易求，我们就得到了一个不错的算法。

  一个自然的想法是，构造 $g(i,j,k)=d(i)d(j)d(k)$，则

  • $g$ 容易计算。$$\begin{aligned} S_g(A,B,C)&=\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^Cd(i)d(j)d(k)\\ &=\left(\sum_{i=1}^Ad(i)\right)\left(\sum_{j=1}^Bd(j)\right)\left(\sum_{k=1}^Cd(k)\right)\\ &=S_d(A)S_d(B)S_d(C) \end{aligned}$$

  ​ $O(N)$ 预处理 $S_d$ 后容易 $O(1)$ 计算。

  • $g$ 和 $f$ 在低阶处相等。

    可以发现 $g(p^k,1,1)=f(p^k,1,1)$，$g(1,p^k,1)=f(1,p^k,1)$，$g(1,1,p^k)=f(1,1,p^k)$。

    为了更直观地说明这一点，构造贝尔级数。

    $$\begin{matrix} \begin{aligned} f(p^i,p^j,p^k)&=d(p^{i+j+k})=i+j+k+1\\ \mathcal F_p(u_1,u_2,u_3)&=\sum_{i=0}^{+\infty}\sum_{j=0}^{+\infty}\sum_{k=0}^{+\infty}f(p^i,p^j,p^k)u_1^iu_2^ju_3^k\\ &=\sum_{i=0}^{+\infty}\sum_{j=0}^{+\infty}\sum_{k=0}^{+\infty}(i+j+k+1)u_1^iu_2^ju_3^k \end{aligned}\\\\ \begin{aligned} g(p^i,p^j,p^k)&=d(p^i)d(p^j)d(p^k)=(i+1)(j+1)(k+1)\\ \mathcal G_p(u_1,u_2,u_3)&=\sum_{i=0}^{+\infty}\sum_{j=0}^{+\infty}\sum_{k=0}^{+\infty}g(p^i,p^j,p^k)u_1^iu_2^ju_3^k\\ &=\sum_{i=0}^{+\infty}\sum_{j=0}^{+\infty}\sum_{k=0}^{+\infty}(i+1)(j+1)(k+1)u_1^iu_2^ju_3^k \end{aligned}\\\\ \mathcal H_p=\mathcal F_p/\mathcal G_p \end{matrix}$$

    $\mathcal F$ 和 $\mathcal G$ 的系数在“三条棱”上是相等的，根据幂级数逆元的知识，可以发现 $\mathcal H$ 在这三条棱上都是 $1$，没有更多的项。即是说：

    $$h(1,1,p^k)=h(1,p^k,1)=h(p^k,1,1)=0$$

    不仅如此，计算发现：

    $$\begin{matrix} h(1,1,p)=h(1,p,1)=h(p,1,1)=-1\\ h(1,1,p^2)=h(1,p^2,1)=h(p^2,1,1)=2 \end{matrix}$$

    除了上述项，其余都是零。

• $h$ 密度的分析

  枚举量就是 $h(1\sim A,1\sim B,1\sim C)$ 的非零项数量，将条件放宽为 $h(x,y,z)$，$xyz\leq ABC\leq N^3$。

  记 $s(n)$ 为满足 $xyz=n$ 的非零 $h(x,y,z)$ 个数，则复杂度为 $O(S_s(N^3))$。

  可以证明 $s$ 是积性函数，且 $s(p)=0$，$s(p^2)=3$。

  • Lemma: 对于数论函数 $F$，找出最小的 $k_0\geq 1$ 满足 $F(p^{k_0})\neq 0$。记 $c$ 为 $F(p^{k_0})$ 的最大值，其余 $F(p^k)=O({\rm poly}(k))$，则 $S_F(n)=O(n^{1/k_0}\log^{c-1}n)$。

  根据这个引理，可以得到 $S_s(n)=O(n^{1/2}\log^2n)$，对应复杂度 $O(N\sqrt{N}\log^2 N)$。

  此界太过宽松，这是因为我们将矩形区域粗暴地放缩成了双曲区域 $xyz\leq ABC$，两者的面积比就是 $\Theta(\log^2N)$。为了去除多余的 log 因子，必须考虑为矩形区域。

  定义三元数论函数 $v(x,y,z)=[h(x,y,z)\neq 0]$，我们的目标是估计 $v$ 的矩阵和。

  令 $v_1$ 满足 $v_1(p,p,1)=v_1(p,1,p)=v_1(1,p,p)=1$，其余都是 $0$，$v_2=v-v_1$。显然，$v\leq v_1*v_2$。

  先考察 $v_1$ 的矩阵和，$v_1(a,b,c)$ 总是可以写作 $v_1(xy,xz,yz)$，枚举 $x,y,z$ 进行计数：

  $$\begin{aligned} S_{v_1}(A,B,C)&=\sum_{x=1}^{+\infty}\sum_{y=1}^{+\infty}\sum_{z=1}^{+\infty}[x\perp y][x\perp z][y\perp z][xy\leq A][xz\leq B][yz\leq C]\\ &\leq \sum_{x=1}^{+\infty}\sum_{y=1}^{+\infty}\sum_{z=1}^{+\infty}[xy\leq A][xz\leq B][yz\leq C]\\ &\sim\sum_{x=1}^{+\infty}\sum_{y=1}^{+\infty}[xy\leq A]\min(B/x,C/y)\\ \end{aligned}$$

  用讨论将 $\min$ 拆掉，$C/y$ 较小的部分是

  $$\begin{aligned} &\sim C\sum_{x=1}^{\sqrt{AB/C}}\sum_{Cx/B\leq y\leq A/x}^{+\infty}1/y\\ &\sim C\sum_{x=1}^{\sqrt{AB/C}}\ln(AB/Cx^2)\\ &\sim C\cdot T\Big(\sqrt{AB/C}\Big) \end{aligned}$$

  根据对称性，另一部分是 $B\cdot T\Big(\sqrt{AC/B}\Big)$。

  其中 $T(n)=\sum_{i=1}^n\ln(n/i)$，这是一个经典形式

  $$\begin{aligned} T(2n)&=\sum_{i=1}^{2n}\ln(2n/i)\\ &=\sum_{i=1}^{n}\big(\ln 2+\ln(n/i)\big)+\sum_{i=n+1}^{2n}\ln(2n/i)\\ &=(\ln2)\cdot 2n+T(n)\\ &=O(n)+T(n) \end{aligned}$$

  解得 $T(n)=O(n)$。代入得 $S_{v_1}(A,B,C)=O(\sqrt{ABC})$。

  考虑 $v_1*v_2$ 的矩阵和

  $$\begin{aligned} S_{v}(N,N,N)&=\sum_{x,y,z\leq N}v_2(x,y,z)S_{v_1}(N/x,N/y,N/z)\\ &=\sum_{x,y,z\leq N}v_2(x,y,z)\sqrt{\dfrac{N^3}{xyz}} \end{aligned}$$

  记 $s_2(n)$ 为满足 $xyz=n$ 的非零 $v_2(x,y,z)$ 个数。$s_2$ 的全零前缀更长，根据引理，$v_2(x,y,z)$ 在 $xyz\leq n$ 的非零个数不超过 $S_{s_2}(n)=\tilde O(n^{1/3})$，密度可以估计为 $\tilde O(n^{-2/3})$。枚举 $T=xyz$，复杂度可以估计为

  $$\begin{aligned} \sum_{T\leq N^3}\sqrt{\dfrac{N^3}{T}}s_2(T)&=\sum_{T\leq N^3}\sqrt{\dfrac{N^3}{T}}\tilde O(T^{-2/3})\\ &= N^{3/2}\sum_{T\leq N^3}\tilde O(T^{-7/6})\\ &= O(N^{3/2}) \end{aligned}$$

  综上，我们证明了枚举量是 $O(N\sqrt{N})$。

代码如下：

#include<cstdio>
#define MaxN 100050
using namespace std;
const int mod=1000000007;
int p[MaxN>>3],tn,d[MaxN];
void sieve(int n)
{
  d[1]=1;
  for (int i=2;i<=n;i++) {
    if (!d[i]) d[p[++tn]=i]=2;
    for (int j=1,t;j<=tn&&(t=i*p[j])<=n;j++) {
      if (i%p[j]==0) {d[t]=2*d[i]-d[i/p[j]];break;}
      d[t]=2*d[i];
    }
  }
  for (int i=1;i<=n;i++) d[i]+=d[i-1];
}
int A,B,C,ret;
void dfs(int x, int y, int z, int k, int now) {
  ret=(ret+1ll*now*d[A/x]*d[B/y]%mod*d[C/z])%mod;
  for (int j=k;j<=tn;j++) {
    int flag=(1ll*x*p[j]<=A)|((1ll*y*p[j]<=B)<<1)|((1ll*z*p[j]<=C)<<2);
    if (flag<=2||flag==4)break;
    if ((flag&3)==3) dfs(x*p[j],y*p[j],z,j+1,-now);
    if ((flag&5)==5) dfs(x*p[j],y,z*p[j],j+1,-now);
    if ((flag&6)==6) dfs(x,y*p[j],z*p[j],j+1,-now);
    if (flag==7) dfs(x*p[j],y*p[j],z*p[j],j+1,2*now);
  }
}
void solve() {
  ret=0;
  dfs(1,1,1,1,1);
  printf("%d\n",(ret+mod)%mod);
}
int main()
{
  sieve(100000);
  int T;
  scanf("%d",&T);
  while(T--){
    scanf("%d%d%d",&A,&B,&C);
    solve();
  }
  return 0;
}

• 进一步拟合

  令 $q(ab,a,b)=H(a,a,1)H(b,1,b)$（实际上等于 $\mu(a)\mu(b)$），其余为零。令 $h_2=h/q$，$f=h_2*q*g$

  $$\begin{matrix} \begin{aligned} S_f(A,B,C)=\sum_{x=1}^A\sum_{y=1}^B\sum_{z=1}^Ch_2(x,y,z)S_{q*g}(\lfloor A/x\rfloor,\lfloor B/y\rfloor,\lfloor C/z\rfloor) \end{aligned}\\\\ \begin{aligned} S_{q*g}(A,B,C)&=\sum_{x,y}q(xy,x,y)S_{g}(\lfloor A/xy\rfloor,\lfloor B/x\rfloor,\lfloor C/y\rfloor )\\ &=\sum_{x,y}\mu(x)\mu(y)S_d(\lfloor A/xy\rfloor)S_d(\lfloor B/x\rfloor)S_d(\lfloor C/y\rfloor)\\ &=\sum_{x}\mu(x)S_d(\lfloor B/x\rfloor)\sum_{y}\mu(y)S_d(\lfloor A/xy\rfloor)S_F(\lfloor C/y\rfloor)\\ &=\sum_{x}\mu(x)S_d(\lfloor B/x\rfloor)W(\lfloor A/x\rfloor,C)\\ \end{aligned} \end{matrix}$$

  其中 $W(A,C)=\sum_{y}\mu(y)S_d(\lfloor A/y\rfloor)S_d(\lfloor C/y\rfloor)$，有用的 $A,C$ 总是整除位置，状态量 $O(N)$。

  $W$ 可以整除分块求解。复杂度 $O(\sqrt{A}+\sqrt{C})$。如此求解所有需要的 $W$，复杂度是

  $$\sum_{x=1}^{\sqrt{N}}\sum_{y=1}^{\sqrt{N}}\sqrt{N/x}+\sqrt{N/y}\sim N\left(\sum_{x=1}^{\sqrt{N}}x^{-1/2}\right)=O(N^{5/4})$$

  预处理 $W$ 之后，$S_{q*g}(A,B,C)$ 可以用 $O(B)$ 枚举计算（甚至不需要整除分块）。根据对称性，$S_{q*g}(A,B,C)=S_{q*g}(A,C,B)$，所以复杂度可以做到 $O\big(\min(B,C)\big)$。

• $h_2$ 密度的分析

  $h$ 在除以 $q$ 之后，进一步使得 $h_2(q,q,1),h_2(q,1,q)$ 变成了零，只有 $h_2(1,q,q)$ 仍然非零。

  对于 $h(x,y,z)$首先枚举可行的 $(y,z)$，再确定 $x$ 的方案数。

  类似讲义中二元的情况，$y,z$ 的每个素数幂必须同时为 $0$ 或同时非零。

  记 $s(n)$ 为满足 $xy=n$ 的合法 $(x,y)$ 对子数，$s$ 是奇函数且 $s(p^k)=k-1$。

  将多项式打表可以发现，$h_2(q^i,q^j,q^k)\neq 0$ 则 $i\leq j+k$，这说明 $h_2(x,y,z)\neq 0$ 则 $x|yz$。对于一个确定的 $(y,z)$，$x$ 的方案数不超过 $d(yz)$。于是我们关心的是 $s'=s\cdot d$ 的前缀和。

  $s'$ 满足 $s'(p^k)=(k+1)(k-1)$，根据引理，$S_{s'}(n)=O(\sqrt{n}\log^2 n)$。实际上 $n=N^2$，非零 $h_2$ 在 $N\times N\times N$ 中的数目为 $S_{s'}(N^2)=O(N\log^2N)$。

• 对每个 $h_2$ 计算 $S_{q*g}$ 的开销

  $$\begin{aligned} &\sum_{x=1}^A\sum_{y=1}^B\sum_{z=1}^C[h_2(x,y,z)\neq 0]\min(B/y,C/z)\\ &=\sum_{y=1}^B\sum_{z=1}^C\min(B/y,C/z)\sum_{x=1}^A[h_2(x,y,z)\neq 0]\\ &\leq \sum_{y=1}^B\sum_{z=1}^C\min(B/y,C/z)(s\cdot d')(yz)\\ &=\sum_{y=1}^B\sum_{z=1}^C(s\cdot d')(yz)\sum_{t\leq B/y,C/z}1\\ &=\sum_{t=1}^{\min(B,C)}\sum_{y=1}^{B/t}\sum_{z=1}^{C/t}(s\cdot d')(yz)\\ &\leq \sum_{t=1}^{\min(B,C)}S_{s\cdot d'}(BC/t^2)\\ &\sim \sum_{t=1}^{\min(B,C)}\dfrac{\sqrt{BC}}{t^2}\log^2N\\ &\sim N\log^3N\\ \end{aligned}$$

综上，我们证明了算法的复杂度为 $O(N^{5/4}+N\log^3N)$，代码如下：

#include<algorithm>
#include<cstring>
#include<cstdio>
#include<cmath>
#define MaxN 100050
#define MaxS 320
using namespace std;
const int mod=1000000007, K=17;
int p[MaxN>>3],tn,pk[MaxN>>3][K+1],
    d[MaxN],mu[MaxN],smu[MaxN],
    H[K+1][K+1][K+1];
void sieve(int n)
{
  d[1]=mu[1]=1;
  for (int i=2;i<=n;i++) {
    if (!d[i]) {
      d[p[++tn]=i]=2;
      mu[i]=-1;
      pk[tn][0]=1;
      for (int j=1;;j++){
        if (1ll*pk[tn][j-1]*i>n)break;
        pk[tn][j]=pk[tn][j-1]*i;
      }
    }
    for (int j=1,t;j<=tn&&(t=i*p[j])<=n;j++) {
      if (i%p[j]==0) {d[t]=2*d[i]-d[i/p[j]];break;}
      d[t]=2*d[i];
      mu[t]=-mu[i];
    }
  }
  for (int i=1;i<=n;i++) {
    d[i]+=d[i-1];
    smu[i]=smu[i-1]+mu[i];
  }
  H[0][0][0]=1;
  H[1][1][0]=H[1][0][1]=H[0][1][1]=-1;
  H[1][1][1]=2;
  for (int i=1;i<=K;i++)
    for (int j=0;j<=K;j++)
      for (int k=0;k<=K;k++) {
        if (j>0) H[i][j][k]+=H[i-1][j-1][k];
        if (k>0) H[i][j][k]+=H[i-1][j][k-1];
        if (i>1&&j>0&&k>0) H[i][j][k]-=H[i-2][j-1][k-1];
      }
}
int getk(int x, int p){
  int ret=0;
  while(x>=p){x/=p;ret++;}
  return ret;
}
int calcS0(int A, int C) {
  if (A<C)swap(A,C);
  int ret=0,l=1,r;
  while(l*l<=A){
    if (mu[l]) ret=(ret+1ll*mu[l]*d[A/l]*d[C/l])%mod;
    l++;
  }
  for (;l<=C;l=r+1){
    r=min(A/(A/l),C/(C/l));
    ret=(ret+1ll*(smu[r]-smu[l-1])*d[A/l]%mod*d[C/l])%mod;
  }
  return ret;
}
int SAB[MaxS*2][MaxS*2],SAC[MaxS*2][MaxS*2],
    savA,savB,savC,sqrtA,sqrtB,sqrtC;
inline int getSAB(int A, int B) {
  int i=(A<=sqrtA) ? A+MaxS : savA/A,
      j=(B<=sqrtB) ? B+MaxS : savB/B;
  return SAB[i][j] ? SAB[i][j] : SAB[i][j]=calcS0(A,B);
}
inline int getSAC(int A, int C) {
  int i=(A<=sqrtA) ? A+MaxS : savA/A,
      j=(C<=sqrtC) ? C+MaxS : savC/C;
  return SAC[i][j] ? SAC[i][j] : SAC[i][j]=calcS0(A,C);
}
int S(int A, int B, int C) {
  int ret=0;
  if (A<C||B<C) {
    int lim=min(A,B);
    for (int x=1;x<=lim;x++)
      if (mu[x])
        ret=(ret+1ll*mu[x]*d[B/x]*getSAC(A/x,C))%mod;
    return ret;
  }
  for (int x=1;x<=C;x++)
    if (mu[x])
      ret=(ret+1ll*mu[x]*d[C/x]*getSAB(A/x,B))%mod;
  return ret;
}
int ret;
void dfs(int a, int b, int c, int k, int now) {
  ret=(ret+1ll*now*S(a,b,c))%mod;
  for (int j=k;j<=tn;j++) {
    if (b<p[j]||c<p[j])break;
    int ka=getk(a,p[j]), kb=getk(b,p[j]), kc=getk(c,p[j]);
    for (int x=0;x<=ka;x++)
      for (int y=1;y<=kb&&y<=x+1;y++)
        for (int z=max(x-y,1);z<=kc&&y+z<=x+2;z++)
          dfs(a/pk[j][x], b/pk[j][y], c/pk[j][z],
              j+1, 1ll*now*H[x][y][z]%mod);
  }
}
void solve(int A, int B, int C) {
  ret=0;
  savA=A;sqrtA=sqrt(A);
  savB=B;sqrtB=sqrt(B);
  savC=C;sqrtC=sqrt(C);
  memset(&SAB[0][0],0,sizeof(SAB));
  memset(&SAC[0][0],0,sizeof(SAC));
  dfs(A,B,C,1,1);
  printf("%d\n",(ret+mod)%mod);
}
int main()
{
  int T,A,B,C;
  sieve(100000);
  scanf("%d",&T);
  while(T--){
    scanf("%d%d%d",&A,&B,&C);
    solve(A,B,C);
  }
  return 0;
}

两个算法的运行时间如下图所示。由于常数较大，新算法在本题的数据范围下并无优势。

蒟蒻语

第 50 道黑题，写篇题解祭一下 /cy

蒟蒻解

$\begin{aligned} &\sum\limits_{i=1}^{A} \sum\limits_{j=1}^{B} \sum\limits_{k=1}^{C} d(i, j, k)\\ &= \sum\limits_{i=1}^{A} \sum\limits_{j=1}^{B} \sum\limits_{k=1}^{C} \sum\limits_{x|i} \sum\limits_{y|j} \sum\limits_{z|k} [\gcd(x, y) = 1][\gcd(y, z) = 1][\gcd(z, x) = 1]\\ &= \sum\limits_{x=1}^{A} \sum\limits_{y=1}^{B} \sum\limits_{z=1}^{C} [\gcd(x, y) = 1][\gcd(y, z) = 1][\gcd(z, x) = 1] \left\lfloor\frac{A}{x}\right\rfloor \left\lfloor\frac{B}{y}\right\rfloor \left\lfloor\frac{C}{z}\right\rfloor\\ &= \sum\limits_{d_1=1}^{A} \mu(d_1)\sum\limits_{d_2=1}^{B} \mu(d_2) \sum\limits_{d_3=1}^{C} \mu(d_3) \sum\limits_{x=1}^{ \left\lfloor\frac{A}{\operatorname{lcm}(d_1,d_2)}\right\rfloor} \left\lfloor\frac{A}{\operatorname{lcm}(d_1,d_2)x}\right\rfloor \sum\limits_{y=1}^{ \left\lfloor\frac{B}{\operatorname{lcm}(d_2,d_3)}\right\rfloor} \left\lfloor\frac{B}{\operatorname{lcm}(d_2,d_3)y}\right\rfloor \sum\limits_{z=1}^{ \left\lfloor\frac{C}{\operatorname{lcm}(d_3,d_1)}\right\rfloor} \left\lfloor\frac{C}{\operatorname{lcm}(d_3,d_1)z}\right\rfloor\\ \end{aligned}$

(设 $Max = \max (A, B, C)$)

考虑哪些三元组会有贡献。

1. $\mu(d_1), \mu(d_2), \mu(d_3)$ 不为 $0$。
2. $\operatorname{lcm}(d_1,d_2) , \operatorname{lcm}(d_2,d_3), \operatorname{lcm}(d_3,d_1) \le Max$

把一个数看成一个点。

首先我们把 $\mu(x)$ 为 $0$ 的点排除在外。

然后把 $\operatorname{lcm}(x, y) \le Max$ 的点连接起来。

然后再做一遍三元环计数就好了，暴力枚举每一个三元环统计答案。为了更快地统计答案，可以先把边定向。

这样的时间复杂度是 $O(m \sqrt m)$ 的。( $m$ 是边数)

事实上，边数最多只有 $821535$ 条 (用蒟蒻劣质的程序算出来)

因此这道题就被解决了

蒟蒻码

#include<bits/stdc++.h>
#define L(i, j, k) for(int i = j; i <= k; i++)
#define R(i, j, k) for(int i = j; i >= k; i--)
#define ll long long
using namespace std;
const int N = 1e5 + 7;
const int M = 1e6 + 7;
const int mod = 1e9 + 7;
bool Prime[N];
int tot, P[N >> 3], mu[N];
void xxs() {
	mu[1] = 1;
	L(i, 2, 1e5)  {
		if(!Prime[i]) P[++tot] = i, mu[i] = -1;
		for(int j = 1; P[j] * i <= 1e5 && j <= tot; j++) {
			Prime[P[j] * i] = 1;
			if(i % P[j] == 0) {
				mu[P[j] * i] = 0;
				break;
			}
			mu[P[j] * i] = -mu[i];
		}
	}
}
vector<int> ve[N];
unordered_map<int, bool> mp[N];
int T, Max, A, B, C, ans, p[N];
int sA[M], sB[M], sC[M], cnt, deg[N];
int gcd(int x, int y) { return x == 0 ? y : gcd(y % x, x); }
ll lcm(int x, int y) { return 1ll * x / gcd(x, y) * y; }
int ta, tb, tc, sum;
void dh(int a, int b, int c) {
	(sum += 1ll * p[a / ta] * p[b / tb] % mod * p[c / tc] % mod) %= mod;
}
void get(int a, int b, int c) {
	sum = 0;
	if(a == b && b == c) dh(A, B, C);
	else if(a == b || b == c || c == a) dh(A, B, C), dh(C, A, B), dh(B, C, A);
	else dh(A, B, C), dh(A, C, B), dh(B, A, C), dh(B, C, A), dh(C, A, B), dh(C, B, A);
	(ans += (mu[a] * mu[b] * mu[c] * sum % mod + mod) % mod) %= mod;
}
int vis[N];
void js() {
	L(i, 1, cnt) {
		int &u = sA[i], &v = sB[i];
		if(deg[u] > deg[v]) swap(u, v);
		else if(deg[u] == deg[v] && u > v) swap(u, v);
		ve[u].push_back(i);
	}
	L(i, 1, Max) {
		for(int j : ve[i]) vis[sB[j]] = sC[j];
		for(int j : ve[i]) for(int k : ve[sB[j]])
			if(vis[sB[k]]) ta = vis[sB[k]], tb = sC[j], tc = sC[k], get(i, sB[j], sB[k]);
		for(int j : ve[i]) vis[sB[j]] = 0;
	}
}
void getans(int x) {
	for(int l = 1, r; l <= x; l = r + 1) {
		r = (x / (x / l));
		(p[x] += 1ll * (r - l + 1) * (x / l) % mod) %= mod;
	}
}
int main() {
	scanf("%d", &T);
	xxs();
	L(i, 1, 1e5) getans(i);
	while(T--) {
		scanf("%d%d%d", &A, &B, &C);
		ans = cnt = 0;
		Max = max(max(A, B), C);
		R(i, Max, 1) for(int j = i; j <= Max; j += i) {
				if(!mu[j]) continue;
				for(int k = j; k <= 1ll * Max * i / j; k += i) {
					if(!mu[k]) continue;
					if(!mp[j][k]) ++cnt, sA[cnt] = j, sB[cnt] = k, sC[cnt] = j / i * k, deg[j] ++, deg[k] ++, mp[j][k] = 1;
				}
			}
		js();
		L(i, 1, Max) ve[i].clear(), mp[i].clear(), deg[i] = 0;
		printf("%d\n", (ans + mod) % mod);
	}
	return 0;
}