HDU - 4675 GCD of Sequence (莫比乌斯反演+组合数学)

题意：给出序列[a1..aN],整数M和k，求对1-M中的每个整数d，构建新的序列[b1...bN],使其满足：
1. \(1 \le bi \le M\)
2. \(gcd(b 1, b 2, …, b N) = d\)
3. 恰好有k个位置 \(bi!=ai\)
求对每个d，有多少种满足条件的序列
分析：对于前两个条件，就是单纯的莫比乌斯反演。
令\(F(d) = [d|gcd(b1...bN)]\)
\(f(d) = [gcd(b1...bN)]=d]\)
则$f(n) = \sum_{x|d}u(\frac{d}{x})F(d) \( 而该处又有限制恰好k个数与原序列不同，则考虑先从原序列是d的倍数的数中选出N-k个数,再在原序列不是d的倍数的位置上随意放置d的倍数，再将原序列中是d的倍数但在第一步中没有被选择的数放置与其不同的d的倍数。 所以\)F(d) =\dbinom{cnt(d)}{n-k} * (\frac{M}{d})^{n-cnt(d)} * (\frac{M}{d}-1)^{k-n+cnt(d)}$
需要预处理出莫比乌斯函数和阶乘逆

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int mod =1e9+7;
const int maxn=300000+5;
bool vis[maxn];
int prime[maxn],mu[maxn];
LL fac[maxn],inv[maxn];

LL qpow(LL a,LL n)
{
    LL res=1;
    while(n){
        if(n&1) res = res*a%mod;
        a = a*a%mod;
        n>>=1;
    }
    return res;
}

void pre()
{
    fac[0] = fac[1] = 1;
    for(int i=2;i<maxn;++i) fac[i] = (fac[i-1]*i)%mod;
    inv[maxn-1] = qpow(fac[maxn-1],mod-2);
    for(int i=maxn-2;i>=0;--i) inv[i] = ( inv[i+1]*(i+1) )%mod; 
}

LL Comb(LL n,LL k)
{
    if(n<k) return 0;
    else if(n==k) return 1;
    else return ((fac[n]*inv[k]%mod)*inv[n-k])%mod;
}

void init_mu(int n){
    int cnt=0;
    mu[1]=1;
    for(int i=2;i<n;i++){
        if(!vis[i]){
            prime[cnt++]=i;
            mu[i]=-1;
        }
        for(int j=0;j<cnt&&i*prime[j]<n;j++){
            vis[i*prime[j]]=1;
            if(i%prime[j]==0)   {mu[i*prime[j]]=0;break;}
            else { mu[i*prime[j]]=-mu[i];}
        }
    }
}

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

int main()
{
    #ifndef ONLINE_JUDGE
        freopen("in.txt","r",stdin);
        freopen("out.txt","w",stdout);
    #endif
    init_mu(maxn);
    pre();
    int N,M,k,tmp;
    while(scanf("%d %d %d",&N,&M,&k)==3){
        memset(cnt,0,sizeof(cnt));
        for(int i=1;i<=N;++i){
            scanf("%d",&tmp);
            cnt[tmp]++;
        }
        for(int i=1;i<=M;++i){              //打表处理d的倍数的个数
            for(int j=2*i;j<=M;j+=i){
                cnt[i]+=cnt[j];
            }
        }
        for(int d=1;d<=M;++d){
            if(cnt[d]<N-k) {
                F[d] = 0;
                continue;
            }
            LL tmp = Comb(cnt[d],N-k) * qpow(M/d,N-cnt[d]) %mod;
            F[d] = tmp * qpow(M/d-1,k-N+cnt[d]) %mod;
        }
        for(int t=1;t<=M;++t){
            LL res=0;
            for(int d = t;d<=M;d+=t){
                res = (res+mu[d/t]*F[d]+mod)%mod;
            }
            if(t==M) printf("%lld\n",res);
            else printf("%lld ",res);
        }
    }
    return 0;
}