P3327 [SDOI2015]约数个数和

https://www.luogu.com.cn/problem/P3327
题意：设 (d(x)) 为 (x) 的约数个数，求：

[sum_{i=1}^nsum_{j=1}^m d(icdot j) ]

第一个莫比乌斯反演，虽然是看着题解做的，大致明白了一下那个公式该咋用

首先证明一个结论：

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

考虑一个质数 (p)，在 (i,j,k) 中分别是 (p^a,p^b,p^c)
若 (cle a)，则将 (k) 中的 (c) 个 (p) 全部“放入”到 (i) 中
若 (c>a)，则将 (c-a) 个 (p) 放入到 (j) 中
然后同样地推广到其他质数
这样构造，使得对于每一种对 (i,j) 的分配（分配出 (x) 和 (y)），都有唯一对应的 (k)；每一个 (k)，也都只有一种分配方式。那么这就是一个一一映射
那么 (i,j) 不可能同时拥有某个质数 (p)，那么就说明了 (x) 和 (y) 互质，因此得证

更一般的一个推广：

[sigma_0(n_1n_2cdots n_k)=sum_{a_1|n_1}sum_{a_2|n_2}cdots sum_{a_k|n_k}prod_{1le i eq jle k} [gcd(a_i,a_j)=1] ]

---

看这个题，对式子变形

[sum_{i=1}^nsum_{j=1}^m d(icdot j) ]

[sum_{i=1}^nsum_{j=1}^msum_{x|i}^nsum_{y|j}^m [gcd(x,y)=1] ]

变换求和顺序，将 (x,y) 换到前面，那为了保证 (i,j) 分别是它们的倍数，就分别有 (lfloorfrac{n}{x} floor,lfloorfrac{m}{y} floor) 种取值，每种贡献为 (1)

[sum_{x=1}^nsum_{y=1}^m [gcd(x,y)=1] lfloorfrac{n}{x} floorlfloorfrac{m}{y} floor ]

于是开始应用公式进行反演

[ ext{设 } f(N)=sum_{x=1}^nsum_{y=1}^m [gcd(x,y)=N] lfloorfrac{n}{x} floorlfloorfrac{m}{y} floor ]

[g(N)=sum_{N|d}f(d) ]

由于当 (d) 较大时，(f(d)=0)，所以上述定义的 (g(N)) 有意义
对其进行变换，因为它就是 (f) 取遍 (N,2N,3N,cdots) 的所有函数值的和，而这些函数值的贡献都没有交集（要求 (gcd(x,y)) 取不同的数），所以整理后再把 (N) 从 (x,y) 中提取出来：

[egin{aligned} g(N) &=sum_{x=1}^nsum_{y=1}^m [N | gcd(x,y)] lfloorfrac{n}{x} floorlfloorfrac{m}{y} floor\ &=sum_{x=1}^{lfloorfrac{n}{N} floor}sum_{y=1}^{lfloorfrac{m}{N} floor} lfloorfrac{n}{iN} floorlfloorfrac{m}{jN} floor\ &=sum_{x=1}^{lfloorfrac{n}{N} floor}leftlfloorfrac{lfloorfrac{n}{N} floor}{x} ight floor sum_{y=1}^{lfloorfrac{m}{N} floor} leftlfloorfrac{lfloorfrac{m}{N} floor}{y} ight floor\ &=s(lfloorfrac{n}{N} floor)cdot s(lfloorfrac{m}{N} floor)\ end{aligned} ]

借助整除分块预处理 (s)
然后根据反演公式，(f(n)=sum_{n|d} mu(frac{d}{n}) f(d))
那么答案就是 (f(1)=sum_{1|d} mu(frac{d}{1}) f(d)=sum_{i=1}^nmu(i)g(i)=sum_{i=1}^n mu(i) s(lfloorfrac{n}{i} floor) s(lfloorfrac{m}{i} floor))

再预处理一个 (mu) 就行了

#include<cstdio>
#include<algorithm>
#include<iostream>
#include<cmath>
#include<map>
#include<iomanip>
#include<cstring>
#define reg register
#define EN puts("")
inline int read(){
	register int x=0;register int y=1;
	register char c=std::getchar();
	while(c<'0'||c>'9'){if(c=='-') y=0;c=getchar();}
	while(c>='0'&&c<='9'){x=x*10+(c^48);c=getchar();}
	return y?x:-x;
}
#define N 50006
int mu[N],prime[N],notprime[N];
long long s[N];
inline void get_mu(){
	mu[1]=1;
	for(reg int i=2;i<=50000;i++){
		if(!notprime[i]) prime[++prime[0]]=i,mu[i]=-1;
		for(reg int j=1;j<=prime[0]&&i*prime[j]<=50000;j++){
			notprime[i*prime[j]]=1;
			if(i%prime[j]) mu[i*prime[j]]=-mu[i];//mu[i]*mu[prime[j]]=-mu[i]
			else{
				mu[i*prime[j]]=0;break;
			}
		}
	}
	for(reg int x=1;x<=50000;x++){
		mu[x]+=mu[x-1];
		for(reg int i=1,nex;i<=x;i=nex+1){
			nex=x/(x/i);
			s[x]+=(long long)(nex-i+1)*(x/i);
		}
	}
}
int n,m;
int main(){
	get_mu();
	int T=read();while(T--){
		n=read();m=read();
		if(n>m) n^=m,m^=n,n^=m;
		long long ans=0;
		for(reg int i=1,nex;i<=n;i=nex+1){
			nex=std::min(n/(n/i),m/(m/i));
			ans+=(long long)(mu[nex]-mu[i-1])*s[n/i]*s[m/i];
		}
		printf("%lld
",ans);
	}
	return 0;
}