BZOJ 2820 YY问题中，如何运用莫比乌斯反演求解GCD？

本文共计588个文字，预计阅读时间需要3分钟。

题目链接：+BZOJ+2820+限制题+描述：神皇YY无聊完数论后给kAc出了这么一道题。给定N,M，1=x=N，1=y=M且gcd(x,y)为质数的(x,y)有多少对？kAc这种题当然不会了，于是向你请教……多对“

题目链接：
​​​BZOJ 2820 权限题​​

Description
神犇YY虐完数论后给傻×kAc出了一题。给定N,M,求1<=x<=N,1<=y<=M且gcd(x,y)为质数的(x,y)有多少对？kAc这种傻×必然不会了，于是向你来请教……

多组输入
Input
第一行一个整数T 表述数据组数
接下来T行，每行两个正整数，表示N, M
Output
T行，每行一个整数表示第i组数据的结果

Sample Input
2
10 10
100 100

Sample Output
30
2791

HINT
T = 10000
N, M <= 10000000

题解：
莫比乌斯反演。

∑isprime(p)∑na=1∑mb=1gcd(a,b)==p(n<m)

=∑isprime(p)∑⌊np⌋a=1∑⌊mp⌋b=1gcd(a,b)==1

=∑isprime(p)∑⌊np⌋a=1∑⌊mp⌋b=1∑d|gcd(a,b)μ(d)

=∑isprime(p)∑⌊np⌋a=1∑⌊mp⌋b=1∑d|a∧d|bμ(d)

=∑isprime(p)∑⌊np⌋d=1μ(d)⌊npd⌋⌊mpd⌋

设 pd=k.
=∑nk=1∑isprime(p)∧p|kμ(kp)⌊nk⌋⌊mk⌋

∑nk=1sum(k)⌊nk⌋⌊mk⌋

代码：

#include<cstdio>
#include<iostream>
using namespace std;
typedef long long ll;
const int N=1e7+5;
int T,n,m,tot,mu[N],g[N],sum[N],prime[N/3];
bool check[N];
void mobius()
{
mu[1]=1;n=1e7+2;
for(int i=2;i<=n;i++)
{
if(!check[i]){
prime[++tot]=i;
mu[i]=-1;
}
for(int j=1;j<=tot&&i*prime[j]<=n;j++)
{
check[i*prime[j]]=1;
if(!(i%prime[j]))
{
mu[i*prime[j]]=0;
break;
}
else
{
mu[i*prime[j]]=-mu[i];
}
}
}
for(int i=1;i<tot;i++)
{
for(int j=1;j*prime[i]<=n;j++)
{
sum[j*prime[i]] += mu[j];
}
}
for(int i=1;i<=n;i++) sum[i] += sum[i-1];
}
ll calc(int n,int m)
{
if(n>m) swap(n,m);
ll ans=0;int pos=0;
for(int i=1;i<=n;i=pos+1)
{
pos=min(n/(n/i),m/(m/i));
ans+=(ll)(n/i)*(m/i)*(sum[pos]-sum[i-1]);
}
return ans;
}
int main()
{
mobius();
scanf("%d",&T);
while(T--){
scanf("%d%d",&n,&m);
printf("%lld\n",calc(n,m));
}
return 0;
}