HDU 4675-GCD of Sequence（莫比乌斯反演+组合数+逆元）

题意：给你n,m,k,然后给你一个长为n的序列，让你只能改变k个数的情况下，使得gcd(a1.....an)=d，d从1到m，问你相应的方案数。

题解：首先我们要知道若干个数的gcd=d的前提是这几个数都至少是d的倍数，不然无法满足情况，齐次因为要保证只能改变k个数，我们可以开一个num数组记录每个数的倍数出现的个数，假如当前枚举到数字i，如果num[i]<n-k,呢答案一定是0，想都不用想，因为改变了k个数后还有些数不是i的倍数。剩下的就是符合条件的情况了，我们可以用组合数求出所有情况，并用莫比乌斯反演算出贡献，然后求和既是答案。

我们设当前i的倍数为cnt[i],有n个数，需要改变k个数，公式为：

                                                    


因为n很大，需要用到逆元，可以参考我的模板：http://blog.csdn.net/haut_ykc/article/details/52136319

一个很坑的点是做题是忘了mu函数可能是-1，然后取模时忘了保证非负，然后wa了一上午。。。。。

#include<map>
#include<stack>
#include<queue>
#include<vector>
#include<string>
#include<math.h>
#include<stdio.h>
#include<iostream>
#include<string.h>
#include<stdlib.h>
#include<algorithm>
#include<functional>
using namespace std;
typedef long long ll;
#define inf 1000000000
#define mod 1000000007
#define maxn 560050
#define lowbit(x) (x&-x)
#define eps 1e-9
ll a[maxn] = { 1,1 }, b[maxn], mu[maxn];
ll f[maxn], cnt, ans[maxn], num[maxn], fac[maxn], inv[maxn];
ll q(ll x, ll y)
{
ll res = 1;
while (y)
{
if (y % 2)
res = res*x%mod;
x = x*x%mod;
y /= 2;
}
return res;
}
void init()
{
ll i, j;mu[1] = 1;
for (i = 2;i<maxn;i++)
{
if (a[i] == 0)
b[++cnt] = i, mu[i] = -1;
for (j = 1;j <= cnt && i*b[j] <= maxn;j++)
{
a[b[j] * i] = 1;
if (i%b[j] == 0)
{
mu[b[j] * i] = 0;
break;
}
else
mu[b[j] * i] = -mu[i];
}
}
fac[0] = inv[0] = 1;
for (i = 1;i <= 400005;i++)
{
fac[i] = fac[i - 1] * i%mod;
inv[i] = q(fac[i], mod - 2)%mod;
}
}
int main(void)
{
init();
ll i, x, k, p, n, m, j;
while (scanf("%lld%lld%lld",&n,&m,&k)!=EOF)
{
memset(num, 0, sizeof(num));
for (i = 1;i <= n;i++)
{
scanf("%lld", &x);
num[x]++;
}
for (i = 1;i <= m;i++)
for (j = i + i;j <= m;j += i)
num[i] += num[j];
for (i = 1;i <= m;i++)
{
if (num[i] < n - k)
{
f[i] = 0;
continue;
}
ll tmp = fac[num[i]] % mod * inv[n - k] % mod*inv[k + num[i] - n] % mod;
tmp = tmp*q(m / i - 1, num[i] + k - n) % mod*q(m / i, n - num[i]) % mod;
f[i] = tmp%mod;
}
for (i = 1;i <= m;i++)
{
ll sum = 0;
for (j = i;j <= m;j += i)
sum = ((sum + mu[j / i] * f[j] + mod) % mod) % mod;
ans[i] = sum;
}
printf("%lld", ans[1]);
for (i = 2;i <= m;i++)
printf(" %lld", ans[i]);
printf("\n");
}
return 0;
}