BZOJ5306 HAOI2018染色（容斥原理+NTT）

容易想到枚举恰好出现S次的颜色有几种。如果固定至少有i种恰好出现S次，那么方案数是C(M,i)·C(N,i*S)·(M-i)^N-i*S·(i*S)!/(S!)ⁱ，设为f(i)。

　　于是考虑容斥，可得恰好i种的答案为Σ(-1)^j-iC(j,i)·f(j) (j=i~min(M,⌊N/S⌋))。因为容斥是一个枚举子集的过程，在算至少i种的方案时，f(j)被计入了C(j,i)次。

　　f显然可以通过预处理阶乘及其逆元线性地算出来。考虑怎么快速算后一部分。注意到模数，NTT没跑了。拆开组合数，可以发现是与j-i有关的式子和与j有关的式子相乘，那么把其中一个翻转一下就是卷积了。

　　容斥好难啊。

#include<iostream> 
#include<cstdio>
#include<cmath>
#include<cstdlib>
#include<cstring>
#include<algorithm>
using namespace std;
int read()
{
    int x=0,f=1;char c=getchar();
    while (c<'0'||c>'9') {if (c=='-') f=-1;c=getchar();}
    while (c>='0'&&c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar();
    return x*f;
}
#define P 1004535809
#define N 10000010
#define M 100010
#define inv3 334845270
int n,m,s,k,t,w[N],f[M*3],a[M*3],fac[N],inv[N],r[M*3],ans=0;
int ksm(int a,int k)
{
    if (k==0) return 1;
    int tmp=ksm(a,k>>1);
    if (k&1) return 1ll*tmp*tmp%P*a%P;
    else return 1ll*tmp*tmp%P;
}
int C(int n,int m){return 1ll*fac[n]*inv[m]%P*inv[n-m]%P;}
void DFT(int n,int *a,int p)
{
    for (int i=0;i<n;i++) if (i<r[i]) swap(a[i],a[r[i]]);
    for (int i=2;i<=n;i<<=1)
    {
        int wn=ksm(p,(P-1)/i);
        for (int j=0;j<n;j+=i)
        {
            int w=1;
            for (int k=j;k<j+(i>>1);k++,w=1ll*w*wn%P)
            {
                int x=a[k],y=1ll*w*a[k+(i>>1)]%P;
                a[k]=(x+y)%P,a[k+(i>>1)]=(x-y+P)%P;
            }
        }
    }
}
void mul(int n,int *a,int *b)
{
    DFT(n,a,3),DFT(n,b,3);
    for (int i=0;i<n;i++) a[i]=1ll*a[i]*b[i]%P;
    DFT(n,a,inv3);
    int inv=ksm(n,P-2);
    for (int i=0;i<n;i++) a[i]=1ll*a[i]*inv%P;
}
int main()
{
#ifndef ONLINE_JUDGE
    freopen("bzoj5306.in","r",stdin);
    freopen("bzoj5306.out","w",stdout);
    const char LL[]="%I64d";
#else
    const char LL[]="%lld";
#endif
    n=read(),m=read(),s=read(),k=min(m,n/s);
    for (int i=0;i<=m;i++) w[i]=read();
    fac[0]=1;for (int i=1;i<=max(n,m);i++) fac[i]=1ll*fac[i-1]*i%P;
    inv[0]=inv[1]=1;for (int i=2;i<=max(n,m);i++) inv[i]=(P-1ll*(P/i)*inv[P%i]%P)%P;
    for (int i=2;i<=max(n,m);i++) inv[i]=1ll*inv[i]*inv[i-1]%P;
    for (int i=0;i<=k;i++)
    f[i]=1ll*C(m,i)*C(n,i*s)%P*ksm(m-i,n-i*s)%P*fac[i*s]%P*ksm(inv[s],i)%P;
    t=1;while (t<=k*2) t<<=1;
    for (int i=0;i<t;i++) r[i]=(r[i>>1]>>1)|(i&1)*(t>>1);
    for (int i=0;i<=k;i++) f[i]=1ll*f[i]*fac[i]%P;
    for (int i=0;i<=k;i++) a[i]=1ll*((i&1)?P-1:1)*inv[i]%P;
    reverse(a,a+k+1);
    mul(t,f,a);
    for (int i=0;i<=k;i++) ans=(ans+1ll*f[i+k]*w[i]%P*inv[i]%P)%P;
    cout<<ans;
    return 0;
}