[题解] 约数个数和

[题解] 约数个数和 数论分块+莫比乌斯反演

题目链接
首先我们需要知道一个结论:$$d(ij)=\sum _{x|i}\sum _{y|j}\left[\gcd (x,y)==1\right]$$
感性理解就是对约数取了个并集。

对于这个结论，我们还有其扩展形式。
$$d(ijk)=\sum _{x|i}\sum _{y|j}\sum _{z|k}[\gcd (x,y)=1][\gcd (y,z)=1][\gcd (x,z)=1]$$

我们对这个式子进行变形:
$$\begin{array}{rl}d(ij)=& \sum _{x|i}\sum _{y|j}\left[\gcd (x,y)==1\right]\\ =& \sum _{x|i}\sum _{y|j}\sum _{k|\gcd (x,y)}\mu (k)\\ =& \sum _{k|i,k|j}\mu (k)\sum _{x|i}\sum _{y|j}1\end{array}$$
有的同学可能会有疑问，为什么这里不让$k$从$1$到$\min (i,j)$开始枚举呢?因为这里的k没有遍历$1\to \min (i,j)$。因为总是存在$k<\min (i,j)$，使得$k\nmid i$或$k\nmid j$。
我们继续化简:
$$\begin{array}{rl}d(ij)=& \sum _{k|i,k|j}\mu (k)\left(\sum _{x|\frac{i}{k}}1\right)\left(\sum _{y|\frac{j}{k}}1\right)\\ =& \sum _{k|i,k|j}\mu (k)d\left(\frac{i}{k}\right)d\left(\frac{j}{k}\right)\end{array}$$
我们把这个式子代入原式:

$$\begin{array}{rl}\sum _{i=1}^n\sum _{j=1}^{m}d(ij)=& \sum _{i=1}^n\sum _{j=1}^m\sum _{k|i,k|j}\mu (k)d\left(\frac{i}{k}\right)d\left(\frac{j}{k}\right)\\ =& \sum _{k=1}^{\min (n,m)}\mu (k)\sum _{k|i}^{i\le n}d\left(\frac{i}{k}\right)\sum _{k|j}^{j\le m}d\left(\frac{j}{k}\right)\\ =& \sum _{k=1}^{\min (n,m)}\mu (k)\left(\sum _{i=1}^{\frac{n}{k}}d(i)\right)\left(\sum _{j=1}^{\frac{m}{k}}d(j)\right)\end{array}$$
我们设$$s(n)=\sum _{i=1}^{n}d(i)$$
即$s$是$d$的前缀和，那么原式可以化为:
$$\sum _{k=1}^{\min (n,m)}\mu (k)s\left(\lfloor \frac{n}{k}\rfloor \right)s\left(\lfloor \frac{m}{k}\rfloor \right)$$
通过交换求和符号，我们减少了某些不必要的重复枚举，降低了时间复杂度。
*注意 对于函数$s$，我们仍然可以对其用数论分块。运用时只需要把它看成一种映射即可。
前面的$\mu$可以预处理出前缀和，后面应用二维数论分块，总的时间复杂度$O\left(n+T\sqrt{n}\right)$。
这里用线性筛维护$d$的时候，我们需要记录每个正整数$i$的最小素因子的次数$a[i]$，便于更新$d$.
维护一个关于数的最小素因子的数组是常用手法。

#include<bits/stdc++.h>
#define ll long long
#define ull unsigned long long
#define pii pair<int,int>
using namespace std;

const double eps = 1e-10;
const double pi = acos(-1.0);
const int maxn = 5e4 + 10;

int T;
int p[maxn],cnt;
bool isp[maxn];
ll mu[maxn],d[maxn];
ll a[maxn];

void getp(){
	isp[1] = 1;
	mu[1] = d[1] = 1;
	a[1] = 0;
	for(int i = 2; i < maxn; i++){
		if(!isp[i]){
			p[++cnt] = i;
			mu[i] = -1;
			d[i] = 2;
			a[i] = 1;
		}
		for(int j = 1; j <= cnt && i * p[j] < maxn; j++){
			isp[i*p[j]] = 1;
			if(i %p[j] == 0){
				mu[i*p[j]] = 0;
				a[i*p[j]] = a[i]+1;
				d[i*p[j]] = d[i]/(a[i]+1)*(a[i]+2);
				break;
			}
			else{
				mu[i*p[j]] = -mu[i];
				a[i*p[j]] = 1;
				d[i*p[j]] = 2*d[i];
			}
		}
	}
	for(int i = 1; i < maxn; i++){
		mu[i] = mu[i-1] + mu[i];
		d[i] = d[i-1] + d[i];
	}
}

void solve(){
	getp();
	scanf("%d",&T);
	while(T--){
		int n,m;
		scanf("%d%d",&n,&m);
		int t = min(n,m);
		int l = 1, r = t + 1;
		ll ans = 0;
		while(l <= t){
			r = min(n/(n/l),m/(m/l));
			ans += d[n/l]*d[m/l]*(mu[r] - mu[l-1]);
			l = r+1;
		}
		printf("%lld\n",ans);
	}
}

int main()
{
	solve();
	return 0;
}