HDU - 4675 GCD of Sequence (容斥组合推公式)

题目描述:

点击打开链接

给一个序列 {an}（ai <= M），问有多少种 {bn} ，满足恰好 k 个位置与 {an} 不同，且 gcd of {bn} = d（1<= d <=M）。

很灵巧的题，我们首先把a中d的倍数的数目统计出来，假设为x个，如果x<=n-k的话答案肯定是不存在的，问题是答案存在的情况怎么求，我们可以确定的是b中有n-x个位置肯定无法从a中直接得到d的倍数，那么对于这n-x个位置，每个位置能选择的数就有(M/d)种，那么这部分的方案就是(M/d)^(n-x)，但是由于有k的限制，a中的x个数肯定是无法都直接搬到b中的，所以要从中选n-k个数搬到b中，这个地方我们可以反过来想，我们要选出k-(n-x)个位置，说明这些位置不能直接用a中的数，那么我们能用的数就只有(M/d-1)种，所以这部分的方案数就是C(x,k-(n-x))*(M/d-1)^(k-(n-x))，最后相乘就是此时的答案，但是你在算1-M的答案肯定是算重了的，比如算2的时候肯定把2的倍数的结果都算了一遍，所以我们可以倒过来处理，从M算到1这样每次算一个数的时候就把已经算出来的这个数倍数的答案减去就行了。

AC代码：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<cmath>
#include<stack>
#include<vector>
#include<set>
#include<queue>
#include<algorithm>
using namespace std;
const int MAXM=1e6;
const long long MOD=1000000007;

long long ans[MAXM+10];
long long cnt[MAXM+10];
long long inv[MAXM+10];
long long fact[MAXM+10];

long long qpow(long long x,long long y)
{
    long long res=1;
    while(y)
    {
        if (y&1) res=(res*x)%MOD;
        x=(x*x)%MOD;
        y=y>>1;
    }
    return res;
}
void init()
{
    fact[0]=1;
    inv[0]=1;
    for (long long i=1;i<=MAXM;i++)
    {
        fact[i]=(i*fact[i-1])%MOD;
        inv[i]=qpow(fact[i],MOD-2);
    }
}
long long comb(long long x,long long y)
{
    return ((fact[x]*inv[y])%MOD)*inv[x-y]%MOD;
}

int n,m,k;

int main()
{
    init();
    while(scanf("%d%d%d",&n,&m,&k)!=EOF)
    {
        memset(cnt,0,sizeof(cnt));
        int x;
        for (int i=0;i<n;i++)
        {
            scanf("%d",&x);
            cnt[x]++;
        }
        for (int i=m;i>=1;i--)
        {
            long long tot=0;
            for (int j=i;j<=m;j=j+i) tot=tot+cnt[j];
            if (n-tot>k)
            {
                ans[i]=0;
                continue;
            }
            ans[i]=qpow(m/i,n-tot)*qpow(m/i-1,k-(n-tot))%MOD;
            ans[i]=ans[i]*comb(tot,k-(n-tot))%MOD;
            for (int j=i+i;j<=m;j=j+i)
                ans[i]=(ans[i]-ans[j]+MOD)%MOD;
        }
        for (int i=1;i<=m;i++)
        {
            if (i==m) printf("%I64d\n",ans[i]);
            else printf("%I64d ",ans[i]);
        }
    }
    return 0;
}