[HAOI2018]染色

[HAOI2018]染色

标签： 容斥原理 NTT 多项式求逆

---

题目链接

提供一种与网上大多数博客中不太一样的做法。

Solution

这题很明显是容斥原理。
如果我们直接计算不考虑重复，那么可以写出式子

[sum_{i=0}^{min(n/s,m) } C_m^i C_n^{is} {(is)! over (s!)^i} w_i (m-i)^{n-is} ]

想法是枚举恰好出现了i种颜色出现s次，然后其他的随便取，但是这样很明显是会重复的。
因为在剩下的随便取的过程中，可能会又出现一些颜色恰好出现了s次。

我们考虑给出现种数i加一个容斥系数(f_i)，显然，假如有j种颜色出现了s次，那么按照上面的式子，j的贡献会被i计算(C_j^i)次。
我们只需要对于每一个(w_i),都有$$w_i=sum_{j=0}^i C_i^j f_j$$
那么这样算就不会重了。
最终的式子是$$sum_{i=0}^{min(n/s,m) } C_m^i C_n^{is} {(is)! over (s!)^i} f_i (m-i)^{n-is} $$

现在唯一的难点就是如何算(f)，直接递推是(O(n^2))的。

注意到是一个卷积的形式，

[sum_{j=0} {f_j over j!} {1 over (i-j)!} = { w_iover i!} ]

令(F(x)=sum_{i=0}^{infty} {f_iover i!}x^i,G(x)=sum_{i=0}^{infty} {1over i!}x^i,W(x)=sum_{i=0}^{infty} {w_i over i!}x^i)

显然有(F(x)G(x)=W(x)),即(F(x)={W(x) over G(x)})

所以只需要分治FFT或者多项式求逆即可。

Code

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cmath>
#include<queue>
#include<stack>
#include<set>
#include<map>
using namespace std;
#define ll long long
#define RG register
#define REP(i,a,b) for(RG int i=(a),_end_=(b);i<=_end_;i++)
#define DREP(i,a,b) for(RG int i=(a),_end_=(b);i>=_end_;i--)
#define EREP(i,a) for(int i=start[(a)];i;i=e[i].next)
inline int read()
{
	int sum=0,p=1;char ch=getchar();
	while(!(('0'<=ch && ch<='9') || ch=='-'))ch=getchar();
	if(ch=='-')p=-1,ch=getchar();
	while('0'<=ch && ch<='9')sum=sum*10+ch-48,ch=getchar();
	return sum*p;
}

const int maxn=3e5+20;
const int mod=1004535809;
const int maxp=1e7;

int n,m,s,g[maxn];
int jc[maxp+20],inv[maxp+20],jcn[maxp+20];

int rev[maxn];
int st[maxn];

inline int power(int a,int b)
{
	int ans=1;
	while(b)
	{
		if(b&1)ans=(ll)ans*a%mod;
		b>>=1;
		a=(ll)a*a%mod;
	}
	return ans;
}

inline void NTT(int *p,int n,int op)
{
	int l=0;
	while((1<<l)<n)l++;
	REP(i,0,n-1)rev[i]=(rev[i>>1]>>1)|((i&1)<<l-1);
	REP(i,1,n-1)if(i<rev[i])swap(p[i],p[rev[i]]);
	for(int i=1;i<n;i<<=1)
	{
		int w=power(3,(mod-1)/(i<<1));
		st[0]=1;REP(j,1,i-1)st[j]=(ll)st[j-1]*w%mod;
		for(int j=0;j<n;j+=i<<1)
		{
			for(int k=0;k<i;k++)
			{
				int x=p[j+k],y=(ll)p[i+j+k]*st[k]%mod;
				p[j+k]=x+y;if(p[j+k]>=mod)p[j+k]-=mod;
				p[i+j+k]=x-y;if(p[i+j+k]<0)p[i+j+k]+=mod;
			}
		}
	}
	if(op==-1)
	{
		int inv=power(n,mod-2)%mod;
		REP(i,0,n-1)p[i]=(ll)p[i]*inv%mod;
		reverse(p+1,p+n);
	}
}

inline void init()
{
	n=read();m=read();s=read();
	REP(i,0,min(m,n/s))g[i]=read();
	jc[0]=jc[1]=jcn[0]=jcn[1]=inv[1]=1;
	REP(i,2,max(n,m))jc[i]=(ll)i*jc[i-1]%mod,inv[i]=(ll)(mod-mod/i)*inv[mod%i]%mod,jcn[i]=(ll)inv[i]*jcn[i-1]%mod;
	
}

int f[maxn],h[maxn],H[maxn];
int A[maxn],B[maxn];

void Inv(int *p,int *q,int n)
{
	if(n==1)return q[0]=p[0],void();
	Inv(p,q,n>>1);
	REP(i,0,(n<<1)-1)A[i]=B[i]=0;
	REP(i,0,n-1)A[i]=q[i],B[i]=p[i];
	NTT(A,n<<1,1);NTT(B,n<<1,1);
	REP(i,0,(n<<1)-1)A[i]=(ll)A[i]*A[i]%mod*B[i]%mod;
	NTT(A,n<<1,-1);
	REP(i,0,n-1)q[i]=((ll)2*q[i]-A[i]+mod)%mod;
}

inline void get_f()
{
	int N=1,lim=min(m,n/s);
	while(N<=lim)N<<=1;
	REP(i,0,lim)g[i]=(ll)g[i]*jcn[i]%mod;
	REP(i,0,lim)h[i]=jcn[i];
	Inv(h,H,N);
	N<<=1;
	NTT(H,N,1);NTT(g,N,1);
	REP(i,0,N-1)f[i]=(ll)H[i]*g[i]%mod;
	NTT(f,N,-1);
	REP(i,0,lim)f[i]=(ll)f[i]*jc[i]%mod;
}

inline int C(int n,int m){ if(n<m)return 0;return (ll)jc[n]*jcn[m]%mod*jcn[n-m]%mod;}

inline void doing()
{
	int ans=0,S=1;
	REP(i,0,min(m,n/s))
	{
		ans=(ans+(ll)C(m,i)*C(n,i*s)%mod*jc[i*s]%mod*S%mod*f[i]%mod*power(m-i,n-i*s))%mod;
		S=(ll)S*jcn[s]%mod;
	}
	printf("%d
",ans);
}

int main()
{
#ifndef ONLINE_JUDGE
	freopen("paint.in","r",stdin);
	freopen("paint.out","w",stdout);
#endif
	init();
	get_f();
	doing();
	return 0;
}