官方题解只能做到 $O(N^{1.5} \log^3 N)$？玩的太差了。记 $N=\max(A,B,C)$。这道题可以做到 $O(N^{1+o(1)})$ 时间、$O(N)$ 空间，不必停留在三元环枚举。更准确地说，下面算法的时间复杂度是 $O(N\log^3N\cdot\Psi(N))$，其中 $\Psi(N)=N^{o(1)}$；这个额外因子的含义和证明会在后面给出。
题目要求计算 $\sum_{i\le A}\sum_{j\le B}\sum_{k\le C}\tau(ijk)$，对 $10^9+7$ 取模，其中 $\tau$ 表示约数个数函数，且 $A,B,C\le10^5$。官方题解将其转化为最小公倍数图上的三元环计数，给出的单组复杂度是 $O(N^{3/2}\log^3N)$。下面采用另一种数论展开，再把一维枚举整体合并。([QOJ][1])
数论展开与整除分块
先定义约数个数的前缀和 $D(x)=\sum_{i=1}^{x}\tau(i)$，其中 $D(0)=0$，并记 $f_X(u)=D(\lfloor X/u\rfloor)$。所有需要的 $D(x)$ 都可以在线性筛求出 $\tau$ 后，取一次前缀和得到。
考虑一个质数 $p$，设它在 $i,j,k$ 中的指数分别为 $\alpha,\beta,\gamma$。它对 $\tau(ijk)$ 的贡献为 $\alpha+\beta+\gamma+1$。
这个贡献可以这样展开：先取 $(\alpha+1)(\beta+1)(\gamma+1)$，减去 $\alpha\beta(\gamma+1)$、$\alpha\gamma(\beta+1)$、$\beta\gamma(\alpha+1)$，最后加上 $2\alpha\beta\gamma$。直接展开多项式即可验证结果确实为 $\alpha+\beta+\gamma+1$。
把这个恒等式对所有质数相乘。每个质数有五种选择：不进行修正；同时从 $i,j$ 中除去一个 $p$，系数为 $-1$；同时从 $i,k$ 中除去一个 $p$，系数为 $-1$；同时从 $j,k$ 中除去一个 $p$，系数为 $-1$；同时从三者中除去一个 $p$，系数为 $2$。
分别把后三种“两两修正”中选中的质数乘积记为 $r,s,t$，把“三者同时修正”中选中的质数乘积记为 $d$。每个质数最多属于其中一类，所以 $drst$ 必须无平方因子。注意，这不仅要求四个数各自无平方因子，还要求它们两两互质。
设 $\omega(d)$ 表示不同质因子的个数。将展开式代入原来的三重求和，并交换求和顺序，得到整个算法的核心恒等式：
$$ F(A,B,C)= \sum_{\substack{d,r,s,t\ge1\\ \mu^2(drst)=1}} 2^{\omega(d)}\mu(r)\mu(s)\mu(t)\, f_A(drs)f_B(drt)f_C(dst). $$
这里，固定 $d,r,s,t$ 后，剩余三个变量已经完全独立。例如，关于 $i$ 的部分就是 $\sum_{drs\mid i,\ i\le A}\tau(i/(drs))=D(\lfloor A/(drs)\rfloor)$。
直接枚举这个式子仍然不够快。接下来利用对称性，只枚举 $r\le s\le t$，再把相应的不同排列一起计算。同时，由于原问题对三个上界对称，不妨先将它们排序为 $A\le B\le C$。
当 $r < s < t$ 时，三个乘积满足 $drs < drt< dst$。把 $r,s,t$ 的六种排列分别代入核心恒等式中的三个 $f$，将得到的六个乘积之和记为 $W_d(r,s,t)$。它们的系数 $2^{\omega(d)}\mu(r)\mu(s)\mu(t)$ 完全相同。
存在非零贡献的充要条件是 $drs\le A$、$drt\le B$、$dst\le C$。这是因为三个乘积与三个上界都已排序，能够匹配当且仅当对应位置满足大小关系。因此，固定 $d,r,s$ 后，只需要考虑
$t\le U=\min(\lfloor B/(dr)\rfloor,\lfloor C/(ds)\rfloor)$。
重复变量需要单独注意。由于 $drst$ 无平方因子，两个相等的变量只能都等于 $1$。所以，除 $r=s=t=1$ 外，只可能出现 $r=s=1< t$，此时应计算三个不同排列，而不是六个。$r=s=t=1$ 的贡献直接累加 $2^{\omega(d)}f_A(d)f_B(d)f_C(d)$ 即可；其余情况统一从 $t=s+1$ 开始枚举。
现在固定 $d,r,s$，并要求 $\mu(drs)\ne0$。与 $t$ 有关的约数前缀和，只涉及以下六个整除商：
$\lfloor A/(drt)\rfloor$、$\lfloor B/(drt)\rfloor$、$\lfloor C/(drt)\rfloor$，以及 $\lfloor A/(dst)\rfloor$、$\lfloor B/(dst)\rfloor$、$\lfloor C/(dst)\rfloor$。
因此，可以按照这六个整除商同时不变的极大区间，对 $t$ 做整除分块。 在一个块 $[L,R]$ 内，$W_d(r,s,t)$ 是常数，不必逐个计算每个 $t$ 的贡献。
具体地，将六个分子分别预先写成 $\lfloor X/(dr)\rfloor$ 和 $\lfloor X/(ds)\rfloor$，其中 $X\in{A,B,C}$。对于其中一个分子 $q$，当前商为 $v=\lfloor q/L\rfloor$。若 $v>0$，它给出的右端点限制是 $\lfloor q/v\rfloor$；若 $v=0$，则不再限制右端点。取所有限制与 $U$ 的最小值即可。
令 $m=drs$。块内剩下唯一需要计算的量是 $\sum_{t=L}^{R}\mu(t)[\gcd(t,m)=1]$：$\mu(t)$ 已经自动处理了 $t$ 含有平方因子的情况，而互质条件保证 $t$ 与 $d,r,s$ 没有公共质因子。
问题因此变成了：如何快速计算带互质限制的莫比乌斯函数区间和？
互质莫比乌斯和与复杂度
定义普通莫比乌斯前缀和 $M(x)=\sum_{i=1}^{x}\mu(i)$，以及限制互质的前缀和 $M_m(x)=\sum_{i=1}^{x}\mu(i)[\gcd(i,m)=1]$。$M(x)$ 可以直接在线性筛后预处理。
我们有下面的恒等式：
$$ M_m(x)= \sum_{\substack{q\le x\\ \operatorname{rad}(q)\mid m}} M\!\left(\left\lfloor\frac{x}{q}\right\rfloor\right), $$
其中 $\operatorname{rad}(q)$ 表示 $q$ 的不同质因子的乘积。换句话说，求和中的 $q$ 只能使用 $m$ 的质因子，但这些质因子的指数可以任意大。
这里不能只枚举 $m$ 的约数。 例如，当 $m=6$ 时，需要考虑所有不超过 $x$ 的 $2^a3^b$，而不仅仅是 $1,2,3,6$。
证明也很直接。令 $h_m(q)=[\operatorname{rad}(q)\mid m]$，则狄利克雷卷积满足 $(\mu*h_m)(t)=\mu(t)[\gcd(t,m)=1]$。对于 $p\mid m$，$h_m(p^e)$ 恒为 $1$，卷积后所有正指数项都消掉；对于 $p\nmid m$，则保留原来的莫比乌斯函数系数。对这个卷积恒等式求前缀和，就得到上述公式。
因此，只要分解 $m$，然后用深度优先搜索枚举这些质因子的幂次，就可以计算 $M_m(x)$。分解使用预处理的最小质因子数组即可，不需要试除。
区间和可以直接同时计算两端：将前缀端点 $(L-1,R)$ 传入递归，每次枚举一个质数的幂，并将两个端点同时整除该质数幂。当没有剩余质数，或者最小的剩余质数已经大于右端点时，返回两个普通 $M$ 值之差。两个端点相等时，直接返回 $0$。
下面严格分析这个做法，而不是把一次互质莫比乌斯和查询当成常数时间。
先计算整除块的总数。固定 $d,r,s$ 后，从 $t=s+1$ 开始，六个整除商的初值都不超过 $N/(drs)$。它们都是单调不增的非负整数，因此共同划分出的块数为 $O(N/(drs))$。
此外，$r\le s\le t$ 且 $dst\le N$，意味着 $ds^2\le N$。即使忽略所有互质限制、无平方因子限制以及更紧的上界，整除块总数仍然不超过
$$ O\!\left( \sum_{d\le N} \sum_{s\le\sqrt{N/d}} \sum_{r\le s} \frac{N}{drs} \right) = O(N\log^3N). $$
再分析每个块的查询代价。记 $\Psi_m(x)=\#\{q\le x:\operatorname{rad}(q)\mid m\}$，以及 $\Psi(N)=\max_{1\le m,x\le N}\Psi_m(x)$。
上述递归的代价为 $O(\Psi_m(x))$。实现时按质数从小到大枚举，并在剩余质数超过右端点时立即结束，那么每个未终止的递归节点至少有两个子节点，而叶子数不超过需要枚举的整数个数，总节点数也就与之同阶。
关键是：当 $m,x\le N$ 时，$\Psi_m(x)$ 一致满足 $N^{o(1)}$ 的上界。
任取固定的 $\varepsilon>0$，利用几何级数，有 $\Psi_m(x)\le x^\varepsilon\prod_{p\mid m}(1-p^{-\varepsilon})^{-1}$。对于充分大的质数 $p$，有 $(1-p^{-\varepsilon})^{-1}\le p^\varepsilon$；不满足这个条件的质数只有有限多个，它们的贡献可以吸收到只依赖于 $\varepsilon$ 的常数 $C_\varepsilon$ 中。于是
$\Psi_m(x)\le C_\varepsilon x^\varepsilon\operatorname{rad}(m)^\varepsilon\le C_\varepsilon N^{2\varepsilon}$。
由于 $\varepsilon$ 可以任意小，便有 $\Psi(N)=N^{o(1)}$。因此，最终时间复杂度为 $O(N\log^3N\cdot\Psi(N))=O(N^{1+o(1)})$。
枚举 $d,r,s$ 本身只有 $O(N\log N)$ 组候选，为每组分解一次 $drs$ 的代价也不会超过上述总界。预处理只需要线性筛 $\mu$、$\tau$、最小质因子以及 $\mu^2(d)2^{\omega(d)}$，再计算 $D$ 和 $M$ 的前缀和。整个算法不需要保存所有询问状态，也不需要建图，空间复杂度为 $O(N)$。
多组数据时，设 $N_{\max}$ 为所有数据的最大上界，则只需一次 $O(N_{\max})$ 预处理，总时间为 $O(N_{\max}+\sum_i N_i^{1+o(1)})$，空间为 $O(N_{\max})$。这是一个经过证明的算法上界，并不依赖于“本题规模下三元环数量不多”之类的实验观察，也不将它宣称为已经证明最优的复杂度。
实现中还可以把很短的块直接枚举。例如，对于长度不超过 $32$ 的块，逐项判断 $\mu(dst)\ne0$ 且 $\mu(rt)\ne0$，满足时累加 $\mu(t)$。由于固定的 $drs$ 已经无平方因子，这两个条件恰好完成所需检查；同时 $dst\le N$、$rt\le N$，可以直接访问筛表，完全不需要计算 $\gcd$。这个优化不改变复杂度上界，但能减少递归常数。
