hdu4675 GCD of Sequence 莫比乌斯+组合数学

/**
题目：hdu4675 GCD of Sequence
链接：http://acm.hdu.edu.cn/showproblem.php?pid=4675
题意：给定n个数的a数组，以及m，k；
构造满足1<=bi<=m，和a数组恰好k个位置ai!=bi的b数组。
输出b数组所有数的gcd分别为1~m的数组个数。

思路：

f(n)表示gcd==n的数组个数。
g(n)表示gcd是n的倍数的数组个数。

f(n) = sigma[n|d]mu[d/n]*g(d);

如何求g(d)呢？

如果没有k，显然是g(d)=(M/d)^n;

可问题是存在k..... 必须满足所有的数都是d的倍数。 且有k个bi与ai不相同。

有M/d个数是d的倍数。 如果a数组有cnt个d的倍数。

那么剩下的n-cnt(如果n-cnt>k那么无解)个数必须变成d的倍数，有(M/d)^(n-cnt)种情况;

还剩下k-(n-cnt)个数需要从a数组cnt个是d的倍数中改变。有C(cnt,k-(n-cnt))*(M/d-1)^(k-(n-cnt));

所以g(d) = (M/d)^(n-cnt)*C(cnt,k-(n-cnt))*(M/d-1)^(k-(n-cnt)); (n-cnt<=k)
    g(d) = 0; (n-cnt>k)

    C(n,m) = n!/((n-m)!*m!)
*/
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <set>
#include <iostream>
#include <vector>
#include <map>
using namespace std;
typedef long long LL;
#define ms(x,y) memset(x,y,sizeof x)
typedef pair<int, int> P;
const LL INF = 1e10;
const int mod = 1e9 + 7;
const int maxn = 3e5 + 100;
int prime[maxn], tot, not_prime[maxn];
int mu[maxn], cnt[maxn];
LL fac[maxn], inv[maxn], f[maxn];
int n, m, k;
void init()
{
    inv[0] = inv[1] = 1;
    for(int i = 2; i < maxn; i++){
        inv[i] = (mod-mod/i)*inv[mod%i]%mod;
    }
    fac[0] = fac[1] = 1;
    for(int i = 2; i < maxn; i++){
        fac[i] = fac[i-1]*i%mod;
        inv[i] = inv[i]*inv[i-1]%mod;
    }
}
LL Pow(LL x,int y)
{
    LL p = 1;
    while(y){
        if(y&1) p = p*x%mod;
        x = x*x%mod;
        y>>=1;
    }
    return p;
}
void mobius()
{
    mu[1] = 1;
    tot = 0;
    for(int i = 2; i < maxn; i++){
        if(!not_prime[i]){
            mu[i] = -1;
            prime[++tot] = i;
        }
        for(int j = 1; prime[j]*i<maxn; j++){
            not_prime[prime[j]*i] = 1;
            if(i%prime[j]==0){
                mu[prime[j]*i] = 0;
                break;
            }
            mu[prime[j]*i] = -mu[i];
        }
    }
}

LL get(int d)
{
    if(n-cnt[d]>k) return 0;
    return Pow((LL)m/d,n-cnt[d])*fac[cnt[d]]%mod*inv[n-k]%mod*inv[k-n+cnt[d]]%mod*Pow((LL)m/d-1,k-n+cnt[d])%mod;
}
int main()
{
    //freopen("YYnoGCD.in","r",stdin);
    //freopen("YYnoGCD.out","w",stdout);
    //freopen("in.txt","r",stdin);
    int T;
    mobius();
    init();
    while(scanf("%d%d%d",&n,&m,&k)!=EOF)
    {
        ms(cnt,0);
        int x;
        for(int i = 1; i <= n; i++){
            scanf("%d",&x);
            cnt[x]++;
        }
        for(int i = 1; i <= m; i++){
            for(int j = 2*i; j <= m; j+=i){
                cnt[i] += cnt[j];
            }
        }
        ms(f,0);

        for(int i = 1; i <= m; i++){
            for(int j = i; j <= m; j+=i){
                f[i] = (f[i]+mu[j/i]*get(j)+mod)%mod;
            }
        }
        for(int i = 1; i < m; i++) printf("%lld ",f[i]);
        printf("%lld\n",f[m]);
    }
    return 0;
}