[HAOI2018] 染色

对一个长度为 (n) 的序列进行染色，有 (m) 种颜色。对一种方案，如果恰好出现 (s) 次的颜色总数为 (k)，则得分为 (W_k)，求所有染色方案得分的总和。(nleq 10^7,mleq 10^5,s leq 150)

Solution

最大有效颜色数为 (lim=min(m,n/s))

设恰好出现 (s) 次的颜色有至少 (i) 种的方案数位 (f[i])，则选出这 (i) 种颜色，并给他们分配位置，剩下的相互独立填入即可，即

[f[i]=C_m^i frac{n!}{(S!)^i (n-iS)!}(m-i)^{n-iS} ]

设 (ans[i]) 表示出现 (s) 次的颜色恰好有 (i) 种的方案数，由容斥原理，

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

为了方便做卷积，变形为

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

不妨令

[A[k]=frac{(-1)^k}{k!} quad quad B[k]=f[k]cdot k! ]

那么卷积就可以描述为

[C[i]=sum_{j-k=i} A[k]B[j] ]

很自然地，设 (D[k]=A[lim-k])，则

[C[i]=sum_{j-k=i} D[lim-k]B[j]=sum_{j+u=lim+i} D[u]B[j]=sum_{j+k=lim+i}B[j]D[k] ]

答案就是

[sum_{i=0}^{lim} ans[i]cdot W_i ]

#include <bits/stdc++.h>
using namespace std;
#define pw(n) (1<<n)
using namespace std;
#define int long long
namespace NTT {
    const int N=4000005;
    const int mod=1004535809,g=3;
    int n,m,bit,bitnum,a[N+5],b[N+5],rev[N+5];
    void getrev(int l){
        for(int i=0;i<pw(l);i++){
            rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
        }
    }
    int fastpow(int a,int b){
        int ans=1;
        for(;b;b>>=1,a=1LL*a*a%mod){
            if(b&1)ans=1LL*ans*a%mod;
        }
        return ans;
    }
    void NTT(int *s,int op){
        for(int i=0;i<bit;i++)if(i<rev[i])swap(s[i],s[rev[i]]);
        for(int i=1;i<bit;i<<=1){
            int w=fastpow(g,(mod-1)/(i<<1));
            for(int p=i<<1,j=0;j<bit;j+=p){
                int wk=1;
                for(int k=j;k<i+j;k++,wk=1LL*wk*w%mod){
                    int x=s[k],y=1LL*s[k+i]*wk%mod;
                    s[k]=(x+y)%mod;
                    s[k+i]=(x-y+mod)%mod;
                }
            }
        }
        if(op==-1){
            reverse(s+1,s+bit);
            int inv=fastpow(bit,mod-2);
            for(int i=0;i<bit;i++)a[i]=1LL*a[i]*inv%mod;
        }
    }
    void solve(vector <int> A,vector <int> B,vector <int> &C) {
        n=A.size()-1;
        m=B.size()-1;
        for(int i=0;i<=n;i++) a[i]=A[i];
        for(int i=0;i<=m;i++) b[i]=B[i];
        m+=n;
        bitnum=0;
        for(bit=1;bit<=m;bit<<=1)bitnum++;
        getrev(bitnum);
        NTT(a,1);
        NTT(b,1);
        for(int i=0;i<bit;i++)a[i]=1LL*a[i]*b[i]%mod;
        NTT(a,-1);
        C.clear();
        for(int i=0;i<=m;i++) C.push_back(a[i]);
    }
}

struct poly {
    const int mod=1004535809;
    vector <int> a;
    void cut(int n) {
        while(a.size()>n) a.pop_back();
    }
    void load(int *x,int n) {
        a.clear();
        for(int i=0;i<n;i++) a.push_back(x[i]);
    }
    poly operator *(int b) {
        poly c=*this;
        for(int i=0;i<a.size();i++) (((c.a[i]*=b)%=mod)+=mod)%=mod;
        return c;
    }
    poly operator *(const poly &b) {
        poly c;
        NTT::solve(a,b.a,c.a);
        return c;
    }
    poly operator +(poly b) {
        int len=max(a.size(),b.a.size());
        a.resize(len);
        b.a.resize(len);
        poly c;
        for(int i=0;i<len;i++) c.a.push_back((a[i]+b.a[i])%mod);
        return c;
    }
    poly operator -(poly b) {
        int len=max(a.size(),b.a.size());
        a.resize(len);
        b.a.resize(len);
        poly c;
        for(int i=0;i<len;i++) c.a.push_back(((a[i]-b.a[i])%mod+mod)%mod);
        return c;
    }
};

const int N = 10000005;
const int M = 2000005;
const int mod = 1004535809;
int n,m,s,w[M],frac[N],D[M],B[M],C[M],lim;

int qpow(int p,int q) {
    return (q&1?p:1) * (q?qpow(p*p%mod,q/2):1) %mod;
}
int inv(int p) {
    return qpow(p,mod-2);
}
void init_frac() {
    frac[0]=1;
    for(int i=1;i<N;i++) frac[i]=frac[i-1]*i%mod;
}
int c(int n,int m) {
    return frac[m]*inv(frac[n])%mod*inv(frac[m-n])%mod;
}

signed main() {
    ios::sync_with_stdio(false);
    init_frac();
    cin>>n>>m>>s;
    lim=min(m,n/s);
    for(int i=0;i<=m;i++) cin>>w[i];
    for(int i=0;i<=lim;i++) {
        D[i]=(((lim-i)&1?-1:1)*inv(frac[lim-i])+mod)%mod;
        B[i]=c(i,m)*frac[n]%mod*qpow(m-i,n-i*s)%mod*frac[i]%mod*
            inv(qpow(frac[s],i))%mod*inv(frac[n-i*s])%mod;
    }
    poly pB,pD;
    pB.load(B,lim+1);
    pD.load(D,lim+1);
    poly pC=pB*pD;
    int ans=0;
    for(int i=0;i<=lim;i++) ans+=pC.a[i+lim]*inv(frac[i])%mod*w[i]%mod, ans+=mod, ans%=mod;
    cout<<ans;
}