bzoj3994/ SDOI2015 约数个数和（莫比乌斯函数

题目描述

设d(x)为x的约数个数，给定N、M，求 ${\sum }_{i=1}^n{\sum }_{j=1}^{m}d\left(ij\right)$i=1∑n​j=1∑m​d(ij)

输入格式

输入文件包含多组测试数据。第一行，一个整数T，表示测试数据的组数。接下来的T行，每行两个整数N、M。

输出格式

T行，每行一个整数，表示你所求的答案。

输入输出样例

输入 #1

2
7 4
5 6

输出 #1

110
121
说明/提示
1<=N, M<=50000

1<=T<=50000

分析

首先有一个结论（别看我我不会证
$d\left(ij\right)=\sum _{x\mid i}\sum _{y\mid j}\left[gcd\right(x,y)==1]$d(ij)=x∣i∑​y∣j∑​[gcd(x,y)==1]
于是
$ans={\sum }_{i=1}^n{\sum }_{j=1}^{m}d\left(ij\right)={\sum }_{i=1}^n{\sum }_{j=1}^m\sum _{x\mid i}\sum _{y\mid j}\left[gcd\right(x,y)==1]={\sum }_{i=1}^n{\sum }_{j=1}^m\sum _{x\mid i}\sum _{y\mid j}\sum _{k\mid x,k\mid y}\mu \left(k\right)={\sum }_{k=1}^n\mu \left(k\right){\sum }_{i=1}^{\lfloor {\displaystyle \frac{n}{k}}\rfloor }\lfloor {\displaystyle \frac{n}{kd}}\rfloor {\sum }_{j=1}^{\lfloor {\displaystyle \frac{m}{k}}\rfloor }\lfloor {\displaystyle \frac{m}{kj}}\rfloor$ans=i=1∑n​j=1∑m​d(ij)=i=1∑n​j=1∑m​x∣i∑​y∣j∑​[gcd(x,y)==1]=i=1∑n​j=1∑m​x∣i∑​y∣j∑​k∣x,k∣y∑​μ(k)=k=1∑n​μ(k)i=1∑⌊kn​⌋​⌊kdn​⌋j=1∑⌊km​⌋​⌊kjm​⌋
记 $g\left(n\right)={\sum }_{i=1}^n\lfloor {\displaystyle \frac{n}{i}}\rfloor$g(n)=i=1∑n​⌊in​⌋，于是
$ans={\sum }_{k=1}^n\mu \left(k\right)\ast g\left(\lfloor {\displaystyle \frac{n}{k}}\rfloor \right)\ast g\left(\lfloor {\displaystyle \frac{m}{k}}\rfloor \right)$ans=k=1∑n​μ(k)∗g(⌊kn​⌋)∗g(⌊km​⌋)，除法分块即可。
总复杂度 $O(n\sqrt{n}+T\sqrt{n})$O(nn ​+Tn ​)

代码如下

#include <bits/stdc++.h>
#define N 50005
#define LL long long
using namespace std;
int g[N], mu[N], s[N], x[N], p[N], cnt, maxn = 50000;
LL ans, z = 1;
int main(){
	int i, j, n, m, t, l, r, T;
	mu[1] = 1;
	for(i = 2; i <= maxn; i++){
		if(!x[i]) x[i] = 1, p[++cnt] = i, mu[i] = -1;
		for(j = 1; j <= cnt; j++){
			t = i * p[j];
			if(t > maxn) break;
			x[t] = 1;
			if(i % p[j] == 0) break;
			mu[t] = -mu[i];
		}
	}
	for(i = 1; i <= maxn; i++) s[i] = s[i - 1] + mu[i];
	for(i = 1; i <= maxn; i++){
		for(l = 1; l <= i; l = r + 1){
			r = i / (i / l);
			g[i] += (r - l + 1) * (i / l);
		}
	} 
	scanf("%d", &T);
	while(T--){
		scanf("%d%d", &n, &m);
		if(n > m) swap(n, m);
		ans = 0;
		for(l = 1; l <= n; l = r + 1){
			r = min(n / (n / l), m / (m / l));
			ans += z * (s[r] - s[l - 1]) * g[n / l] * g[m / l];
		}
		printf("%lld\n", ans);
	}
	return 0;
}