[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_i\over i!}
\]
令\(F(x)=\sum_{i=0}^{\infty} {f_i\over i!}x^i,G(x)=\sum_{i=0}^{\infty} {1\over 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\n",ans);
}
int main()
{
#ifndef ONLINE_JUDGE
freopen("paint.in","r",stdin);
freopen("paint.out","w",stdout);
#endif
init();
get_f();
doing();
return 0;
}
