【莫比乌斯反演】GCD1

还是不要脸的先给链接

hz2016评测《《点击访问

caioj《《点击访问

这题是莫比乌斯反演的模板题

 

  只要让F(t)=满足gcd(x,y)%t==0的数对个数 f(t)满足gcd(x,y)=t的数对个数，则F(t)和f(t)就存在莫比乌斯反演的关系了。显然F(t)=(b/t)*(d/t) 因为如果gcd(x,y)=1,则gcd(x?k,y?k)=k，所以我们把b和d同时除以k， 得到的f(1)再去重就是答案。令lim=min(b/k,d/k),就根据我给出的莫比乌斯反演第二个公式
 
然后就可以来做我们第一题莫比乌斯反演啦。

我这还有一篇莫比乌斯反演的推论，莫比乌斯反演推导《《点击访问
#include<map>
#include<queue>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
#define Maxprime 15485863
#define Maxx 1000000
#define Maxn 1000000
#define mes(x,y) memset(x,y,sizeof(x));
#define mpy(x,y) memcpy(x,y,sizeof(x))
#define INF 2147483647
using namespace std;
int T,a,b,c,d,k,prime[Maxn+1],U[Maxn+1];//U表示系数 
long long ans1,ans2;
bool v[Maxn+1];
void get_Mobius(){
    mes(v,false);
    U[1]=1;
    prime[0]=0;
    for(int i=2;i<=Maxx;i++){
        if(v[i]==false){
            prime[++prime[0]]=i;
            U[i]=-1;//本身是个质数 需要利用G(1)来消去F(1)所以系数为-1 
        }
        for(int j=1;(j<=prime[0])&&(i*prime[j]<=Maxx);j++){
            if(i*prime[j]>Maxx)break;
            v[i*prime[j]]=true;
            if(i%prime[j]==0){
                U[i*prime[j]]=0;//除开质数的两种情况 
                break;
            }
            else{
                U[i*prime[j]]=-U[i];//多了一个数的乘积，与前一个相反 
            }
        }
    }
}
/*
只要让F(t)=满足gcd(x,y)%t==0的数对个数 
f(t)满足gcd(x,y)=t的数对个数，则F(t)和f(t)就存在莫比乌斯反演的关系了。显然F(t)=(b/t)*(d/t) 
因为如果gcd(x,y)=1,则gcd(x?k,y?k)=k，所以我们把b和d同时除以k，
得到的f(1)再去重就是答案。令lim=min(b/k,d/k),就根据我给出的莫比乌斯反演第二个公式 
*/
int main(){
    get_Mobius();
    scanf("%d",&T);
    while(T--){
        scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
        if(k==0){
            printf("0\n");
            continue;
        }
        b/=k;a/=k;
        d/=k;c/=k;
        if(b>d)swap(b,d);//只要for一遍i~最小的就好 
        ans1=ans2=0;
        for(int i=1;i<=b;i++)ans1+=(long long)U[i]*(b/i)*(d/i);//带入 G(n) 的公式有 F(1)=sigma(U[i]*(a/i)*(b/i))(1<=i<=N)
        for(int i=1;i<=b;i++)ans2+=(long long)U[i]*(b/i)*(b/i);//x,y相同数的状况 
        ans1-=ans2/2;
        printf("%lld\n",ans1);
    }
    return 0;
}
 

查看原文：http://hz2016.tk/blog/?p=37