【CodeChef-EASYEX】Easy exam（斯特林数）（FFT）

传送门

---

题解：

设$x_{i,j}$表示第$i$次是否扔出了$j$，则我们要求的就是$$\prod _{i=1}^L(\sum _{j=1}^n{x}_{j,i}{)}^F$$

的期望。

考虑这个式子完全展开后的某一项，假设有$z$个$x_{j,i}$，则系数为$S_{F,z}z!$，所有的有$z$个项的系数和为$S_{F,z}{n}^{\underline{z}}$

然而考虑到选择$x_{j,i_1},x_{j,i_2}$，是不可能的，贡献为$0$。

设$f_{i,j}$表示在前$i$个$(\sum x)f$中，选出$j$个不同的$x$的方案数（不考虑选的顺序）。

发现转移对应一个子集划分，可以直接用卷积算出来。

支持一个多项式快速幂就行了。

---

代码：

#include<bits/stdc++.h>
#define ll long long
#define re register
#define cs const

using std::cerr;
using std::cout;

cs int mod=2003;
inline int power(int a,int b){
	int res=1;for(;b;b>>=1,a=(ll)a*a%mod)(b&1)&&(res=(ll)a*res%mod);
	return res;
}

cs double PI=acos(-1);
struct cp{
	double x,y;
	cp(){}
	cp(double _x,double _y):x(_x),y(_y){}
	friend cp operator+(cs cp &a,cs cp &b){return cp(a.x+b.x,a.y+b.y);}
	friend cp operator-(cs cp &a,cs cp &b){return cp(a.x-b.x,a.y-b.y);}
	friend cp operator*(cs cp &a,cs cp &b){return cp(a.x*b.x-a.y*b.y,a.x*b.y+b.x*a.y);}
	friend cp operator/(cs cp &a,int n){return cp(a.x/n,a.y/n);}
};

cs int SIZE=1<<15|1;
int r[SIZE];
inline void FFT(cp *A,int len,int typ){
	for(int re i=0;i<len;++i)if(i<r[i])std::swap(A[i],A[r[i]]);
	for(int re i=1;i<len;i<<=1){
		cp wn=cp(cos(PI/i),typ*sin(PI/i));
		for(int re j=0;j<len;j+=i<<1){
			cp w(1,0);
			for(int re k=0;k<i;++k,w=w*wn){
				cp x=A[j+k],y=A[j+k+i]*w;
				A[j+k]=x+y;
				A[j+k+i]=x-y;
			}
		}
	}
	if(typ==-1)for(int re i=0;i<len;++i)A[i]=A[i]/len;
}

inline int mul(int *a,int al,int *b,int bl,int *c){
	static cp A[SIZE],B[SIZE];
	int l=1,cl=al+bl;
	while(l<=cl)l<<=1;
	for(int re i=0;i<l;++i)r[i]=r[i>>1]>>1|((i&1)?l>>1:0);
	for(int re i=0;i<l;++i)A[i]=B[i]=cp(0,0);
	for(int re i=0;i<=al;++i)A[i].x=a[i];
	for(int re i=0;i<=bl;++i)B[i].x=b[i];
	FFT(A,l,1),FFT(B,l,1);
	for(int re i=0;i<l;++i)A[i]=A[i]*B[i];
	FFT(A,l,-1);
	cl=std::min(cl,mod+5);
	for(int re i=0;i<=cl;++i)c[i]=((ll)floor(A[i].x+0.5))%mod;
	return cl;
}

cs int N=1e3+7;
int S[N][N];
int a[SIZE],b[SIZE];

inline void solve(){
	int n,k,l,f;scanf("%d%d%d%d",&n,&k,&l,&f);
	if(k%mod==0){
		puts("0");
		return ;
	}
	for(int re i=0;i<=f;++i)a[i]=b[i]=S[f][i];
	int al=f,bl=f;
	for(int re t=l-1;t;){
		if(t&1)bl=mul(a,al,b,bl,b);
		if(t>>=1)al=mul(a,al,a,al,a);
	}
	int A=1,ik=power(k,mod-2),ans=0;
	for(int re i=1;i<=std::min(n,l);++i)
	A=(ll)A*(n-i+1)%mod*ik%mod;
	for(int re i=l;i<=std::min(n,l*f);++i){
		if(!A)break;
		ans=(ans+(ll)A*b[i])%mod;
		A=(ll)A*(n-i)%mod*ik%mod;
	}
	cout<<ans<<"\n";
}

signed main(){
#ifdef zxyoi
	freopen("easyex.in","r",stdin);
#endif
	S[0][0]=1;
	for(int re i=1;i<N;++i)
	for(int re j=1;j<=i;++j)S[i][j]=(S[i-1][j-1]+j*S[i-1][j])%mod;
	int T;scanf("%d",&T);
	while(T--)solve();
	return 0;
}