「题解」[SDOI2015]约数个数和

link
又到了喜闻乐见的推式子环节

本题唯一的难点-引理：\(d(ij)=\sum_{u|i}\sum_{v|j}[\gcd(u,v)=1]\)。
证明：对于 \(ij\) 因数 \(k\)，令 \(c\) 为 \(k\) 唯一分解后某一个 \(prime\) 的指数，\(a\) 为 \(i\) 的指数，\(b\) 为 \(j\) 的指数。我们 钦定 \(a\) 优先贡献 \(c\)，剩下不够的用 \(b\) 补。那么一个 \(k\) 对应唯一的贡献情况，一种贡献情况也对应唯一的一个 \(k\)。

对于 \(c\leqslant a\) 情况，我们又 钦定 给甲方（\(i\)）乘上 \(prime^c\) 的价值。对于 \(c>a\)，钦定 给乙方（\(j\)）乘上 \(prime^{c-a}\) 的价值。这样两边的价值肯定互质，且 \(k\) 也是一一对应，不重不漏。容易发现两个价值就是原式的含义。得证。

但是 OI 题不是证明题。考虑如果我们原先不知道这个式子，怎么看出来。首先看如果 \(ij\) 有个因数为 \(6\)，那我可以 \(i=2,j=3\)，也可以 \(i=3,j=2\)，这可不好搞。

既然反正这个 \(6\) 都只能算一次，我们就只看成 \(2\) 和 \(3\) 是 \(i\) 提供的贡献。那么这样 \(j\) 可不认啊，你 \(i\) 有可能自己都不够这么多，凭什么提供这么多贡献。但是只给 \(j\) 也不行。

于是就自然想到 \(i,j\) 以一个规则提供贡献，很容易想到 \(i\) 先补完（情况 1），不够的 \(j\) 来补（情况 2）。

然后你想情况1的 \(j\) 不提供贡献和情况2的 \(i\) 全部提供完是一个东西。那么可以将情况2的 \(i\) 也视作没提供贡献。那么两边选的数就互质。那么稍微延伸一下这个式子就 yy 出来了。

然后就比较简单了。
\(\sum_{i=1}^n\sum_{j=1}^m\sum_{u|i}\sum_{v|j}[\gcd(u,v)=1]\)
推式子常用技巧：把因数转化成倍数。?
\(=\sum_{u=1}^n\sum_{v=1}^m \lfloor \frac n u \rfloor \lfloor \frac m v \rfloor[\gcd(u,v)=1]\)
此时两种思考，第一种是 Mobius 反演，第二种是直接利用 Mobius 函数的性质。这里展开讲第二种。
\(=\sum_{u=1}^n\sum_{v=1}^m \lfloor \frac n u \rfloor \lfloor \frac m v \rfloor\sum_{w|\gcd(u,v)}\mu_w\)
\(=\sum_{w=1}^{\min(n,m)}\mu_w\sum_{u=1}^{n/w}\sum_{v=1}^{m/w}\lfloor \frac n {wu} \rfloor \lfloor \frac m {wv} \rfloor\)
然后这道题就做完了。容易想到将 \(n/w\) 以及 \(m/w\) 整除分块，预处理前缀和搞一搞即可。

#include 
#include 
#include 
#include 
#define LL long long
using namespace std;
const int MAXN = 5e4 + 5;
int prime[MAXN / 10], cnt, mu[MAXN], n, m;
LL p[MAXN];
bool vis[MAXN];
void primes() {
	mu[1] = 1;
	for(int i = 2; i <= MAXN - 5; i ++) {
		if(vis[i] == 0) prime[++ cnt] = i, mu[i] = -1;
		for(int j = 1; j <= cnt && prime[j] * i <= MAXN - 5; j ++) {
			vis[i * prime[j]] = 1;
			if(i % prime[j] == 0) {
				mu[i * prime[j]] = 0; break;
			}
			mu[i * prime[j]] = -mu[i];
		}
	}
	for(int i = 1; i <= MAXN - 5; i ++) mu[i] += mu[i - 1];
}
int T;
LL Calc(int x, int y) {
	LL res = 0;
	for(int l = 1, r; l <= x && l <= y; l = r + 1) {
		r = min(x / (x / l), y / (y / l));
		res += 1LL * (mu[r] - mu[l - 1]) * p[x / l] * p[y / l];
	}	
	return res;
}
LL calc(int x) {
	LL res = 0;
	for(int l = 1, r; l <= x; l = r + 1) {
		r = x / (x / l); res += (r - l + 1) * 1LL * (x / l);
	}
	return res;
}
int main() {
	primes();
	for(int i = 1; i <= 50000; i ++) p[i] = calc(i);
	for(scanf("%d", &T); T; T --) {
		scanf("%d%d", &n, &m); printf("%lld\n", Calc(n, m));
	}
	return 0;
}