[HAOI2018]染色 容斥

对于一个颜色集合$X$，设$X$中所有颜色均出现$S$次（不考虑其他颜色）的方案数为$f_X$，设X中所有颜色均出现$S$次（不存在其他颜色出现S次）的方案数为$g_X$

可以得到：
$$f_X=\frac{N!}{(S!{)}^{|X|}\ast (N-S\cdot |X|)!}\cdot (M-|X|{)}^{N-S\cdot |X|}$$
$$f_X=\sum _{X\subseteq Y}{g}_Y$$
对第二条式子容斥，得到：
$$g_X=\sum _{X\subseteq Y}(-1{)}^{|Y|-|X|}{f}_Y$$
因为$f_x,g_x$只于$|X|$有关，所以：
$$g_X=\sum _{j\ge |X|}(-1{)}^{j-|X|}(\genfrac{}{}{0ex}{}{m-|X|}{j-|X|})f_j$$
$$ans=\sum _{i=0}^m{w}_i\left(\genfrac{}{}{0ex}{}{m}{i}\right)g_i$$

算法瓶颈在于容斥
$$g_i=\sum _{j\ge i}(-1{)}^{j-i}(\genfrac{}{}{0ex}{}{m-i}{j-i})f_j$$
$$=\sum _{j\ge i}(-1{)}^{j-i}(\genfrac{}{}{0ex}{}{m-i}{j-i})f_j$$
$$=\sum _{j\ge i}(-1{)}^{j-i}\frac{(m-i)!}{(j-i)!(m-j)!}{f}_j$$
$$\frac{g_i}{(m-i)!}=\sum _{j\ge i}(-1{)}^{j-i}\cdot \frac{1}{(j-i)!}\frac{f_j}{(m-j)!}$$
$$\frac{g_{m-i}}{i!}=\sum _{j\le i}(-1{)}^{i-j}\cdot \frac{1}{(i-j)!}\frac{f_{m-j}}{j!}$$
NTT即可

#include <bits/stdc++.h>
#define N (1<<18)
#define P 1004535809
#define rint register int
using namespace std;
typedef long long ll;

inline int qpow(int a,int x) {
 int res = 1;
 for (;x;x>>=1) {
  if (x & 1) res = 1ll * res * a % P;
  a = 1ll * a * a % P;
 } return res;
}
inline void dft(int *a,int f) {
 for (rint i=1;i<N;i++) {
  rev[i] = (rev[i>>1] >> 1) | ((i&1) ? (N>>1) : 0);
  if (i > rev[i]) swap(a[i],a[rev[i]]);
 }
 for (rint i=1;i<N;i<<=1) {
  int wn = qpow(3,(P-1)/(i<<1)), w, X, Y;
  if (f == -1) wn = qpow(wn,P-2);
  for (rint j=0;j<N;j+=(i<<1)) {
   w = 1;
   for (rint k=0;k<i;k++,w = 1ll * w * wn % P) {
    X = a[j+k]; Y = 1ll * w * a[j+i+k] % P;
    a[j+k] = (X+Y) % P;  a[j+i+k] = (P+X-Y) % P;
   }
  }
 }
  if (f == -1) {
  int invn = qpow(N,P-2);
  for (rint i=0;i<N;i++) a[i] = 1ll * a[i] * invn % P;
 }
}
int main() {
 for (rint i=1;i<N;i++) rev[i] = (rev[i>>1] >> 1) | ((i&1) << 17);
 scanf("%d%d%d",&n,&m,&s);
 for (rint i=0;i<=m;i++) scanf("%d",w+i);
 fac[0] = fac[1] = inv[0] = inv[1] = 1;  mx = max(n,max(m,s));
 
 for (rint i=2;i<=mx;i++) inv[i] = P - (1ll * (P/i) * inv[P % i] % P);
 for (rint i=2;i<=mx;i++) fac[i] = 1ll * fac[i-1] * i % P, inv[i] = 1ll * inv[i-1] * inv[i] % P;
 t = min(m,n/s);
 for (rint i=0;i<=t;i++) {
  f[i] = 1ll * fac[n] * qpow(inv[s],i) % P * inv[n-i*s] % P * qpow(m-i,n-i*s) % P;
 }
 
 for (rint i=0;i<=m;i++) {
  a[i] = 1ll * f[m-i] * inv[i] % P;
  b[i] = (i & 1) ? P - inv[i] : inv[i];
 }
 dft(a,1);  dft(b,1);
  dft(c,-1);
 
 for (rint i=0;i<=t;i++) {
  ans = (1ll * c[m-i] * fac[m-i] % P * w[i] % P * C(m,i) + ans) % P;
 }
 printf("%d\n",ans);
}