洛咕 P4491 [HAOI2018]染色

显然颜色数量不会超过(lim=min(m,n/S))

考虑容斥，计算恰好出现了(S)次的颜色有至少(i)种的方案数(f[i])，钦定(i)种颜色正好放(S)种

有(m)种颜色选(i)种，所以乘一个(C_m^i)

然后这n个位置分成(i+1)个部分：被钦定的(i)种颜色，每个有(S)个；剩下的(m-i)种颜色，一共(n-iS)个。先看作是可重的全排列数，那么方案就有(frac{n!}{(S!)^i(n-iS)!})种。前(i)各部分都是只有一种颜色，后面部分每个有(m-i)种取法，所以还有一个((m-i)^{n-iS})

综上，(f[i]=C_m^icdot frac{n!}{(S!)^i(n-iS)!}cdot(m-i)^{n-iS})

接下来就是答案，恰好出现了(S)次的颜色有正好(i)种的方案数(ans[i])

用容斥，(ans[i]=sum_{j=i}^{lim}(-1)^{j-i}C_j^if[j])

那个组合数很麻烦，拆开

(ans[i]=sum_{j=i}^{lim}(-1)^{j-i}frac{j!}{i!(j-i)!}f[j])

(ans[i]cdot i!=sum_{j=i}^{lim}frac{(-1)^{j-i}}{(j-i)!}f[j]cdot j!)

这就可以直接用NTT做了，如果不知道怎么做的可以先写zjoi2014 力

#include<bits/stdc++.h>
#define il inline
#define vd void
#define mod 1004535809
typedef long long ll;
il int gi(){
    int x=0,f=1;
    char ch=getchar();
    while(!isdigit(ch)){
        if(ch=='-')f=-1;
        ch=getchar();
    }
    while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
    return x*f;
}
il int pow(int x,int y){
    int ret=1;
    while(y){
        if(y&1)ret=1ll*ret*x%mod;
        x=1ll*x*x%mod;y>>=1;
    }
    return ret;
}
#define inv(a) pow((a),mod-2)
const int G=3,iG=inv(G);
int fact[10000001],W[100010];
int A[262147],B[262147],rev[262147];
il int C(int n,int m){
    if(n<m)return 0;
    return 1ll*fact[n]*inv(1ll*fact[m]*fact[n-m]%mod)%mod;
}
il vd ntt(int*A,int n,int t){
    for(int i=0;i<n;++i)if(rev[i]>i)std::swap(A[rev[i]],A[i]);
    for(int o=1;o<n;o<<=1){
        int W=pow(t?G:iG,(mod-1)/(o<<1));
        for(int*p=A;p!=A+n;p+=o<<1)
            for(int i=0,w=1;i<o;++i,w=1ll*w*W%mod){
                int t=1ll*w*p[i+o]%mod;
                p[i+o]=(p[i]-t+mod)%mod,p[i]=(p[i]+t)%mod;
            }
    }
    if(!t){
        int invN=inv(n);
        for(int i=0;i<n;++i)A[i]=1ll*invN*A[i]%mod;
    }
}
int main(){
    int n=gi(),m=gi(),s=gi();
    for(int i=0;i<=m;++i)W[i]=gi();
    int LIM=std::max(m,n);
    fact[0]=1;for(int i=1;i<=LIM;++i)fact[i]=1ll*fact[i-1]*i%mod;
    int lim=std::min(m,n/s);
    int N=1,lg=0;while(N<(lim+1)<<1)N<<=1,++lg;
    for(int i=0;i<=lim;++i)A[i]=1ll*fact[i]*C(m,i)%mod*fact[n]%mod*pow(m-i,n-i*s)%mod*inv(1ll*pow(fact[s],i)*fact[n-i*s]%mod)%mod;
    for(int i=0;i<=lim;++i){
        B[i]=inv(fact[lim-i]);
        if((lim-i)&1)B[i]=mod-B[i];
    }
    for(int i=0;i<N;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<lg-1);
    ntt(A,N,1),ntt(B,N,1);
    for(int i=0;i<N;++i)A[i]=1ll*A[i]*B[i]%mod;
    ntt(A,N,0);
    int ans=0;
    for(int i=0;i<=lim;++i)ans=(ans+1ll*W[i]*A[lim+i]%mod*inv(fact[i])%mod)%mod;
    printf("%d
",ans);
    return 0;
}