Solution

别人讲的很清楚怎么 dp 了。这里讨论一下为什么这样容斥。

就拿 xzz 给的式子来说吧：$f[i]=C_m^i\cdot \frac{n!}{(S!)^i(n-iS)!}\cdot(m-i)^{n-iS}$

xzz 是这么定义状态的：$f[i]$ 是至少有 $i$ 种颜色出现的方案数。那么就很有趣了，按照上面给的式子，这个会计算重复。比如你钦定 1 这个颜色出现了 S 次，得到状态1122，然后你钦定 2 这个颜色也出现了 S 次，就会有状态计算重复。

然后更有趣的是下一个状态：$ans[i]=\sum_{j=i}^{lim}(-1)^{j-i}C_j^if[j]$

这个恰好 $i$ 次的状态定的没错。

我弄了很久觉得没什么感性的方法可以理解这个方程，于是就画了一个三个集合的文氏图。

假设有 ABC 三种颜色，我们要求出现恰好一次的颜色的总方案数，那么我们在集合图上对应着就是三个集合 没有任何交集的并集部分 。

那么我们实际上会做如下操作：首先三个集合的计数分别 $+\binom{1}{1}$，那么现在就是两个集合的交会被计算 2 次，然后三个集合的交被计算 3 次。

之后两个集合的交计数分别 $-\binom{2}{1}$，那么现在两个集合的交皆被计算 0 次，然后三个集合的交被计算 - 3 次。

之后在给三个集合的交就是 $+\binom{3}{1}$。很神奇的得到了我们要的答案。

考虑一个 $j$ 个集合的交，它会被是它子集的包含 $k$ 个集合的交的每次整体操作重复贡献 $\binom{j}{k}$ 次。那么事实上如果我们要求 $i$ 个集合的交的单独部分，则以此方案下我们需要证明：

$$\sum \limits_{k=i}^j (-1)^{k-i}\binom{k}{i} \binom{j}{k}=[j=i]$$

而这个分类讨论一下奇偶性就很好证明了，即我们可以证明如果这么做确实能得到这样的方案。但是至于这么做的组合意义是什么我还没有理解...

然后就随便 NTT 搞下就行了...

话说为什么我的代码每次常数都那么大...

Code

#include<bits/stdc++.h>
#define MOD (1004535809)
using namespace std;

const int N=1e5+10, V=1e7+10;
typedef int ll;
int fac[V], inv[V], finv[V], f[N<<2], g[N<<2], r[N<<2];

ll qpow(ll x, ll y) {
    x%=MOD; ll ret=1;
    while(y) {
        if(y&1) ret=1ll*ret*x%MOD;
        x=1ll*x*x%MOD, y>>=1;
    }
    return ret;
}

int n;
void build(int m) {
    int l=0;
    for(n=1; n<=m; n<<=1) l++;
    for(int i=1; i<n; i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
}
int pl[]={1, -1};
ll C(int x, int y) { return 1ll*fac[x]*finv[y]%MOD*finv[x-y]%MOD; }
void dft(ll x[], int d) {
    for(int i=1; i<n; i++) if(r[i]>i) swap(x[i], x[r[i]]);
    ll t, w, o;
    for(int i=1; i<n; i<<=1) {
        o=qpow(3, d*(MOD-1)/(i<<1)+MOD-1);
        for(int j=0; j<n; j+=(i<<1)) {
            w=1;
            for(int k=0; k<i; k++, w=1ll*w*o%MOD)
                t=1ll*x[i+j+k]*w%MOD, x[i+j+k]=(x[j+k]-t+MOD)%MOD,
                (x[j+k]+=t)%=MOD;
        }
    }
}

int main() {
    ios::sync_with_stdio(false), cin.tie(0);
    fac[0]=finv[0]=fac[1]=finv[1]=inv[1]=1;
    for(int i=2;i<V;i++) fac[i]=1ll*fac[i-1]*i%MOD, inv[i]=(-1ll*MOD/i*inv[MOD%i]%MOD+MOD)%MOD;
    for(int i=2;i<V;i++) finv[i]=1ll*inv[i]*finv[i-1]%MOD;
    int n, m, s;
    cin>>n>>m>>s;
    for(int i=0; i<=m && n-i*s>=0; i++) {
        f[i]=(pl[i&1]*finv[i]+MOD)%MOD,
        g[i]=1ll*C(m, i)*qpow(finv[s], i)%MOD*fac[n]%MOD*finv[n-i*s]%MOD*qpow(m-i, n-i*s)%MOD*fac[i]%MOD;
    }
    build(m<<1), reverse(g, g+m+1);
    dft(f, 1), dft(g, 1); for(int i=0; i<::n; i++) f[i]=1ll*f[i]*g[i]%MOD;
    dft(f, -1); ll inv=qpow(::n, MOD-2); for(int i=m; i>=0; i--) f[i]=1ll*f[i]*inv%MOD;
    ll w, ans=0;
    for(int i=0;i<=m && n>=1ll*s*i;i++) {
        cin>>w;
        (ans+=1ll*w*f[m-i]%MOD*finv[i]%MOD)%=MOD;
    }
    cout<<ans<<endl;
    return 0;
}