HAOI2018染色——容斥

最新推荐文章于 2025-01-31 15:42:28 发布

原创最新推荐文章于 2025-01-31 15:42:28 发布 · 226 阅读

1 ·

CC 4.0 BY-SA版权

容斥原理同时被 2 个专栏收录

9 篇文章

订阅专栏

分治fft

3 篇文章

订阅专栏

本文介绍了一种利用NTT（Number Theoretic Transform）优化组合计数问题的方法，特别是针对颜色方案计数的问题。通过分治FFT和二项式反演技术，实现了高效的计算方案。代码示例展示了具体的实现细节。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

题目大意

loj

思路

设 $f_i$ 表示至少出现了i种颜色的方案数

$\begin{aligned} f_i&={m \choose i}\times \frac{(s\times i)!}{(s!)^{i}}\times {n\choose s\times i}\times (m-i)^{n-s\times i}\\ f_i&={m \choose i}\times \frac{n!}{(s!)^{i}\times (n-s\times i)!}\times (m-i)^{n-s\times i}\\ \end{aligned}$

设 $g_i$ 表示恰好出现了i种颜色的方案数，不难得到 $g_i$ 的容斥式子。

$\begin{aligned} g_i=&f_i-\sum_{j=i+1}^{m}g_{j}\times {j\choose i}\\ g_i=&f_i-\frac{\sum_{j=i+1}^{m}\times \frac{g_j\times j!}{(j-i)!}}{i!} \end{aligned}$

上面是一个经典的分治FFT形式的式子，使用分治FFT可以做到 $n\log^2n$ 。

其实可以更加优化，上面的式子同样也是一个经典的二项式反演。

$\begin{aligned} f_i&=\sum_{j=i}^{m}{j\choose i}\times g_{j}\\ g_i&=\sum_{j=i}^{m}(-1)^{j-i}\times {j\choose i}\times f_{j} \end{aligned}$

直接用NTT一遍求出g的值即可。

/*=======================================
 * Author : ylsoi
 * Time : 2019.3.22
 * Problem : loj2527_ntt
 * E-mail : ylsoi@foxmail.com
 * ====================================*/
#include<bits/stdc++.h>

#define REP(i,a,b) for(int i=a,i##_end_=b;i<=i##_end_;++i)
#define DREP(i,a,b) for(int i=a,i##_end_=b;i>=i##_end_;--i)
#define debug(x) cout<<#x<<"="<<x<<" "
#define fi first
#define se second
#define mk make_pair
#define pb push_back
typedef long long ll;

using namespace std;

void File(){
	freopen("loj2527_ntt.in","r",stdin);
	freopen("loj2527_ntt.out","w",stdout);
}

template<typename T>void read(T &_){
	_=0; T f=1; char c=getchar();
	for(;!isdigit(c);c=getchar())if(c=='-')f=-1;
	for(;isdigit(c);c=getchar())_=(_<<1)+(_<<3)+(c^'0');
	_*=f;
}

string proc(){
	ifstream f("/proc/self/status");
	return string(istreambuf_iterator<char>(f),istreambuf_iterator<char>());
}

const int maxn=1e7+10;
const int maxm=1e5+10;
const int mod=1004535809;
int n,m,S,w[maxm],ans;
int fac[maxn],ifac[maxn];

int qpow(int x,int y=mod-2){
	int ret=1; x%=mod;
	while(y){
		if(y&1)ret=1ll*ret*x%mod;
		x=1ll*x*x%mod;
		y>>=1;
	}
	return ret;
}

void inc(int &x,int y){
	x+=y;
	if(x>=mod)x-=mod;
	else if(x<0)x+=mod;
}

int C(int x,int y){
	if(x<0 || y<0 || x<y)return 0;
	return 1ll*fac[x]*ifac[y]%mod*ifac[x-y]%mod;
}

namespace Poly{
	int om[maxm<<2],dn[maxm<<2],lim,cnt;
	void dft(int *A,int ty){
		if(ty==-1)reverse(A+1,A+lim);
		REP(i,0,lim-1)if(i<dn[i])swap(A[i],A[dn[i]]);
		for(int len=1;len<lim;len<<=1){
			int w=om[len<<1];
			for(int L=0;L<lim;L+=len<<1){
				int wk=1;
				REP(i,L,L+len-1){
					int u=A[i],v=1ll*A[i+len]*wk%mod;
					A[i]=(u+v)%mod;
					A[i+len]=(u-v)%mod;
					wk=1ll*wk*w%mod;
				}
			}
		}
		if(ty==-1){
			int inv=qpow(lim);
			REP(i,0,lim-1)A[i]=(1ll*A[i]*inv%mod+mod)%mod;
		}
	}
}

void init(){
	read(n),read(m),read(S);
	REP(i,0,m)read(w[i]);
	fac[0]=1;
	int lim=max(n,m);
	REP(i,1,lim)fac[i]=1ll*fac[i-1]*i%mod;
	ifac[lim]=qpow(fac[lim]);
	DREP(i,lim-1,0)ifac[i]=1ll*ifac[i+1]*(i+1)%mod;
}

int f[maxm<<2],g[maxm<<2];

using namespace Poly;

