题目地址
看到这题目就不想做了系列，出题人是不是都不知道杜教筛是什么东西啊，他家杜教筛可以预处理优化到(O(n^frac{2}{3}))
先吓唬你一下，我们要求：
我们令(x = gcd(i, j), y = gcd(i, k), z = gcd(j, k))
于是我们有：
展开第一陀恶心人的柿子，可以得到：
（证明可以通过唯一分解定理展开）
于是我们惊喜地发现，我们只需要求(()还原(x, y, z))：
我们把柿子分开写，得到：
因为(gcd(i, j))与k无关（其他同理），于是我们可以提出k，变成p个(sum_{i=1}^nsum_{j=1}^m gcd(i, j)^2)相加，故有：
设(F(n, m) = sum_{i=1}^nsum_{j=1}^m gcd(i, j)^2)
原式(=p imes F(n, m) + n imes F(p, m) + m imes F(n, p))
所以问题转化为如何快速的求出(F(n, m))（不妨设(n ≤ m)）
(i, j)同时除以d得：
提前枚举K得：
跟据反演套路，我们令(T = d * k)
令(g(x) = x * x, mu(x) = mu(x))
于是(F(n, m) = sum_{T = 1}^n left lfloor frac{n}{T} ight floor imes left lfloor frac{m}{T} ight floor imes(g(x) imes mu(x)))
发现后面那一堆就是一个卷积的形式，而且是两个积性函数的卷积，于是我们可以用线性筛筛出，筛法见于神之怒加强版
总结一下求法：
我们记(prim[i])为第i个质数，(low[i])为i质因数分解后最小的质因子的指数次幂，(f[x])为(g(x))卷(mu(x))，(mu[x])为(mu(x))
类似于线性筛，我们分三种情况讨论：
(case1:) (i)是质数
我们可以直接算出来：(low[i] = i,; f[i] = i imes i - 1,; mu[i] = -1)
(case2:i \% prim[j] != 0)：这种情况是满足积性函数定义的，即(i, prim[j])互质，令(p = i imes prim[j])，所以我们有(low[p] = prim[j],; mu[p] = -mu[i],; f[i] = f[i] imes f[prim[j]])
(case3:i \% prim[j] == 0)：这种情况不满足积性函数的定义，于是我们考虑用我们刚刚维护的(low[i])
考虑到(i / low[i],prim[j])是互质的，那我们可不可以类比(case2)通过(f[i / low[i]], f[prim[j]])来推出(f[p])呢？
假设(a, b)互质，我们有：(f[a imes b] = f[a] imes f[b])
所以我们也有：(f[p] = f[i] imes f[prim[j]] = f[i / low[i]] imes f[prim[j] imes low[i]])
根据定义，(i / low[i])没有(prim[j])这个约数，所以(i / low[i],; prim[j] imes low[i])互质
i肯定有(prim[j])这个质因数，而且根据线性筛的性质，(prim[j])肯定是i的最小质因数，故我们有(mu[p] = 0,; low[p] = low[i] imes prim[j])
如果(low[i] < i)，证明(low[i] imes prim[j] <= i)，故我们可以直接用(f[p] = f[i / low[i]] imes f[prim[j] imes low[i]])
如果(low[i] == i)，证明i以前肯定没有被筛过，若(low[i] = prim[j] ^ k)，则k肯定为质数，那么我们就可以直接用：(f[p] = (i * prim[j]) ^ 2 - i ^ 2)
(ps:)
写完后看题解才发现：这道题根于神之怒不一样，并不需要记录(low[i])，对于(case3)，(f[p] = f[i] imes prim[j] ^2)即可
(Code:)
#include<bits/stdc++.h>
using namespace std;
#define il inline
#define re register
il int read() {
re int x = 0, f = 1; re char c = getchar();
while(c < '0' || c > '9') { if(c == '-') f = -1; c = getchar();}
while(c >= '0' && c <= '9') x = x * 10 + c - 48, c = getchar();
return x * f;
}
#define rep(i, s, t) for(re int i = s; i <= t; ++ i)
#define mod 1000000007
#define maxn 20000007
il int Pow(int x) {return 1ll * x * x % mod;}
int cnt, prim[maxn], f[maxn], mu[maxn], low[maxn];
bool vis[maxn];
il void init() {
f[1] = mu[1] = 1;
rep(i, 2, 2e7) {
if(!vis[i]) f[i] = Pow(i) - 1, mu[i] = -1, low[i] = i, prim[++ cnt] = i;
for(re int j = 1; j <= cnt; ++ j) {
if(i * prim[j] > 2e7) break;
int p = i * prim[j];
vis[p] = 1;
if(i % prim[j] == 0) {
low[p] = low[i] * prim[j];
if(low[i] == i) f[p] = (Pow(1ll * i * prim[j] % mod) - Pow(i) + mod) % mod;
else f[p] = 1ll * f[i / low[i]] * f[low[i] * prim[j]] % mod;
break;
}
mu[p] = -mu[i], low[p] = prim[j], f[p] = 1ll * f[i] * f[prim[j]] % mod;
}
}
rep(i, 1, 2e7) f[i] = (f[i - 1] + f[i]) % mod;
}
il int F(int n, int m) {
if(n > m) swap(n, m);
int ans = 0;
for(re int l = 1, r; l <= n; l = r + 1) {
r = min(n / (n / l), m / (m / l));
ans = (ans + 1ll * (1ll * (n / l) * (m / l) % mod) * (f[r] - f[l - 1]) % mod) % mod;
}
return (ans + mod) % mod;
}
il int get(int n, int m, int p) {
return ((1ll * p * F(n, m) % mod + 1ll * n * F(p, m) % mod) % mod + 1ll * m * F(n, p)) % mod;
}
il void outit() {
for(re int i = 1, T = read(); i <= T; ++ i) printf("%d
", get(read(), read(), read()));
}
int main() {
init(), outit();
return 0;
}