公约数

题目地址
看到这题目就不想做了系列，出题人是不是都不知道杜教筛是什么东西啊，他家杜教筛可以预处理优化到(O(n^frac{2}{3}))

先吓唬你一下，我们要求：

[sum_{i=1}^nsum_{j=1}^msum_{k=1}^pgcd(icdot j,icdot k,jcdot k) imes gcd(i,j,k) imes left(frac{gcd(i,j)}{gcd(i,k) imes gcd(j,k)}+frac{gcd(i,k)}{gcd(i,j) imes gcd(j,k)}+frac{gcd(j,k)}{gcd(i,j) imes gcd(i,k)} ight) ]

我们令(x = gcd(i, j), y = gcd(i, k), z = gcd(j, k))

于是我们有：

[sum_{i=1}^nsum_{j=1}^msum_{k=1}^pgcd(icdot j,icdot k,jcdot k) imes gcd(i,j,k) imes left(frac{x}{y imes z}+frac{y}{x imes z}+frac{z}{x imes y} ight) ]

[sum_{i=1}^nsum_{j=1}^msum_{k=1}^pgcd(icdot j,icdot k,jcdot k) imes gcd(i,j,k) imes frac{x ^ 2 + y ^ 2 + z ^ 2}{x imes y imes z} ]

展开第一陀恶心人的柿子，可以得到：

[gcd(icdot j,icdot k,jcdot k) = frac{x imes y imes z}{gcd(i, j, k)} ]

（证明可以通过唯一分解定理展开）

于是我们惊喜地发现，我们只需要求(()还原(x, y, z))：

[sum_{i=1}^nsum_{j=1}^msum_{k=1}^p gcd(i, j)^2 +gcd(j, k) ^ 2 + gcd(i, k) ^ 2 ]

我们把柿子分开写，得到：

[sum_{i=1}^nsum_{j=1}^msum_{k=1}^p gcd(i, j)^2+sum_{i=1}^nsum_{j=1}^msum_{k=1}^pgcd(j, k) ^ 2 +sum_{i=1}^nsum_{j=1}^msum_{k=1}^p gcd(i, k) ^ 2 ]

因为(gcd(i, j))与k无关（其他同理），于是我们可以提出k，变成p个(sum_{i=1}^nsum_{j=1}^m gcd(i, j)^2)相加，故有：

[p imes sum_{i=1}^nsum_{j=1}^m gcd(i, j)^2+n imes sum_{j=1}^msum_{k=1}^pgcd(j, k) ^ 2 +m imes sum_{i=1}^nsum_{k=1}^p gcd(i, k) ^ 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)）

[F(n, m) = sum_{d = 1}^nsum_{i=1}^nsum_{j=1}^m d^2 imes [gcd(i, j) == d] ]

(i, j)同时除以d得：

[F(n, m) = sum_{d = 1} ^ nsum^{n / d}_{i = 1}sum^{m / d}_{j = 1} d ^ 2 imes [gcd(i, j) == 1] ]

[F(n, m) = sum_{d = 1} ^ nd ^ 2sum^{n / d}_{i = 1}sum^{m / d}_{j = 1}sum_{k|gcd(i, j)} mu(k) ]

提前枚举K得：

[F(n, m) = sum_{d = 1} ^ n d ^ 2sum_{k = 1}^{n / d}mu(k)sum^{n / (d * k)}_{i = 1}sum^{m / (d * k)}_{j = 1} ]

[F(n, m) = sum_{d = 1} ^ n d ^ 2sum_{k = 1}^{n / d}mu(k)left lfloor frac{n}{d * k} ight floor imes left lfloor frac{m}{d * k} ight floor ]

跟据反演套路，我们令(T = d * k)

[F(n, m) = sum_{T = 1}^n left lfloor frac{n}{T} ight floor imes left lfloor frac{m}{T} ight floorsum_{d | T} d ^ 2 imes mu(frac{T}{d}) ]

令(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;
}