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