int main(){
//	File();
	init();
	REP(i,0,m)if(S*i<=n)f[i]=1ll*C(m,i)*fac[S*i]%mod*
		qpow(ifac[S],i)%mod*C(n,S*i)%mod*qpow(m-i,n-S*i)%mod*fac[i]%mod;
	REP(i,0,m)g[i]=1ll*(i&1 ? -1 : 1)*ifac[i]%mod;
	reverse(g,g+m+1);

	lim=1,cnt=0;
	while(lim<=m+m)lim<<=1,++cnt;
	if(!cnt)cnt=1;
	om[lim]=qpow(3,(mod-1)/lim);
	for(int i=lim>>1;i;i>>=1)
		om[i]=1ll*om[i<<1]*om[i<<1]%mod;
	REP(i,0,lim-1)dn[i]=dn[i>>1]>>1|((i&1)<<(cnt-1));

	dft(f,1),dft(g,1);
	REP(i,0,lim-1)g[i]=1ll*g[i]*f[i]%mod;
	dft(g,-1);

	REP(i,0,m)inc(ans,1ll*g[i+m]*w[i]%mod*ifac[i]%mod);
	printf("%d\n",ans);

	return 0;
}

/*=======================================
 * Author : ylsoi
 * Time : 2019.3.21
 * Problem : loj2527
 * E-mail : ylsoi@foxmail.com
 * ====================================*/
#include<bits/stdc++.h>

#define REP(i,a,b) for(int i=a,i##_end_=b;i<=i##_end_;++i)
#define DREP(i,a,b) for(int i=a,i##_end_=b;i>=i##_end_;--i)
#define debug(x) cout<<#x<<"="<<x<<" "
#define fi first
#define se second
#define mk make_pair
#define pb push_back
typedef long long ll;

using namespace std;

void File(){
	freopen("loj2527.in","r",stdin);
	freopen("loj2527.out","w",stdout);
}

template<typename T>void read(T &_){
	_=0; T f=1; char c=getchar();
	for(;!isdigit(c);c=getchar())if(c=='-')f=-1;
	for(;isdigit(c);c=getchar())_=(_<<1)+(_<<3)+(c^'0');
	_*=f;
}

string proc(){
	ifstream f("/proc/self/status");
	return string(istreambuf_iterator<char>(f),istreambuf_iterator<char>());
}

const int maxn=1e7+10;
const int maxm=1e5+10;
const int mod=1004535809;
int n,m,S,w[maxm],ans;
int fac[maxn],ifac[maxn];

int qpow(int x,int y=mod-2){
	int ret=1; x%=mod;
	while(y){
		if(y&1)ret=1ll*ret*x%mod;
		x=1ll*x*x%mod;
		y>>=1;
	}
	return ret;
}

int C(int x,int y){
	if(x<0 || y<0 || x<y)return 0;
	return 1ll*fac[x]*ifac[y]%mod*ifac[x-y]%mod;
}

void inc(int &x,int y){
	x+=y;
	if(x>=mod)x-=mod;
	else if(x<0)x+=mod;
}

namespace Poly{
	int om[maxm<<2],dn[maxm<<2],lim,cnt;
	void dft(int *A,int ty){
		if(ty==-1)reverse(A+1,A+lim);
		REP(i,0,lim-1)if(i<dn[i])swap(A[i],A[dn[i]]);
		for(int len=1;len<lim;len<<=1){
			int w=om[len<<1];
			for(int L=0;L<lim;L+=len<<1){
				int wk=1;
				REP(i,L,L+len-1){
					int u=A[i],v=1ll*A[i+len]*wk%mod;
					A[i]=(u+v)%mod;
					A[i+len]=(u-v)%mod;
					wk=1ll*wk*w%mod;
				}
			}
		}
		if(ty==-1){
			int inv=qpow(lim);
			REP(i,0,lim-1)A[i]=1ll*A[i]*inv%mod;
		}
	}
	void multi(int *A,int *B,int *C,int la,int lb){
		lim=1,cnt=0;
		while(lim<=la+lb)lim<<=1,++cnt;
		REP(i,0,lim-1){
			dn[i]=dn[i>>1]>>1|((i&1)<<(cnt-1));
			if(i>la)A[i]=0;
			if(i>lb)B[i]=0;
		}
		dft(A,1),dft(B,1);
		REP(i,0,lim-1)C[i]=1ll*A[i]*B[i]%mod;
		dft(C,-1);
	}
}

void math_init(){
	fac[0]=1;
	int lim=max(n,m);
	REP(i,1,lim)fac[i]=1ll*fac[i-1]*i%mod;
	ifac[lim]=qpow(fac[lim],mod-2);
	DREP(i,lim-1,0)ifac[i]=1ll*ifac[i+1]*(i+1)%mod;

	using namespace Poly;
	lim=1;
	while(lim<=m+m)lim<<=1;
	om[lim]=qpow(3,(mod-1)/lim);
	for(int i=lim>>1;i;i>>=1)
		om[i]=1ll*om[i<<1]*om[i<<1]%mod;
}

int f[maxm<<2],g[maxm<<2],ap[maxm<<2],bp[maxm<<2];

void divide(int l,int r){
	if(l==r){
		g[l]=(f[l]-1ll*g[l]*ifac[l]%mod)%mod;
		inc(ans,1ll*g[l]*w[l]%mod);
		return;
	}

	using namespace Poly;
	int mid=(l+r)>>1;

	divide(mid+1,r);

	REP(i,0,r-mid-1)ap[i]=1ll*g[i+mid+1]*fac[i+mid+1]%mod;
	REP(i,0,r-l-1)bp[i]=ifac[i+1];
	reverse(bp,bp+r-l);
	multi(ap,bp,ap,r-mid-1,r-l-1);
	REP(i,l,mid)inc(g[i],ap[i-mid-1+r-l]);

	divide(l,mid);
}

int main(){
	//File();

	read(n),read(m),read(S);
	REP(i,0,m)read(w[i]);

	math_init();

	REP(i,0,m)if(S*i<=n)
		f[i]=1ll*C(m,i)*fac[n]%mod*qpow(ifac[S],i)%mod*ifac[n-S*i]%mod*qpow(m-i,n-S*i)%mod;

	divide(0,m);

	printf("%d\n",ans);

	return 0;
}