P3327 [SDOI2015] 约数个数和

P3327 [SDOI2015] 约数个数和

大意

给定 \(n, m\)，求 \(\sum_{i=1}^n \sum_{j=1}^m d(ij)\)，其中 \(d(x)\) 表示 \(x\) 的约数个数。

思路

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

将该恒等式代入原式：

\[\text{Ans} = \sum_{i=1}^n \sum_{j=1}^m \sum_{x|i} \sum_{y|j} [\gcd(x, y) = 1] \]

通过更换求和次序，先枚举 \(x\) 和 \(y\)。注意到在 \(1 \sim n\) 中，\(x\) 作为因子的次数为 \(\lfloor \frac{n}{x} \rfloor\)，同理 \(y\) 出现的次数为 \(\lfloor \frac{m}{y} \rfloor\)：

\[\text{Ans} = \sum_{x=1}^n \sum_{y=1}^m \lfloor \frac{n}{x} \rfloor \lfloor \frac{m}{y} \rfloor [\gcd(x, y) = 1] \]

利用莫比乌斯反演的性质 \([\gcd(x, y) = 1] = \sum_{k|\gcd(x, y)} \mu(k)\) 代入：

\[\text{Ans} = \sum_{x=1}^n \sum_{y=1}^m \lfloor \frac{n}{x} \rfloor \lfloor \frac{m}{y} \rfloor \sum_{k|x, k|y} \mu(k) \]

再次变换求和顺序，先枚举 \(k\)：

\[\text{Ans} = \sum_{k=1}^{\min(n, m)} \mu(k) \left( \sum_{x=1}^{\lfloor n/k \rfloor} \lfloor \frac{n}{kx} \rfloor \right) \left( \sum_{y=1}^{\lfloor m/k \rfloor} \lfloor \frac{m}{ky} \rfloor \right) \]

令 \(S(N) = \sum_{i=1}^N \lfloor \frac{N}{i} \rfloor\)，则最终公式为：

\[\text{Ans} = \sum_{k=1}^{\min(n, m)} \mu(k) S(\lfloor \frac{n}{k} \rfloor) S(\lfloor \frac{m}{k} \rfloor) \]

代码

#include<iostream>
using namespace std;

const int MAXN = 50005;
long long sum[MAXN];
int mu[MAXN], pr[MAXN], tot = 0;
long long S[MAXN], d[MAXN];
bool vis[MAXN];

void init(int n){
	mu[1] = 1;
	for(int i = 2;i <= n;i ++){
		if(!vis[i]){
			pr[++ tot] = i;
			mu[i] = -1;
		}
		for(int j = 1;j <= tot && pr[j] * i <= n;j ++){
			vis[pr[j] * i] = 1;
			if(i % pr[j] == 0){
				mu[i * pr[j]] = 0;
				break;
			}
			mu[i * pr[j]] = -mu[i];
		}
	}
	for(int i = 1;i <= n;i ++){
		sum[i] = sum[i - 1] + mu[i];
	}
	for(int i = 1;i <= n;i ++){
		for(int j = i;j <= n;j += i){
			d[j] ++;
		}
	}
	for(int i = 1;i <= n;i ++){
		S[i] = S[i - 1] + d[i];
	}
}

void solve(){
	int n, m;
	cin >> n >> m;
	if(n > m) swap(n, m);
	long long ans = 0;
	for(int l = 1, r;l <= n;l = r + 1){
		r = min(n / (n / l), m / (m / l));
		ans += (long long)(sum[r] - sum[l - 1]) * S[n / l] * S[m / l];
	}
	cout << ans << '\n';
}

int main(){
	
	init(MAXN - 4);
	
	int t; cin >> t;
	while(t --){
		solve();
	}
	return 0;
}