洛谷P2257 YY的gcd 莫比乌斯反演+整除分块

题目链接：
洛谷

莫比乌斯反演+整除分块$qwq$
如果不会莫比乌斯反演，珂以看我的博客qwq
做这道题之前珂以先看HDU1695和我的题解qwq

先把答案表示出来：
$ans=\Sigma {}_{i=1}^N\Sigma {}_{j=1}^M[gcd(i,j)为质数]$
发现“为质数”不太好转化，珂以考虑把$gcd(i,j)$提前枚举：
$=\Sigma {}_{p为质数}\Sigma {}_{i=1}^N\Sigma {}_{j=1}^M[gcd(i,j)==p]$
同$hdu1695$的推法，这里珂以把$i,j$同时除以$p$：
$=\Sigma {}_{p为质数}\Sigma {}_{i=1}^{N/p}\Sigma {}_{j=1}^{M/p}[gcd(i,j)==1]$
令$f(n)$表示$gcd(i,j)==n$的个数，$F(n)$表示$n|gcd(i,j)$的个数。
则$F(n)=\Sigma {}_{n|d}f(d)$
做一下莫比乌斯反演得$f(n)=\Sigma {}_{n|d}\mu (\frac{d}{n})F(d)=\Sigma {}_{n|d}\mu (\frac{d}{n})\lfloor \frac{N}{d}\rfloor \lfloor \frac{M}{d}\rfloor$
所以
$ans=\Sigma {}_{p为质数}\Sigma {}_{p|d}\mu (\frac{d}{p})F(d)$
令$\frac{d}{p}=t$，把式子改为枚举$t$：
$ans=\Sigma {}_{p为质数}\Sigma {}_{t=1}^{min(n,m)/p}\mu (t)F(tp)=\Sigma {}_{p为质数}\Sigma {}_{t=1}^{min(n,m)/p}\mu (t)\lfloor \frac{N}{tp}\rfloor \lfloor \frac{M}{tp}\rfloor$
令$T=tp$，则
$ans=\Sigma \Sigma \mu (\frac{T}{t})\lfloor \frac{N}{T}\rfloor \lfloor \frac{M}{T}\rfloor$
$=\Sigma {}_{T=1}^{min(n,m)}\lfloor \frac{N}{T}\rfloor \lfloor \frac{M}{T}\rfloor \Sigma {}_{t|T}\mu (\frac{T}{t})$
然后整除分块+前缀和就珂以了qwq

上代码：

#include<stdio.h>
#include<algorithm>
#define re register int
using namespace std;
typedef long long ll;
int read() {
	re x=0,f=1;
	char ch=getchar();
	while(ch<'0' || ch>'9') {
		if(ch=='-')	f=-1;
		ch=getchar();
	}
	while(ch>='0' && ch<='9') {
		x=10*x+ch-'0';
		ch=getchar();
	}
	return x*f;
}
const int Size=10000005;
ll n,m,sum[Size],f[Size];
int tot,p[Size],prime[Size/10]; 
bool vis[Size];
void getp(int maxn) {
	p[1]=1;
	for(re i=2; i<=maxn; i++) {
		if(!vis[i]) {
			p[i]=-1;
			prime[++tot]=i;
		}
		for(re j=1; j<=tot && i*prime[j]<=maxn; j++) {
			vis[i*prime[j]]=true;
			if(i%prime[j]==0)	break;
			p[i*prime[j]]=-p[i];
		}
	}
	for(re i=1; i<=tot; i++) {
		for(re j=1; j*prime[i]<=maxn; j++) {
			f[j*prime[i]]+=p[j];
		}
	}
	for(re i=1; i<=maxn; i++) {
//		sum[i]+=sum[i-1];
		sum[i]=sum[i-1]+f[i];
	}
}
int main() {
	getp(10000000);
	int T=read();
	while(T--) {
		n=read();
		m=read();
		int maxn=min(n,m),lst;
		ll ans=0;
		for(re i=1; i<=maxn; i=lst+1) {
			lst=min(n/(n/i),m/(m/i));
			ans+=(ll)(n/i)*(m/i)*(sum[lst]-sum[i-1]);
		}
		printf("%lld\n",ans);
	}
	return 0;
}