P3327 [SDOI2015]约数个数和

P3327 [SDOI2015]约数个数和

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

则我们原式化为

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

莫比乌斯反演：

(f(d)=sum_{i=1}^{n}sum_{j=1}^{m}[gcd(i,j)=d])

(F(n)=sum_{n|d}f(d))

(Rightarrow f(n)=sum_{n|d}mu(lfloorfrac{d}{n} floor)F(d))

根据(mu)函数性质

(Ans=sum_{i=1}^{n}sum_{j=1}^{m}sum_{x|i}sum_{y|j}sum_{d|gcd(x,y)}mu(d))

推(du)公(liu)式的地方来了：

(Ans=sum_{i=1}^{n}sum_{j=1}^{m}sum_{x|i}sum_{y|j}sum_{d=1}^{min(n,m)}mu(d)*[d|gcd(x,y)])

(Ans=sum_{d=1}^{min(n,m)}mu(d)sum_{i=1}^{n}sum_{j=1}^{m}sum_{x|i}sum_{y|j}[d|gcd(x,y)])

(Ans=sum_{d=1}^{min(n,m)}mu(d)sum_{x=1}^{n}sum_{y=1}^{m}[d|gcd(x,y)]lfloorfrac{n}{x} floorlfloorfrac{m}{y} floor)

(Ans=sum_{d=1}^{min(n,m)}mu(d)sum_{x=1}^{lfloorfrac{n}{d} floor}sum_{y=1}^{{lfloorfrac{m}{d} floor}}lfloorfrac{n}{dx} floorlfloorfrac{m}{dy} floor)

(Ans=sum_{d=1}^{min(n,m)}mu(d)(sum_{x=1}^{lfloorfrac{n}{d} floor}lfloorfrac{n}{dx} floor)(sum_{y=1}^{{lfloorfrac{m}{d} floor}}lfloorfrac{m}{dy} floor))

#include<cstring>
#include<iostream>
#include<algorithm>
#include<string>
#include<cstdio>
using namespace std;
typedef long long LL;
const LL maxn=50000+10;
inline int Read(){
	LL x=0,f=1; char c=getchar();
	while(c<'0'||c>'9'){
		if(c=='-') f=-1; c=getchar();
	}
	while(c>='0'&&c<='9'){
		x=(x<<3)+(x<<1)+c-'0',c=getchar();
	}
	return x*f;
}
int T;
int prime[maxn],mu[maxn],sum[maxn];
LL a[maxn];
bool visit[maxn];
inline void F_phi(LL N){
	mu[1]=1;
    int tot=0;
	for(int i=2;i<=N;++i){
		if(!visit[i]){
			mu[i]=-1,
			prime[++tot]=i;
		}
		for(int j=1;j<=tot&&i*prime[j]<=N;++j){
			visit[i*prime[j]]=true;
		    if(i%prime[j]==0)
		        break;
		    else
		        mu[i*prime[j]]=-mu[i];
		}
	}
	for(int i=1;i<=N;++i)
	    sum[i]=sum[i-1]+1ll*mu[i];
	for(int i=1;i<=N;++i){
		LL num(0);
		for(int l=1,r;l<=i;l=r+1){
			r=i/(i/l);
			num=num+1ll*(r-l+1)*1ll*(i/l);
		}
		a[i]=num;
	}
}
int main(){
	T=Read();
	F_phi(50000);
	while(T--){
		int n=Read(),m=Read();
		if(n>m)
		    swap(n,m);
		LL ans(0);
		int N=n;
		for(int l=1,r;l<=N;l=r+1){
			r=min(n/(n/l),m/(m/l));
			ans=ans+1ll*(sum[r]-sum[l-1])*(1ll*a[n/l])*(1ll*a[m/l]);
		}
		printf("%lld
",ans);
	}
	return 0;
}/*
2
7 4
5 6

110
121
*/