BZOJ 2820 YY的GCD （莫比乌斯反演）

题目链接：
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) ∑ n a=1 ∑ m b=1 gcd(a,b)==p(n<m) 

=∑ isprime(p) ∑ ⌊np ⌋ a=1 ∑ ⌊mp ⌋ b=1 gcd(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 .
=∑ n k=1 ∑ isprime(p)∧p|k μ(kp )⌊nk ⌋⌊mk ⌋ 

∑ n k=1 sum(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;}