bzoj4305 数列的GCD

传送门
由于要求所有的$d$的答案，我们考虑容斥。
用$f(d)$表示$gcd=d$时的答案。用$F(d)$表示$d|gcd$的答案之和。那么我们有：
$$F(n)=\sum _{n|d}f(d)\iff f(n)=\sum _{n|d}\mu (\frac{d}{n})F(d)$$
上式为莫比乌斯反演的一种形式。
那么现在考虑求$F(d)。$令$cn{t}_d$表示$a$数组中有多少个$a_i$为$d$的倍数。
现在每个数都必须是$d$的倍数。也就是说对于那些与$a$不同的$b$，有$\frac{m}{d}$种选择。
现在我们要求有$k$个不相同，也就是说有$N-k$个是相同的。
那么这$N-k$个相同的一定从$cn{t}_d$中选（$d$的倍数）。选不够贡献就是$0$。
相同的选完之后，选不相同的。那么有$cn{t}_d-(N-k)$个只能选$(\frac{m}{d}-1)$个。因为不能和$a$相同。剩下的就可以随便选了。于是我们得到$F$的表达式：
$$F(d)=\left(\genfrac{}{}{0ex}{}{cn{t}_d}{N-k}\right)\ast (\frac{m}{d}-1{)}^{cn{t}_d-(N-k)}\ast (\frac{m}{d}{)}^{N-cn{t}_d}$$
那么我们可以先$O(mlog(m))$处理出$cnt$数组，然后在$O(log)$的时间内求出$F$（快速幂）。
然后对于每个$f(i)$，枚举$i$的倍数。于是$\frac{m}{1}+\frac{m}{2}+\frac{m}{3}+...$，约为$O(mlog(m))$。

写代码的时候又犯了一个sb错误。$++vis[x]$写成了$vis[x]=1$，无语。
看来还是变量名字取得不好啊，以后要取恰当一点。。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod=1e9+7,maxn=3e5+10;
int fac[maxn],ifac[maxn];
int n,m,k,cnt[maxn],vis[maxn],x,mx;
int P[maxn],mu[maxn],mark[maxn],tot=0;
inline int add(int x,int y){return x+y>=mod?x+y-mod:x+y;}
inline int dec(int x,int y){return x-y<0?x-y+mod:x-y;}
inline int mul(int x,int y){return (ll)x*y%mod;}
inline int quickpow(int a,int b,int ret=1){for(;b;b>>=1,a=mul(a,a))if(b&1)ret=mul(ret,a);return ret;}
inline int C(int n,int m){
	if(m>n) return 0;
	return mul(fac[n],mul(ifac[m],ifac[n-m]));
}
inline int F(int d){
	if(n-k>cnt[d]) return 0;
	return mul(C(cnt[d],n-k),mul(quickpow((m/d)-1,cnt[d]-n+k),quickpow(m/d,n-cnt[d])));
}
inline void linear_sieves(int N=maxn){
	mark[0]=mark[1]=mu[1]=1;
	for(int i=2;i<N;++i){
		if(!mark[i]) P[++tot]=i,mu[i]=-1;
		for(int j=1;j<=tot&&i*P[j]<N;++j){
			mark[i*P[j]]=1;
			if(i%P[j]) mu[i*P[j]]=-mu[i];
			else break;
		}
	}fac[0]=fac[1]=1;
	for(int i=2;i<N;++i) fac[i]=mul(fac[i-1],i);
	ifac[N-1]=quickpow(fac[N-1],mod-2);
	for(int i=N-2;i>=0;--i) ifac[i]=mul(ifac[i+1],i+1);
}
inline int read(){
	int x=0;char ch=getchar();
	while(!isdigit(ch)) ch=getchar();
	while(isdigit(ch)) x=(x<<3)+(x<<1)+ch-'0',ch=getchar();
	return x;
}
int main(){
	linear_sieves(),n=read(),m=read(),k=read();
	for(int i=1;i<=n;++i) ++vis[x=read()],mx=max(mx,x);
	for(int i=1;i<=mx;++i)
		for(int j=i;j<=mx;j+=i)
			cnt[i]+=vis[j];
	for(int i=1;i<=m;++i){
		int ans=0;
		for(int j=i,t=1;j<=m;j+=i,++t){
			if(mu[t]==0) continue;
			if(mu[t]==1) ans=add(ans,F(j));
			if(mu[t]==-1)ans=dec(ans,F(j));
		}printf("%d ",ans);
	}
}