YY模拟：绝对伏特加（数学期望）

题意：

$\text{Alan}$ 在玩骰子游戏，$\text{Alan}$ 会玩$n$ 轮骰子，每轮的数值在$[1,K]$ 中随机出现。记$a_i$ 表示$n$ 轮投掷中，数值$i$ 出现的次数，求$a^F_1*a^F_2* ... * a^F_L$ 的期望。答案对$2003$取模。

题解：
记 $x_{i,j}$表示第$i$轮结果是否为$j$，那么：

$$ans = E ((x_{1,1}+...+x_{n,1})^F*...*(x_{1,L}+...+x_{n,L})^F)$$

展开多项式，考虑每个单项式的期望的和。 对于一个特定的单项式 ，是在每一个括号里面取$F$项乘起来得到的，所以共有$F*L$ 项 ，我们把从第$i$个括号中取的第$j$个数$x_{b,i}$的$b$放入$A_{i,j}$中，就得到了一个大小为$F*L$的$A$矩阵。所以每个单项式唯一对应一个$A$矩阵。

对于一个特定的$A$矩阵，如果$\exists_{i,j}满足 A_{i,a} = A_{j,b},i \not = j$，那么这个单项式的期望为0，因为一轮骰子不可能出现两个数。

否则这个单项式的期望为$\frac{1}{K^t}$,$t$为$A_{i,j}$的种类数。
枚举$t$ , 问题变为了我们在$F*L$的矩阵中填入$[1,n]$ ，且不同种类数为$t$的方案数。

按行考虑，我们先不考虑颜色种类和彼此的顺序，最后乘上$\binom{n}{t} t!$即可,设$dp_{i,j}$ 表示一行第$i$位出现$j$种数的方案数，则有$dp_{i,j}=dp_{i-1,j-1}+dp_{i-1,j}*j$。
最后倍增做$L$次背包（多项式乘法），然后枚举第几项即可。 注意对$2003$取模，后面的项因为要乘$t!$被消掉的，所以只需要存$2003$位。

时间复杂度：$O(F^2 + 2003^2 \log L )$

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N=1e3+50, mod=2003;
inline int add(int x,int y) {return (x+y>=mod) ? (x+y-mod) : (x+y);}
inline int dec(int x,int y) {return (x-y<0) ? (x-y+mod) : (x-y);}
inline int mul(int x,int y) {return (LL)x*y%mod;}
inline int power(int a,int b,int rs=1) {for(;b;b>>=1,a=mul(a,a)) if(b&1) rs=mul(rs,a); return rs;}
int n,K,L,F;
int fac[mod],dn[mod],dp[N][N],G[mod],f[mod],ipw[mod];
inline void mul(int *a,int *b,int *c) {
    static int tp[mod],tp2[mod];
    for(int i=0;i<mod;i++) tp[i]=a[i], tp2[i]=b[i], c[i]=0;
    for(int i=0;i<mod;i++)
        for(int j=0;j<=i;j++)
            c[i]=add(c[i],mul(tp[j],tp2[i-j]));
}
int main() {
    cin>>n>>K>>L>>F;
    const int invk=power(K,mod-2);
    ipw[0]=fac[0]=dn[0]=1;
    for(int i=1;i<mod;i++) ipw[i]=mul(ipw[i-1],invk), fac[i]=mul(fac[i-1],i);
    for(int i=1;i<mod;i++) dn[i]=mul(dn[i-1],n-i+1);
    dp[0][0]=1;
    for(int i=1;i<=F;i++)
        for(int j=1;j<=i;j++)
            dp[i][j]=add(mul(dp[i-1][j],j),dp[i-1][j-1]);
    for(int i=1;i<=F;i++) G[i]=dp[F][i];
    f[0]=1;
    for(int b=L;b;b>>=1,mul(G,G,G))
        if(b&1) mul(f,G,f);
    int ans=0;
    for(int i=1;i<mod && i<=n;i++)
        ans=add(ans,mul(dn[i],mul(ipw[i],f[i])));
    cout<<ans;
}