Statement：U119598 福利题
考虑将 \(x, y, z\) 唯一分解。
记我们当前讨论的质因子为 \(p\)，记 \(x_i, y_i, z_i\) 为 \(x,y,z\) 的唯一分解中 \(p\) 项对应的指数，则原式中六个项的 \(\text{lcm}\) 的唯一分解中 \(p\) 的指数为
\[\max \left\{ \begin{array}{ll} \max(x_i,z_i) + \min(y_i,z_i) - \min(x_i,y_i) \\ \max(x_i,y_i) + \min(x_i,z_i) - \min(y_i,z_i) \\ \max(y_i,z_i) + \min(x_i,y_i) - \min(x_i,z_i) \\ \max(x_i,z_i) +\min(x_i,y_i) - \min(y_i,z_i) \\ \max(x_i,y_i) + \min(y_i,z_i) - \min(x_i,z_i)\\\max(y_i,z_i) + \min(x_i,z_i) - \min(x_i,y_i) \end{array} \right.
\]
不妨设 \(x_i \leq y_i \leq z_i\)，则上式最大值为 \(z_i + y_i - x_i\)
回归原式，这个指数的形式等同于 \(\frac{xyz}{\gcd(x,y,z)^2}\)，因此原式可以被化为
\[\sum_{x = 1} ^ n \sum_{y = 1} ^ n \sum_{z = 1} ^ n \frac{xyz}{\gcd(x,y,z) ^ 2}
\]
提出 \(\gcd\)，得
\[\sum_{d = 1} ^ n d \sum_{i = 1} ^ {n/d} \sum_{j = 1} ^ {n/d}\sum_{k = 1} ^ {n/d}[\gcd(i, j, k) = 1]ijk
\]
\[\sum_{d = 1} ^ n d \sum_{p = 1} ^ {n / d}\mu(p)\sum_{i = 1} ^ {n/dp} \sum_{j = 1} ^ {n/dp}\sum_{k = 1} ^ {n/dp}p^3ijk
\]
转而记 \(L = dp\)，记 \(S(n) = \sum \limits_{i = 1} ^ n\sum \limits_{j = 1} ^ n\sum \limits_{k = 1} ^ n ijk = (\frac{n(n + 1)}{2})^3\) 得
\[\sum_{L = 1} ^ n S(\frac{n}{L})L \sum_{p | L} \mu(p)p^2
\]
对 \(\sum_{p | L} \mu(p)p^2\) 线性筛即可 \(\mathcal{O}(n)\) 解决。
后记：理论上界即使 int128 也是会炸的，但是本题中常数好像挺小，最大答案仅 \(10^{35}\) 量级。
#include <bits/stdc++.h>
using namespace std;
using i64 = long long;
using i128 = __int128;
constexpr int N = 1e6 + 10;
inline void write(i128 x){
if(x > 9) write(x / 10);
cout << (char)(x % 10 + '0');
}
vector<int> pr;
i64 f[N]; int vis[N], low[N];
void sieve(){
f[1] = 1;
for(int i = 2; i < N; i++){
if(!vis[i]) pr.push_back(i), low[i] = i, f[i] = 1 - 1ll * i * i;
for(int j = 0; j < pr.size() && pr[j] * i < N; j++){
vis[pr[j] * i] = 1;
if(i % pr[j] == 0){
low[pr[j] * i] = low[i] * pr[j];
if(low[i] == i){
f[pr[j] * i] = f[i];
} else {
f[pr[j] * i] = f[i / low[i]] * f[low[pr[j] * i]];
}
break;
}
low[pr[j] * i] = pr[j];
f[pr[j] * i] = f[i] * f[pr[j]];
}
}
}
inline i128 S(int n){
i128 t = 1ll * n * (n + 1) / 2;
return t * t * t;
}
void solve(){
int n; cin >> n;
i128 ans = 0;
for(int i = 1; i <= n; i++){
ans += S(n / i) * i * f[i];
}
write(ans);
}
int main(){
ios::sync_with_stdio(false);
cin.tie(nullptr); cout.tie(nullptr);
sieve();
solve();
}