【数论】[SDOI2015]约数个数和

传送门：
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;
}