传送门
思路还是很巧妙的。
显然能对答案产生贡献的最多颜色为(min{m,n/S})
设(f(i))为出现次数为(S)的颜色种类至少有(i)种的涂色方案
那么
[f(i)=inom{m}{i}frac{n!}{(S!)^i(n-iS)!}(m-i)^{n-iS}
]
然后让我们分析一下这个式子，(inom{m}{i})很显然
重点是(frac{n!}{(S!)^i(n-iS)!})是怎么推出来的
我们发现对于那些恰好出现了(S)次的颜色
首先这(iS)个位置的选法是(inom{n}{iS}=frac{n!}{(n-iS)!(iS)!})
然后对于这(i)种颜色的排布，我们可以列出它的方案数是等于
[frac{(iS)!}{(iS-S)!S!}frac{(iS-S)!}{(iS-2S)!S!}frac{(iS-2S)!}{(iS-3S)!S!}...\
frac{(iS)!}{S!}frac{1}{S!}frac{1}{S!}...\
frac{(iS)!}{(S!)^i}
]
然后就可以得出来一共的方案数为(frac{n!}{(S!)^i(n-iS)!})
然后对于剩下的((n-iS))个位置，剩下的((m-i))种颜色随便放，就是((m-i)^{n-iS})
于是(f(i))就推出来了
然后用一个很简单的容斥计算答案
[ans(i)=sum_{j=i}^{min(m,n/s)}(-1)^{j-i}inom{j}{i}f(j)
]
看到组合数就要拆
[ans(i)=sum_{j=i}^{min(m,n/s)}(-1)^{j-i}frac{j!}{i!(j-i)!}f(j)\
i!*ans(i)=sum_{j=i}^{min(m,n/s)}frac{(-1)^{j-i}}{(j-i)!}f(j)j!\
]
设
[A(x)=frac{(-1)^x}{x!}\
B(x)=f(x)x!
]
那么
[i!*ans(i)=sum_{j=i}^{min(m,n/s)}A(j-i)B(j)\
i!*ans(i)=sum_{j=0}^{min(m,n/s)-i}A(j)B(j+i)\
]
然后显然就是将(B)给reverse一下
[i!*ans(i)=sum_{j=0}^{min(m,n/s)-i}A(j)B(min(m,n/s)-j-i)\
]
显然就是一个卷积的形式了，NTT做完之后再reverse回来
设(D)是做完之后的多项式
[ans=sum_{i=0}^{min(m,n/s)}frac{w(i)D(i)}{i!}
]
代码：
#include<cstdio>
#include<iostream>
#include<algorithm>
using namespace std;
void read(int &x) {
char ch; bool ok;
for(ok=0,ch=getchar(); !isdigit(ch); ch=getchar()) if(ch=='-') ok=1;
for(x=0; isdigit(ch); x=x*10+ch-'0',ch=getchar()); if(ok) x=-x;
}
const int maxn=5e5+10,maxm=1e7+10,mod=1004535809,g=3,gi=334845270;
#define rg register
int tot,r[maxm],len,n,m,ans,s,w[maxm],N,f[maxm],A[maxm],B[maxm],fac[maxm],inv[maxm];
int mul(int x,int y){return 1ll*x*y-1ll*x*y/mod*mod;}
int add(int x,int y){return x+y>=mod?x+y-mod:x+y;}
int del(int x,int y){return x-y<0?x-y+mod:x-y;}
int mi(int a,int b){
int ans=1;
while(b){
if(b&1)ans=mul(ans,a);
b>>=1,a=mul(a,a);
}
return ans;
}
void ntt(int *a,int f){
for(rg int i=0;i<N;i++)if(r[i]>i)swap(a[i],a[r[i]]);
for(rg int i=1;i<N;i<<=1){
int wn=mi(f?g:gi,(mod-1)/(i<<1));
for(rg int j=0;j<N;j+=(i<<1)){
int w=1;
for(rg int k=0;k<i;k++){
int x=a[j+k],y=mul(w,a[j+k+i]);
a[j+k]=add(x,y),a[j+k+i]=del(x,y),w=mul(w,wn);
}
}
}
if(f)return ;int inv=mi(N,mod-2);
for(rg int i=0;i<N;i++)a[i]=mul(a[i],inv);
}
int C(int x,int y){return mul(fac[x],mul(inv[y],inv[x-y]));}
int main()
{
read(n),read(m),read(s);N=min(n/s,m);
fac[0]=inv[0]=1;int now=max(n,m);
for(rg int i=1;i<=now;i++)fac[i]=mul(fac[i-1],i);
inv[now]=mi(fac[now],mod-2);
for(rg int i=now-1;i;i--)inv[i]=mul(inv[i+1],i+1);
for(rg int i=0;i<=m;i++)read(w[i]);
for(rg int i=0;i<=N;i++)
f[i]=mul(C(m,i),mul(fac[n],mul(inv[n-i*s],mul(mi(mi(fac[s],i),mod-2),mi(m-i,n-i*s)))));
for(rg int i=0;i<=N;i++)A[i]=mul((i&1?mod-1:1),inv[i]),B[i]=mul(f[i],fac[i]);
reverse(B,B+N+1);tot=N*2;for(N=1;N<=tot;N<<=1)len++;
for(rg int i=0;i<N;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(len-1));
ntt(A,1),ntt(B,1);
for(rg int i=0;i<N;i++)A[i]=mul(A[i],B[i]);
ntt(A,0);N=tot>>1;reverse(A,A+N+1);
for(rg int i=0;i<=N;i++)ans=add(ans,mul(w[i],mul(A[i],inv[i])));
printf("%d
",ans);
}