[题解 LuoguP4491 [HAOI2018]染色

传送门

神仙计数题 Orz

先令(F[k])表示出现次数恰好为(S)次的颜色恰好有(k)中的方案数，那么

[Ans=sumlimits_{i=0}^mW_iF[i] ]

怎么求(F[k])呢？一个naive的想法，我指定哪(k)种颜色恰好染(S)次，然后剩下的(n-kS)个位置用剩下的(m-k)种颜色随便染

相当与对一个有(k+1)种元素的集合（前(k)中元素分别有(S)个，最后一种元素有(n-kS)个）做可重复集合全排列

这个方案数是(frac{n!}{(S!)^k(n-kS)!})的

剩下的(n-kS)个位置又可以用剩下的颜色随便染，所以还要乘上((m-k)^{n-kS})

所以我们得到了一个柿子(inom{m}{k} imesfrac{n!}{(S!)^k(n-kS)!} imes (m-k)^{n-kS})

从这个柿子也可以看出，合法的颜色种数不会超过(lim=min(m,lfloor n/S floor))

然后发现他假了，首先这样得到的并不是恰好有(k)种的方案数，其次在后面随便染色的时候我们还可能算上重复的方案，它根本不能称为方案数。

但不要放弃希望，我们把上面的柿子叫做(G[k])好了。考虑怎么用(F)算(G)，又有一个naive的想法枚举合法的颜色有多少个得到这个柿子

(G[k]=sum_{i=k}^{lim} F[i])

但上面我们说过了，(G[k])可能会算上重复的方案数。

对于一个恰好有(i)种颜色满足要求的方案，不难发现在(G[k])中会被重复计算(inom{i}{k})次（先在(i)个颜色中指定(k)个然后剩下的随便染，这样可能会染到同一种方案）

所以

[G[k]=sumlimits_{i=k}^{lim} inom{i}{k} F[i] ]

发现右边出现了组合数，这恰恰是一个可以二项式反演的柿子，于是反演一波得到

[F[k]=sumlimits_{i=k}^{lim} (-1)^{i-k}inom{i}{k} G[i] ]

这样就可以得到一个平方的算法了，考虑优化。

先把组合数炸开来

[F[k]=sumlimits_{i=k}^{lim} (-1)^{i-k}frac{i!}{k!(i-k)!}G[i] ]

然后

[F[k] imes k!=sumlimits_{i=k}^{lim}frac{(-1)^{i-k}}{(i-k)!} imes i!G[i] ]

令(A[i]=i!G[i],B[i]=frac{(-1)^i}{i!})，有

[F[k] imes k!=sumlimits_{i=k}^{lim} A[i]B[i-k] ]

这个柿子已经有点可以卷积了，随便翻转(A)或者(B)都可以，下面翻转(A)好了，令翻转(A[0...lim])后得到(A^T)，那么

[F[k] imes k!=sumlimits_{i=k}^{lim} A^T[lim-i]B[i-k] ]

将(A)与(B)卷积，然后(F[k] imes k!)就是卷积的第(lim-k)项。

#include <bits/stdc++.h>
using namespace std;
#define rep(i,a,b) for (int i=(a);i<(b);++i)
#define per(i,a,b) for (int i=(a)-1;i>=(b);--i)
#define mp make_pair
#define pb push_back
#define fi first
#define se second
#define all(x) (x).begin(),(x).end()
#define SZ(x) ((int)(x).size())
typedef double db;
typedef long long ll;
typedef pair<int,int> PII;
typedef vector<int> VI;

const int NN=1e7+10,N=3e5+10,P=1004535809;
inline int add(int x,int y) {return (x+=y)>=P?x-P:x;}
inline int sub(int x,int y) {return (x-=y)<0?x+P:x;}
inline int normal(int x) {return x<0?x+P:x;}
inline int fpow(int x,int y) {
    int ret=1; for(;y;y>>=1,x=1ll*x*x%P)
        if(y&1) ret=1ll*ret*x%P;
    return ret;
}
const int gn=3,ign=fpow(gn,P-2);

int fac[NN],ifac[NN],inv[NN];
inline int getC(int n,int r) {return 1ll*fac[n]*ifac[n-r]%P*ifac[r]%P;}

namespace Poly {
    int rev[N];
    inline void init(int n) {
        rep(i,0,n) rev[i]=rev[i>>1]>>1|((i&1)?n>>1:0);
    }

    void ntt(int *f,int n,int flg) {
        rep(i,0,n) if(rev[i]<i) swap(f[i],f[rev[i]]);
        for(int len=2,k=1;len<=n;len<<=1,k<<=1) {
            int wn=fpow(flg==1?gn:ign,(P-1)/len);
            for(int i=0;i<n;i+=len)
                for(int j=i,w=1;j<i+k;j++,w=1ll*w*wn%P) {
                    int tmp=1ll*w*f[j+k]%P;
                    f[j+k]=sub(f[j],tmp),f[j]=add(f[j],tmp);
                }
        }
        if(flg==-1) {
            int invn=fpow(n,P-2);
            rep(i,0,n+1) f[i]=1ll*f[i]*invn%P;
        }
    }
}
using Poly::ntt;

void init(int n) {
    ifac[0]=ifac[1]=fac[0]=fac[1]=inv[1]=1;
    rep(i,2,n+1) {
        inv[i]=1ll*inv[P%i]*(P-P/i)%P;
        ifac[i]=1ll*ifac[i-1]*inv[i]%P;
        fac[i]=1ll*fac[i-1]*i%P;
    }
}

int A[N],B[N],W[N];

int main() {
#ifdef LOCAL
    freopen("a.in","r",stdin);
#endif
    int n,m,s; scanf("%d%d%d",&n,&m,&s);
    rep(i,0,m+1) scanf("%d",&W[i]);
    init(max(max(n,m),s));
    int lim=min(m,n/s);
    rep(i,0,lim+1) {
        A[i]=1ll*fac[i]*getC(m,i)%P*fac[n]%P*fpow(ifac[s],i)%P
                *ifac[n-i*s]%P*fpow(m-i,n-i*s)%P;
        B[i]=(i&1)?P-ifac[i]:ifac[i];
    }
    reverse(A,A+lim+1);
    int limit=1; while(limit<=2*lim) limit<<=1; Poly::init(limit);
    ntt(A,limit,1),ntt(B,limit,1);
    rep(i,0,limit) A[i]=1ll*A[i]*B[i]%P;
    ntt(A,limit,-1);
    int ans=0;
    rep(i,0,m+1) ans=add(ans,1ll*W[i]*A[lim-i]%P*ifac[i]%P);
    printf("%d
",ans);
    return 0;
}