【BZOJ 2820】YY的GCD

【题目】

传送门

题目描述：

神犇 YY 虐完数论后给傻$\times$ kAc 出了一题：

给定 $n,m$，求 $1\le x\le n$，$1\le y\le m$ 且 $\gcd (x,y)$ 为质数的 $(x,y)$ 有多少对。

kAc 这种傻$\times$必然不会了，于是向你来请教$\dots \dots$

输入格式：

第一行一个整数 $T$，表示数据组数。

接下来 $T$ 行，每行两个正整数，表示 $n,m$。

输出格式：

$T$ 行，每行一个整数表示第 $i$ 组数据的结果

样例数据：

输入
2
10 10
100 100

输出
30
2791

提示：

对于 $100\%$ 的数据，$T\le 10000$，$n,m\le 10000000$

【分析】

题目要求的是这个东西：

$$\sum _{i=1}^n\sum _{j=1}^m[\gcd (i,j)\mathrm{i}\mathrm{s}\mathrm{a}\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{m}\mathrm{e}]$$

如果用 $\mathbb{P}$ 表示素数集的话，上式就是

$$\sum _{p\in \mathbb{P}}\sum _{i=1}^n\sum _{j=1}^m[\gcd (i,j)=p]$$

我们定义

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

$$F(x)=\sum _{x|d}f(d)$$

显然，我们可以得到

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

用莫比乌斯反演得到

$$f(x)=\sum _{x|d}\mu (\frac{d}{x})F(d)$$

那么，最后的答案就是

$$\sum _{p\in \mathbb{P}}f(p)=\sum _{p\in \mathbb{P}}\sum _{p|d}\mu (\frac{d}{p})\lfloor \frac{n}{d}\rfloor \lfloor \frac{m}{d}\rfloor$$

我们把 $d$ 改为枚举 $p$ 的几倍，即

$$\sum _{p\in \mathbb{P}}\sum _{d=1}^{\lfloor \frac{min(n,m)}{p}\rfloor }\mu (d)\lfloor \frac{n}{d\cdot p}\rfloor \lfloor \frac{m}{d\cdot p}\rfloor$$

令 $D=d\cdot p$，交换一下枚举顺序，即

$$\sum _{D=1}^{min(n,m)}\sum _{p|D,p\in \mathbb{P}}\mu (\frac{D}{p})\lfloor \frac{n}{D}\rfloor \lfloor \frac{m}{D}\rfloor$$

继续化简，得到

$$\sum _{D=1}^{min(n,m)}\lfloor \frac{n}{D}\rfloor \lfloor \frac{m}{D}\rfloor (\sum _{p|D,p\in \mathbb{P}}\mu (\frac{D}{p}))$$

令 $g(D)=\sum _{p|D,p\in \mathbb{P}}\mu (\frac{D}{p})$。

那么，我们可以先线性筛出 $\mu$ 的函数值，然后枚举每个质数以及每个质数的倍数把 $g$ 暴力更新出来。

剩下的就可以整除分块做了。

【代码】

#include<cstdio>
#include<cstring>
#include<algorithm>
#define N 10000005
#define ll long long
using namespace std;
bool mark[N];
int prime[N],mu[N],g[N];
void linear_sieves()
{
	int i,j,sum=0;
	memset(mark,true,sizeof(mark));
	mark[0]=mark[1]=false,mu[1]=1;
	for(i=2;i<N;++i)
	{
		if(mark[i])  prime[++sum]=i,mu[i]=-1;
		for(j=1;j<=sum&&i*prime[j]<N;++j)
		{
			mark[i*prime[j]]=false;
			if(i%prime[j])  mu[i*prime[j]]=-mu[i];
			else  {mu[i*prime[j]]=0;break;}
		}
	}
	for(i=1;i<=sum;++i)
	  for(j=1;j*prime[i]<N;++j)
	    g[j*prime[i]]+=mu[j];
	for(i=1;i<N;++i)  g[i]+=g[i-1];
}
ll solve(int n,int m)
{
	int i,j;ll ans=0;
	if(n>m)  swap(n,m);
	for(i=1;i<=n;i=j+1)
	{
		j=min(n/(n/i),m/(m/i));
		ans+=1ll*(g[j]-g[i-1])*(n/i)*(m/i);
	}
	return ans;
}
int main()
{
	int n,m,T;
	scanf("%d",&T);
	linear_sieves();
	while(T--)
	{
		scanf("%d%d",&n,&m);
		printf("%lld\n",solve(n,m));
	}
	return 0;
}