HDU 4675 GCD of Sequence （莫比乌斯反演）

题意：给$n$个数，$a_1,a_2,...,a_n$，其中$1<=a_i<=M$，找出所有的满足条件的序列$b_1,b_2,...,b_n$使得$gcd(b_1,b_2,...,b_n)=d$，其中$1<=b_i<=M$，且序列{$b_i$}中有且仅有$k$个$b_i!=a_i$。

$f(d)$表示$d=gcd(b_1,b_2,...,b_n)$的所有满足的组合的总数。
$F(d)$表示$d|gcd(b_1,b_2,...,b_n)$的所有满足的组合的总数。

题目要求，求出每一个$f(x)$值。

$f(x)=\sum_{x|d}\mu({d\over{x}})*F(d)$。

只要求出所有的$F(d)$即可。

用$cnt(d)$表示原序列中$d$的倍数的个数。

那么$F(d)={cnt(d)\choose{n-k}}*({M\over{d}}-1)^{cnt(d)+k-n}*({M\over{d}})^{n-cnt(d)}$。

算一下复杂度，预处理出所有$x!$即，$x!$的逆元，那么求$n\choose{x}$的复杂度为$o(1)$，快速幂的复杂度为$o(logn)$，那么求$F(d)$的复杂度为$o(nlogn)$。

求一次$f(x)$的复杂度为$o({n\over{x}})$求所有$f(x)$的复杂度为$o(\sum_{x=1}^n{n\over{x}})=o(nlogn)$。

所以总的复杂度为$o(nlogn)$。

下面是代码：

#include <bits/stdc++.h>
#define LL long long
#define FOR(i,x,y)  for(int i = x;i < y;++ i)
#define IFOR(i,x,y) for(int i = x;i > y;-- i)

using namespace std;

const int maxn = 300010;
const int Mod = 1000000007;

LL fac[maxn],inv[maxn];

int prime[maxn],mu[maxn];
bool check[maxn];

void Mobius(){
    memset(check,false,sizeof(check));
    prime[0] = 0;
    mu[1] = 1;
    FOR(i,2,maxn){
        if(!check[i]){
            prime[++prime[0]] = i;
            mu[i] = -1;
        }
        FOR(j,1,prime[0]+1){
            if(i*prime[j] >= maxn)  break;
            check[i*prime[j]] = true;
            if(i%prime[j]){
                mu[i*prime[j]] = -mu[i];
            }
            else{
                mu[i*prime[j]] = 0;
                break;
            }
        }
    }
}

void Gcd(LL a,LL b,LL& d,LL& x,LL& y){
    if(!b)  {d = a;x = 1;y = 0;}
    else{Gcd(b,a%b,d,y,x);y -= x*(a/b);}
}

LL Inv(LL a,LL n){
    LL d,x,y;
    Gcd(a,n,d,x,y);
    return d == 1 ? (x+n)%n : -1;
}

void init(){
    Mobius();
    fac[0] = 1;
    inv[0] = 1;
    FOR(i,1,maxn){
        fac[i] = fac[i-1]*i;
        fac[i] %= Mod;
        inv[i] = Inv(fac[i],Mod);
    }
}

LL quickpow(LL a,LL n,LL m){
    LL ans=1;
    while(n){
        if(n&1) ans = (ans*a)%m;
        a = (a*a)%m;
        n>>=1;
    }
    return ans;
}

LL F[maxn];
int cnt[maxn];

int n,m,k;

void work(){
    memset(cnt,0,sizeof(cnt));
    int num;
    FOR(i,0,n){
        scanf("%d",&num);
        cnt[num] ++;
    }
    FOR(i,1,m+1){
        for(int j = i+i;j <= m;j += i){
            cnt[i] += cnt[j];
        }
    }
    FOR(d,1,m+1){
        if(cnt[d] < n-k)    {F[d] = 0;continue;}
        LL t = (inv[n-k]*inv[cnt[d]+k-n])%Mod;
        t = (fac[cnt[d]]*t)%Mod;
        F[d] = t*((quickpow((m/d)-1,cnt[d]+k-n,Mod)*quickpow(m/d,n-cnt[d],Mod))%Mod)%Mod;
    }
    LL res;
    FOR(x,1,m+1){
        res = 0;
        for(int d = x;d <= m;d += x){
            res += mu[d/x]*F[d];
            res = (res + Mod)%Mod;
        }
        if(x == m){
            printf("%I64d\n",res);
        }
        else
            printf("%I64d ",res);
    }
}

int main()
{
    //freopen("test.in","r",stdin);
    init();
    while(~scanf("%d%d%d",&n,&m,&k)){
        work();
    }
    return 0;
}