P3327 [SDOI2015]约数个数和

Label

莫比乌斯反演

Description

设$d(x)$为$x$的约数个数，给定$T$组$n,m(1\le T,n,m\le 5\times 10^4)$，对于每组$n,m$，求

​ $$\sum _{i=1}^n\sum _{j=1}^{m}d(ij)$$

Solution

为了将$d$变化为可利用反演结论化简的式子，我们首先引入一个结论：

​ $$d(ij)=\sum _{x|i}\sum _{y|j}[gcd(i,j)=1]$$ $(1)$

该结论中出现了$[gcd(i,j)=1]$的形式，我们考虑进一步变换：

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

$$=\sum _{x|i}\sum _{y|j}\sum _{p|gcd(x,y)}\mu (p)$$

$$=\sum _{p|i\wedge p|j}\mu (p)\sum _{x|i}\sum _{y|j}[p|gcd(x,y)]$$

$$=\sum _{p|i\wedge p|j}\mu (p)\sum _{\frac{x}{p}|\frac{i}{p}}\sum _{\frac{y}{p}|\frac{j}{p}}[1|gcd(\frac{x}{p},\frac{y}{p})]$$ $(2)$

$$=\sum _{p|i\wedge p|j}\mu (p)\sum _{\frac{x}{p}|\frac{i}{p}}\sum _{\frac{y}{p}|\frac{j}{p}}1$$

$$=\sum _{p|i\wedge p|j}\mu (p)\sum _{x|\frac{i}{p}}1\sum _{y|\frac{j}{p}}1$$

$$=\sum _{p|i\wedge p|j}\mu (p)d(\frac{i}{p})d(\frac{j}{p})$$

$$\therefore \sum _{i=1}^n\sum _{j=1}^{m}d(ij)$$

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

$$=\sum _{i=1}^n\sum _{j=1}^m=\sum _{p|i\wedge p|j}\mu (p)d(\frac{i}{p})d(\frac{j}{p})$$

$$=\sum _{p=1}^{min(n,m)}\mu (p)\sum _{i=1}^n\sum _{j=1}^m[p|gcd(i,j)]d(\frac{i}{p})d(\frac{j}{p})$$ $(3)$

$$=\sum _{p=1}^{min(n,m)}\mu (p)\sum _{i=1}^n\sum _{j=1}^m[1|gcd(\frac{i}{p},\frac{j}{p})]d(\frac{i}{p})d(\frac{j}{p})$$

$$=\sum _{p=1}^{min(n,m)}\mu (p)\sum _{i=1}^{\lfloor \frac{n}{p}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{p}\rfloor }d(i)d(j)$$

$$=\sum _{p=1}^{min(n,m)}\mu (p)S(\lfloor \frac{n}{p}\rfloor )S(\lfloor \frac{m}{p}\rfloor )$$

（其中$S(x)={\sum }_{k=1}^x(d(k))$）

算法时间复杂度$O(n+T\sqrt{n})$

总结

注释：

（1）积累结论：$d(ij)={\sum }_{x|i}{\sum }_{y|j}[gcd(i,j)=1]$。

（2）与将$gcd(x,y)=d$转化为$gcd(\frac{x}{d},\frac{y}{d})$的套路类似，遇到形似$[x|y]$时，可考虑将其转化为$[1|\frac{y}{x}]$，当$x$整除$y$时此真值式恒为1。

（3）求和符号交换时，需考虑合适的将被交换的求和符号对应的范围转化的方式：有时考虑所有合法值的取值范围，有时转化为真值式$[...]$表示更加简洁。

Code

#include<cstdio>
#include<iostream>
#define ri register int
#define ll long long
using namespace std;

const int MAXN=5e4;
int T,N,M,n,cnt;
ll ans,mu[MAXN+20],prime[MAXN+20],sum[MAXN+20],d[MAXN+20];
bool notprime[MAXN+20];

void Mobius() 
{
	mu[1]=1,notprime[1]=true;
	for(ri i=2;i<=MAXN;++i)
	{
		if(!notprime[i]) prime[++cnt]=i,mu[i]=-1;
		for(ri j=1;j<=cnt&&i*prime[j]<=MAXN;++j)
		{
			notprime[i*prime[j]]=true;
			if(i%prime[j]==0) break;
			else mu[i*prime[j]]=-mu[i];
		}
	}
	for(ri i=1;i<=MAXN;++i) sum[i]=sum[i-1]+mu[i];
}

void GetD()
{
	for(ri i=1;i<=MAXN;++i)
		for(ri j=i;j<=MAXN;j+=i) ++d[j];
	for(ri i=1;i<=MAXN;++i) d[i]=d[i-1]+d[i];
}

int main()
{
	std::ios::sync_with_stdio(false);
	cin>>T;
	Mobius(); GetD();
	for(ri op=1;op<=T;++op)
	{
		cin>>N>>M;
		n=min(N,M),ans=0LL;
		for(ri l=1,r;l<=n;l=r+1)
		{
			r=min(N/(N/l),M/(M/l));
			ans+=(sum[r]-sum[l-1])*d[N/l]*d[M/l];
		}
		cout<<ans<<'\n';
	}
	return 0;
}