P3327 [SDOI2015]约数个数和

洛谷题目链接

挺简单的数论题

主要要知道一个式子：$$d(i\times j)=\sum\limits_{x|i} \sum\limits_{y|j} [gcd(x,y)=1]$$

具体证明其实挺简单的，这里给出一个大佬的博客

那么我们需要求的式子就是：$$\sum\limits_{i=1}^N \sum\limits_{j=1}^M \sum\limits_{x|i} \sum\limits_{y|j} [gcd(x,y)=1]$$

分组一下：$$\sum\limits_{i=1}^N [x|i] \sum\limits_{j=1}^M [y|j] [gcd(x,y)=1]$$

按照套路改为枚举\(x,y\)：$$\sum\limits_{x=1}^N \left\lfloor \frac{N}{x} \right\rfloor \sum\limits_{y=1}^M \left\lfloor \frac{M}{y} \right\rfloor [gcd(x,y)=1]$$

套路的反演也行，不过这里用的是莫比乌斯函数的性质：$$\sum\limits_{x=1}^N \left\lfloor \frac{N}{x} \right\rfloor \sum\limits_{y=1}^M \left\lfloor \frac{M}{y} \right\rfloor \sum\limits_{d|gcd(x,y)} \mu (d)$$

改为枚举\(d\)并提前：$$\sum\limits_{d=1}^N \sum\limits_{x=1}^N \sum\limits_{y=1}^M [d|gcd(x,y)] \mu (d) \left\lfloor \frac{N}{x} \right\rfloor \left\lfloor \frac{M}{y} \right\rfloor$$

\[\sum\limits_{d=1}^N \sum\limits_{x=1}^{\left\lfloor \frac{N}{d} \right\rfloor} \sum\limits_{y=1}^{\left\lfloor \frac{M}{d} \right\rfloor} \mu (d) \left\lfloor \frac{N}{dx} \right\rfloor \left\lfloor \frac{M}{dy} \right\rfloor \]

分组：$$\sum\limits_{d=1}^N \mu (d) \sum\limits_{x=1}^{\left\lfloor \frac{N}{d} \right\rfloor} \left\lfloor \frac{N}{dx} \right\rfloor \sum\limits_{y=1}^{\left\lfloor \frac{M}{d} \right\rfloor} \left\lfloor \frac{M}{dy} \right\rfloor$$

我们只要求出\(\sum\limits_{i=1}^n \left\lfloor \frac{n}{i} \right\rfloor\)就行

接下来是美滋滋的代码时间~~~

#include<iostream>
#include<cstdio>
#include<cstring>
#define N 50007
#define int long long
using namespace std;
int T,n,m,cnt;
int mu[N],sum[N],num[N],prime[N];
bool isp[N];
void Get()
{
	mu[1]=1;
	for(int i=2;i<=N-7;++i)
	{
		if(!isp[i])
		{
			prime[++cnt]=i;
			mu[i]=-1;
		}
		for(int j=1;j<=cnt&&(i*prime[j])<=N-7;++j)
		{
			isp[i*prime[j]]=1;
			if(i%prime[j]==0)
				break;
			else
				mu[i*prime[j]]=-mu[i];
		}
	}
	for(int i=1;i<=N-7;++i)
		sum[i]=sum[i-1]+mu[i];
	for(int i=1;i<=N-7;++i)
	{
		int res=0;
		for(int l=1,r;l<=i;l=r+1)
		{
			r=i/(i/l);
			res+=(r-l+1)*(i/l);
		}
		num[i]=res;
	}
}
signed main()
{
	Get();
	scanf("%lld",&T);
	while(T--)
	{
		scanf("%lld%lld",&n,&m);
		int ans=0;
		if(n>m)
			swap(n,m);
		for(int l=1,r;l<=n;l=r+1)
		{
			r=min(n/(n/l),m/(m/l));
			ans+=(sum[r]-sum[l-1])*num[n/l]*num[m/l];
		}
		printf("%lld\n",ans);
	}
	return 0;
}