SPOJ 4491 PGCD - Primes in GCD Table(莫比乌斯反演)

来源:互联网 发布:域名购买后要做什么 编辑:程序博客网 时间:2024/05/16 07:38

Description
给定N, M,求1<=x<=N, 1<=y<=M且gcd(x, y)为质数的(x, y)有多少对
Input
第一行一个整数T 表述数据组数接下来T行,每行两个正整数,表示N, M
(T<=10000,N,M<=10000000)
Output
T行,每行一个整数表示第i组数据的结果
Sample Input
2
10 10
100 100
Sample Output
30
2791
Solution
这里写图片描述
Code

#include<cstdio>#include<iostream>#include<cstring>#include<algorithm>using namespace std;#define maxn 11111111typedef long long ll;bool check[maxn];int prime[maxn],mu[maxn],sum[maxn];void Moblus(int n){    memset(check,0,sizeof(check));    mu[1]=1;    int tot=0;    for(int i=2;i<=n;i++)    {        if(!check[i])        {            prime[tot++]=i;            mu[i]=-1;        }        for(int j=0;j<tot;j++)        {            if(i*prime[j]>n)break;            check[i*prime[j]]=1;            if(i%prime[j]==0)            {                mu[i*prime[j]]=0;                break;            }            else mu[i*prime[j]]=-mu[i];        }    }}//1<=x<=a,a<=y<=b,f(a,b)表示gcd(x,y)=1的对数 ll f(int a,int b){    if(a>b)swap(a,b);    ll ans=0;    for(int i=1,next=0;i<=a;i=next+1)    {        next=min(a/(a/i),b/(b/i));        ans+=1ll*(a/i)*(b/i)*(sum[next]-sum[i-1]);    }    return ans;}int main(){    Moblus(maxn-10);    memset(sum,0,sizeof(sum));    for(int i=0;prime[i]<maxn-10;i++)        for(int j=1;j*prime[i]<maxn-10;j++)            sum[j*prime[i]]+=mu[j];    for(int i=2;i<maxn-10;i++)sum[i]+=sum[i-1];    int T,n,m;    scanf("%d",&T);    while(T--)    {        scanf("%d%d",&n,&m);        ll ans=f(n,m);        printf("%lld\n",ans);    }    return 0;}
0 0