【数论】[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^{eta_k})

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

而对于右式，(p_k) 如果想要产生贡献，那么只可能是：(x, y) 不同时存在因数 (p_k)，注意到对应的贡献数即为 (alpha_k + eta_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} floor lfloor frac{M}{y} floor [n | (x, y)] \)

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

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

记 (h(x) = sum_{i=1}^x lfloor frac{x}{i} floor) ，(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} floor)h(lfloor frac{M}{d} floor) ]

由整除分块知，(h(lfloor frac{N}{d} floor)h(lfloor frac{M}{d} floor)) 取值最多为 (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;
}