「SDOI 2015」约数个数和「莫比乌斯反演」

题意

设$d(x)$为$x$的约数个数，求${\sum }_{i=1}^n{\sum }_{j=1}^{m}d(ij)$。

题解

首先有个公式：$$d(ij)=\sum _{x|i}\sum _{y|j}[gcd(x,y)=1]$$

证明大致如下：

先考虑这种情况：$i=p^a$,$j=p^b$：

乘积$ij$的约数个数为$a+b+1$个

而${\sum }_{x|i}{\sum }_{y|j}[gcd(x,y)=1]$也为$a+b+1$，因为方案为让$i$取$0$，$j$取$[0,b]$，或$i$取$[0,a]$，$j$取$0$.

然后考虑把每一位合起来，也是一个道理，$i$的次数分别取$0\text{or}[0,a_i]$,$j$的次数分别取$0\text{or}[0,b_i]$

然后开始推导本题解法：

${\sum }_{i=1}^n{\sum }_{j=1}^m{\sum }_{x|i}{\sum }_{y|j}[gcd(x,y)=1]$

先枚举因数：

${\sum }_{x=1}^n{\sum }_{y=1}^m\lfloor \frac{m}{x}\rfloor \lfloor \frac{n}{y}\rfloor [gcd(x,y)=1]$

令$f(u)={\sum }_{x=1}^n{\sum }_{y=1}^m\lfloor \frac{m}{x}\rfloor \lfloor \frac{n}{y}\rfloor [gcd(x,y)=u]$

因此令：

$$g(u)=\sum _{u|d}f(d)=\sum _{x=1}^n\sum _{y=1}^m\lfloor \frac{m}{x}\rfloor \lfloor \frac{n}{y}\rfloor [u|gcd(x,y)]$$

把$u$提出来：

$$g(u)=\sum _{u|d}f(d)=\sum _{x=1}^{\lfloor \frac{n}{u}\rfloor }\sum _{y=1}^{\lfloor \frac{m}{u}\rfloor }\lfloor \frac{m}{xu}\rfloor \lfloor \frac{n}{yu}\rfloor$$

令$sum(n)={\sum }_{i=1}^n\lfloor \frac{n}{i}\rfloor$，则上式：

$$g(u)=sum(\lfloor \frac{n}{u}\rfloor )sum(\lfloor \frac{m}{u}\rfloor )$$

然后莫比乌斯反演：

因为$$g(u)=\sum _{u|d}f(d)$$

所以$$f(u)=\sum _{u|d}\mu (\frac{d}{u})g(d)$$

$$f(u)=\sum _{u|d}\mu (\frac{d}{u})g(d)$$

所以$$f(1)=\sum _{d=1}^n\mu (d)g(d)$$

$$f(1)=\sum _{d=1}^n\mu (d)sum(\lfloor \frac{n}{d}\rfloor )sum(\lfloor \frac{m}{d}\rfloor )$$

于是$O(n\sqrt{n})$预处理$sum$，查询$O(\sqrt{n})$。

#include <cstdio>
#include <ios>

typedef long long LL;

const int N = 5e4 + 10;

LL sum[N], mu[N], mus[N];
int pr[N], tot;
bool tag[N];

LL calc1(int n) {
	LL ans = 0;
	for(int i = 1, j; i <= n; i = j + 1) {
		j = n / (n / i);
		ans += n / i * 1ll * (j - i + 1);
	}
	return ans;
}

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

LL calc2(int n, int m) {
	LL ans = 0;
	for(int i = 1, j; i <= n; i = j + 1) {
		j = std :: min(n / (n / i), m / (m / i));
		ans += (mus[j] - mus[i - 1]) * sum[n / i] * sum[m / i];
	}
	return ans;
}

int main() {
	init(50000);
	int t, n, m;
	for(scanf("%d", &t); t --; ) {
		scanf("%d%d", &n, &m);
		if(n > m) n ^= m ^= n ^= m;
		printf("%lld\n", calc2(n, m));
	} 
	return 0;
}