hdu4675 GCD of Sequence

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=4675

题意：

　　给定一个长度为n的序列a，且 1<=a[i]<=m,求分别有多少个序列b,使得GCD(b[1],b[2],...b[n])=x (1<=x<=m),且正好有k个b[i]！=a[i]。

分析：

　　莫比乌斯反演，主要是确定F(x)。

　　用F(x)表示gcd为x的倍数的方案数，f(x)表示gcd为x的方案数。

　　先考虑F(d)怎么计算。可以把a数组中的数分成两类，第一类是必须对应下标不等的，即a[i]不是d的倍数)，其他的就是第二类。

　　假设第二类的数量是p，第一类的数量就是n−p，因为要选择k个不同的，第一类必须不同，所以需要在第二类中选择k−n+p个，而b数组中每一个数都有[m/p]种选择。

　　所以最终的结果就是F(d)=C(p，k-n+p) ∗ ( ([m/d]−1)^(k−n+p) )∗( [m/d]^(n−p) )，然后暴力计算每一个f(d)就好了。总的复杂度是n(logn)。

代码：

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>

using namespace std;
const int maxn=300010;
const long long mod=1e9+7;

int n,m,k;
int a[maxn];
int mu[maxn];
int vis[maxn];
int prime[maxn];
int cnt;
int num[maxn];
long long jie[maxn],ni[maxn];
long long F[maxn];
long long res[maxn];

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

long long power(long long a,long long n,long long m)
{
    long long ans=1,tmp=a%m;
    while(n)
    {
        if(n&1)
            ans=ans*tmp%m;
        tmp=tmp*tmp%m;
        n=n/2;
    }
    return ans;
}

long long C(long long n,long long m)
{
    if(n==0)
        return 1;
    return jie[n]*ni[m]%mod*ni[n-m]%mod;
}

int main()
{
    init();
    ni[0]=1;
    jie[0]=1;
    for(int i=1;i<maxn;i++)
    {
        jie[i]=jie[i-1]*i%mod;
        ni[i]=power(jie[i],mod-2,mod);
    }
    while(~scanf("%d%d%d",&n,&m,&k))
    {
        for(int i=1;i<=n;i++)
            scanf("%d",&a[i]);
        memset(num,0,sizeof(num));
        for(int i=1;i<=n;i++)
            num[a[i]]++;
        for(int i=1;i<=m;i++)
        {
            long long p=0;
            for(int j=1;j*i<=m;j++)
                p+=num[i*j];
            if(k-n+p<0)
                F[i]=0;
            else
                F[i]=C(p,k-n+p)*power(m/i-1,k-n+p,mod)%mod*power(m/i,n-p,mod)%mod;
        }
        for(int i=1;i<=m;i++)
        {
            long long ans=0;
            for(int j=1;i*j<=m;j++)
            {
                ans+=mu[j]*F[i*j];
                ans=(ans%mod+mod)%mod;
            }
            if(i==1)
                printf("%lld",ans);
            else
                printf(" %lld",ans);
        }
        printf("
");
    }
    return 0;
}