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;
}