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)。

代码：

  1 #include<iostream>
  2 #include<cstdio>
  3 #include<cmath>
  4 #include<cstring>
  5
  6 using namespace std;
  7 const int maxn=300010;
  8 const long long mod=1e9+7;
  9
 10 int n,m,k;
 11 int a[maxn];
 12 int mu[maxn];
 13 int vis[maxn];
 14 int prime[maxn];
 15 int cnt;
 16 int num[maxn];
 17 long long jie[maxn],ni[maxn];
 18 long long F[maxn];
 19 long long res[maxn];
 20
 21 void init()
 22 {
 23     memset(vis,0,sizeof(vis));
 24     mu[1]=1;
 25     cnt=0;
 26     for(int i=2;i<maxn;i++)
 27     {
 28         if(!vis[i])
 29         {
 30             prime[cnt++]=i;
 31             mu[i]=-1;
 32         }
 33         for(int j=0;j<cnt&&i*prime[j]<maxn;j++)
 34         {
 35             vis[i*prime[j]]=1;
 36             if(i%prime[j])
 37                 mu[i*prime[j]]=-mu[i];
 38             else
 39             {
 40                 mu[i*prime[j]]=0;
 41                 break;
 42             }
 43         }
 44     }
 45 }
 46
 47 long long power(long long a,long long n,long long m)
 48 {
 49     long long ans=1,tmp=a%m;
 50     while(n)
 51     {
 52         if(n&1)
 53             ans=ans*tmp%m;
 54         tmp=tmp*tmp%m;
 55         n=n/2;
 56     }
 57     return ans;
 58 }
 59
 60 long long C(long long n,long long m)
 61 {
 62     if(n==0)
 63         return 1;
 64     return jie[n]*ni[m]%mod*ni[n-m]%mod;
 65 }
 66
 67 int main()
 68 {
 69     init();
 70     ni[0]=1;
 71     jie[0]=1;
 72     for(int i=1;i<maxn;i++)
 73     {
 74         jie[i]=jie[i-1]*i%mod;
 75         ni[i]=power(jie[i],mod-2,mod);
 76     }
 77     while(~scanf("%d%d%d",&n,&m,&k))
 78     {
 79         for(int i=1;i<=n;i++)
 80             scanf("%d",&a[i]);
 81         memset(num,0,sizeof(num));
 82         for(int i=1;i<=n;i++)
 83             num[a[i]]++;
 84         for(int i=1;i<=m;i++)
 85         {
 86             long long p=0;
 87             for(int j=1;j*i<=m;j++)
 88                 p+=num[i*j];
 89             if(k-n+p<0)
 90                 F[i]=0;
 91             else
 92                 F[i]=C(p,k-n+p)*power(m/i-1,k-n+p,mod)%mod*power(m/i,n-p,mod)%mod;
 93         }
 94         for(int i=1;i<=m;i++)
 95         {
 96             long long ans=0;
 97             for(int j=1;i*j<=m;j++)
 98             {
 99                 ans+=mu[j]*F[i*j];
100                 ans=(ans%mod+mod)%mod;
101             }
102             if(i==1)
103                 printf("%lld",ans);
104             else
105                 printf(" %lld",ans);
106         }
107         printf("\n");
108     }
109     return 0;
110 }