传送门：
https://www.luogu.com.cn/problem/P3327
https://www.acwing.com/problem/content/1360/

莫比乌斯反演 + 整除分块

分析

首先，我们给出一个结论：

\[d(ij) = \sum_{x|i} \sum_{y|j} [(x, y) = 1] \]

证明：

设 \(i\) 的分解式为 \(i = \prod p_k^{\alpha_k}\) ，类似地，\(j = \prod p_k^{\beta_k}\)

对于质数 \(p_k\) ，\(ij\) 对应的因数个数为 \(\alpha_k + \beta_k + 1\)

而对于右式，\(p_k\) 如果想要产生贡献，那么只可能是：\(x, y\) 不同时存在因数 \(p_k\)，注意到对应的贡献数即为 \(\alpha_k + \beta_k + 1\)

因此，由乘法原理，所求证式成立。

为了防止混淆，我们将题面中的 \(n, m\) 记为 \(N, M\)

由上述结论，题目的式子转化为：

\[\sum_{i=1}^N \sum_{j=1}^M \sum_{x|i} \sum_{y|j} [(x, y) = 1] \]

莫比乌斯反演：若 \(F(n) = \sum_{n|d}f(d)\)，则有 \(f(n) = \sum_{n|d} \mu(\frac{d}{n})F(d)\)

于是我们设 \(f(n) = \sum_{i=1}^N \sum_{j=1}^M \sum_{x|i} \sum_{y|j} [(x, y) = n]\) ，对应的 \(F(n) = \sum_{i=1}^N \sum_{j=1}^M \sum_{x|i} \sum_{y|j} [n | (x, y)]\)

下面对 \(F(n)\) 进行变换：

\(F(n) \\= \sum_{i=1}^N \sum_{j=1}^M \sum_{x|i} \sum_{y|j} [n | (x, y)] \\ = \sum_{x=1}^N \sum_{y=1}^M \lfloor \frac{N}{x} \rfloor \lfloor \frac{M}{y} \rfloor [n | (x, y)] \\\)

设 \(x' = \lfloor \frac{x}{n} \rfloor ~, y'=\lfloor \frac{y}{n} \rfloor ~ , N'=\lfloor \frac{N}{n} \rfloor ~ , M'=\lfloor \frac{M}{n} \rfloor\) ，进而有

\(F(n) \\= \sum_{x=1}^{N'} \sum_{y=1}^{M'} \lfloor \frac{N'}{x'} \rfloor \lfloor \frac{M'}{y'} \rfloor \\ = \sum_{x=1}^{N'} \lfloor \frac{N'}{x'} \rfloor \sum_{y=1}^{M'} \lfloor \frac{M'}{y'} \rfloor\)

记 \(h(x) = \sum_{i=1}^x \lfloor \frac{x}{i} \rfloor\) ，\(h(x)\) 可用整除分块处理。

由莫比乌斯反演，\(f(n) = \sum_{n|d} \mu(\frac{d}{n})F(d)\) ，我们只需求 \(f(1)\)

下面对 \(f(1)\) 进行变换：

\[f(1) = \sum_{d} \mu({d})F(d) = \sum_{d} \mu({d})h(\lfloor \frac{N}{d} \rfloor)h(\lfloor \frac{M}{d} \rfloor) \]

由整除分块知，\(h(\lfloor \frac{N}{d} \rfloor)h(\lfloor \frac{M}{d} \rfloor)\) 取值最多为 \(\sqrt{N} + \sqrt{M}\) 块，而对于每一块，可以用前缀和 \(O(1)\) 处理出对应的 \(\mu\) 值，因此整体的复杂度为 \(O(T(\sqrt{N} + \sqrt{M}))\)

细节见代码：

#include<bits/stdc++.h>
using namespace std;

const int S=50050;

int primes[S], cnt;
bool vis[S];
int mu[S], sum[S];
int h[S];

int g(int b, int l){
	return b/(b/l);
}

void init(){
	// init sum[]
	mu[1]=1;
	for(int i=2; i<S; i++){
		if(!vis[i]) primes[cnt++]=i, mu[i]=-1;
		for(int j=0; i*primes[j]<S; j++){
			vis[i*primes[j]]=true;
			if(i%primes[j]==0) break;
			mu[i*primes[j]]=-mu[i];
		}
	}
	for(int i=1; i<S; i++) sum[i]=sum[i-1]+mu[i];
	
	// init h[]
	for(int x=1; x<S; x++){
		int &v=h[x];
		for(int l=1, r; l<=x; l=r+1){
			r=min(x, g(x, l));
			v+=x/l*(r-l+1);
		}
	}
}

int solve(int N, int M){
	int res=0;
	
	int n=min(N, M);
	for(int l=1, r; l<=n; l=r+1){
		r=min(n, min(g(N, l), g(M, l)));
		res+=(sum[r]-sum[l-1])*h[N/l]*h[M/l];
	}
	return res;
}

int main(){
	int T; cin>>T;
	init();
	while(T--){
		int N, M; cin>>N>>M;
		cout<<solve(N, M)<<endl;
	}
	return 0;
}