LuoguP5176 公约数 题解

题意：
求 (1000) 次

[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) ]

(n、m、p) 规模均为 (2e7)。

---

由于我蒻的一匹， 所以我看了题解。

题目的第一步是想到（或是证出来）一个结论， 如下：

[gcd(i*j,i*k,j*k) = frac{gcd(i,j)*gcd(i,k)*gcd(j,k)}{gcd(i,j,k)} ]

怎么证这个结论呢？

首先就是将式子里的所有项都写成质因数分解的形式。

[gcd(i*j,i*k,j*k) = sum_{p=1}^{infty} p_i^{min(c_{i1}+c_{i2},c_{i1}+c_{i3},c_{i2}+c_{i3})} ]

[gcd(i,j)*gcd(i,k)*gcd(j,k) = sum_{p=1}^{infty} p_i^{min(c_{i1},c_{i2})+min(c_{i1},c_{i3})+min(c_{i2},c_{i3})} ]

[gcd(i,j,k) = sum_{p=1}^{infty} p_i^{min(c_{i1},c_{i2},c_{i3})} ]

于是

[gcd(i*j,i*k,j*k) = frac{gcd(i,j)*gcd(i,k)*gcd(j,k)}{gcd(i,j,k)} ]

(Leftarrow min(c_{i1}+c_{i2},c_{i1}+c_{i3},c_{i2}+c_{i3}) = min(c_{i1},c_{i2})+min(c_{i1},c_{i3})+min(c_{i2},c_{i3}) - min(c_{i1},c_{i2},c_{i3}))

似乎还没有规定 (c_{x1}) 、 (c_{x2}) 、 (c_{x3}) 分别对应 (i) 、 (j) 、 (k) 中的哪一个， 但是无论 (c_{x1}) 、 (c_{x2}) 、 (c_{x3}) 分别对应 (i) 、 (j) 、 (k) 中的哪一个，
(min(c_{i1}+c_{i2},c_{i1}+c_{i3},c_{i2}+c_{i3})) 的值
和 (min(c_{i1},c_{i2})+min(c_{i1},c_{i3})+min(c_{i2},c_{i3}) - min(c_{i1},c_{i2},c_{i3})) 的值
都不会改变（因为它们都是轮换式）。

故规定 (c_{i1} leq c_{i2} leq c_{i3}), 原等式转为

[c_{i1} + c_{i2} = c_{i1} + c_{i1} + c_{i2} - c_{i1} ]

故结论成立。

---

所以

[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) ]

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

展开

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

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

把变量换得整齐点。

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

然后题目就转化为了分别求三个同样形式的式子。

那么怎么求(sum_{i=1}^nsum_{j=1}^m gcd^2(i,j))？

套路上场~~~

设(f(n) = n^2)， 则

(sum_{i=1}^nsum_{j=1}^m gcd^2(i,j) = sum_{i=1}^nsum_{j=1}^m f(gcd(i,j)))

要找出(1 * g = f)， 则(g = f * mu)， 可以 (O(nlog log n)) 卷出来， 于是

[sum_{i=1}^nsum_{j=1}^m f(gcd(i,j)) ]

[= sum_{i=1}^nsum_{j=1}^m sum_{d|gcd(i,j)} g(d) ]

[= sum_{d=1}^{min(n,m)} g(d) lfloor frac{n}d floor lfloor frac{m}d floor ]

分块求即可。

Luogu数据吸氧才能AC的丑陋代码
(其实函数 (g) 可以线性筛出来， 但是我不想写qwq)

#include<bits/stdc++.h>
using namespace std;
const int maxn = 20000001;
const int mod = 1e9 + 7;
#define li long long

int pr_cnt, prime[maxn], v[maxn];
li g[maxn];
void euler(int n) {
	g[1] = 1ll;
	for(register int i=2; i<=n; ++i) {
		if(!v[i]) {
			v[prime[++pr_cnt] = i] = i;
			g[i] = 1ll * i * i % mod;
		}
		for(int j=1; j<=pr_cnt; ++j) {
			if(prime[j] > n/i || prime[j] > v[i]) break;
			v[prime[j] * i] = prime[j];
			g[prime[j] * i] = 1ll * g[i] * g[prime[j]] % mod;
		}
	}
	for(int i=1; i<=pr_cnt; ++i)
		for(int j=n/prime[i]; j>0; --j) {
			g[prime[i] * j] -= g[j];
			g[prime[i] * j] %= mod;
		}
	for(int i=1; i<=n; ++i) {
		g[i] += g[i-1];
		g[i] %= mod;
	}
}

li calc(li n, li m) {
	li len = min(n, m);
	li res = 0ll;
	for(register li i=1,j; i<=len; i=j+1) {
		j = min(n/(n/i), m/(m/i));
		j = min(j, len);
		res += 1ll * (g[j]-g[i-1]) * (n/i) % mod * (m/i) % mod;
		res %= mod;
	}
	return res;
}

int main()
{
	euler(20000000);
	int T; cin >> T; while(T--) {
		li n,m,p;
		scanf("%lld%lld%lld", &n,&m,&p);
		li ans = 0ll;
		ans = 1ll * p * calc(n,m) % mod;
		ans += 1ll * n * calc(m,p) % mod;
		ans %= mod;
		ans += 1ll * m * calc(n,p) % mod;
		ans = (ans%mod+mod)%mod;
		cout << ans << '
';
	}
	return 0;
}