hdu 1695 GCD 莫比乌斯反演入门

GCD

题意:输入5个数a,b,c,d,k;(a = c = 1, 0 < b,d,k <= 100000);问有多少对a <= p <= b, c <= q <= d使得gcd(p,q) = k;

注:对于(p,q)和(q,p)只算一次；

思路:由于遍历朴素求两个数的gcd的时间复杂度为O(n^2*log(n))，朴素算法遍历搜索在判断累加，所以效率很低；

资料   NanoApe's Blog   ACdreamers

莫比乌斯反演:利用整与分之间的可逆来由整体利用容斥原理得到分量的值；这就是用容易求解的F[n]通过莫比乌斯公式(函数mu[])得到分量f[n]的值；

如本题中:F[n]表示公倍数是n的倍数的个数，f[n]表示公倍数为n的个数；即F[d] = Σf[n] ,(d|n);

同样每个F[d](d|n)中都含有f[n]的分量，所以可以使用容斥原理来求解；对应下面两个公式；

公式1:

公式2：

细节:这道题b,d并不相等，由于只是组合不是排列，这需要两次求之后去重；开始直接线性筛法把mu[]预处理出来即可；

坑:以前线性筛素数，一直用的p[j] < MAXN/i没事，今天各种WA...打击真大。。之后改成p[j]*i < MAXN也没做1LL*处理都没事。。

时间复杂度为:O(n)

#include<iostream>
#include<cstdio>
#include<cstring>
#include<string.h>
#include<algorithm>
#include<vector>
#include<cmath>
#include<stdlib.h>
#include<time.h>
#include<stack>
#include<set>
#include<map>
#include<queue>
using namespace std;
#define rep0(i,l,r) for(int i = (l);i < (r);i++)
#define rep1(i,l,r) for(int i = (l);i <= (r);i++)
#define rep_0(i,r,l) for(int i = (r);i > (l);i--)
#define rep_1(i,r,l) for(int i = (r);i >= (l);i--)
#define MS0(a) memset(a,0,sizeof(a))
#define MS1(a) memset(a,-1,sizeof(a))
#define MSi(a) memset(a,0x3f,sizeof(a))
#define inf 0x3f3f3f3f
#define lson l, m, rt << 1
#define rson m+1, r, rt << 1|1
typedef pair<int,int> PII;
#define A first
#define B second
#define MK make_pair
typedef __int64 ll;
template<typename T>
void read1(T &m)
{
    T x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
    m = x*f;
}
template<typename T>
void read2(T &a,T &b){read1(a);read1(b);}
template<typename T>
void read3(T &a,T &b,T &c){read1(a);read1(b);read1(c);}
const int MAXN = 100007;
int mu[MAXN],vis[MAXN],p[MAXN];
void mobius()
{
    mu[1] = 1;//定义
    rep0(i,2,MAXN){
        if(vis[i] == 0){
            mu[i] = -1;
            p[++p[0]] = i;
        }
        for(int j = 1;j <= p[0] && 1LL*p[j]*i < MAXN;j++){ // p[j] < MAXN/i WA了
            vis[i*p[j]] = 1;
            if(i%p[j]) mu[i*p[j]] = -mu[i];
            else{
                mu[i*p[j]] = 0;
                break;
            }
        }
    }
}
int main()
{
    int T,kase = 1,a,b,c,d,k;
    mobius();
    read1(T);
    while(T--){
        ll ans = 0,tmp = 0;
        read2(a,b);read3(c,d,k);
        if(k == 0){
            printf("Case %d: %I64d\n",kase++,ans);
            continue;
        }
        b /= k,d /= k;//**这样更简便
        if(b > d) swap(b,d);
        rep1(i,1,b) ans += 1LL*mu[i]*(b/i)*(d/i);
        rep1(i,1,b) tmp += 1LL*mu[i]*(b/i)*(b/i);
        printf("Case %d: %I64d\n",kase++,ans - tmp/2);
    }
    return 0;
}