ZJOI模拟 绝对伏特加【数学期望+组合数学+生成函数】

题目描述：

$Alan$在玩骰子游戏，$Alan$会玩$n$轮骰子，每轮的数值在$[1,K]$中随机出现。记$a _i$表示$n$轮投掷中，数值$i$出现的次数，求$a_1^F*a_2^F*……a_L^F$的期望。答案对2003取模。
$1\le n,k\le 10^9,L*F\le 50000$

解题思路：

设$x_{i,j}$表示第$j$轮扔到$i$的概率，那么：

$$ans=E((x_{1,1}+…+x_{1,n})^F*...*(x_{L,1}+...+x_{L,n})^F)$$
我们知道单个$x_{i,j}=\frac{1}{k}$，而且和的期望等于期望的和，但只有在各个事件独立时积的期望才等于期望的积，而式子中$x_{i,a}$和$x_{j,a}$显然不独立，所以不能直接算。

考虑将式子拆开，即

$$ans=\sum E(\Pi_{i=1}^L\Pi_{j=1}^Fx_{i,A_{i,j}}) ，A_{i,j}\in[1,n]$$
注意到如果存在$i\ne j,A_{i,a}=A_{j,b}$，即第$A_{i,a}$轮既扔到$i$，又扔到$j$，这显然是不可能的，所以该项值为0；否则记所有$A_{i,j}$中出现的不同值个数为$t$，那么这$t$个事件独立，该项值为$\frac{1}{k^t}$。

考虑枚举在$[1,L*F]$枚举$t$，计算不存在$i\ne j,A_{i,a}=A_{j,b}$的方案数，那么

$$ans=\sum\limits_{t=1}^{L*F}\binom{n}{t}\frac{t!}{k^t}G(t)$$
$G(t)$就表示已经选出了$t$种值且第一次出现顺序下分配给$L*F$个$A_{i,j}$的方案数。

用更直观的语言描述：有$L$个有编号的盒子，每个盒子有$F$个有编号的格子放球，$L ∗ F$ 个球的编号要在$[1,n]$ 中，不同编号有$t$种，且任意两个盒子不能放有相同编号的球，求方案数。

假设一个盒子分到$i$种编号，由于我们已经枚举了编号出现的位置，所以编号之间没有顺序，相当于要把$F$个格子分配到$i$个集合中，方案数显然是第二类斯特林数$S_F^i$ 。

$L$个盒子一共有$t$种编号，考虑生成函数$g(x)=\sum S_F^ix^i$，那么$G(t)$即为$g(x)^L$中$t$次项系数，记为$g(x)^L[x^t]$。

所以

$$ans=\sum\limits_{t=1}^{L*F}\binom{n}{t}\frac{t!}{k^t}g(x)^L[x^t]=\sum\limits_{t=1}^{L*F}\frac{n!}{(n-t)!k^t}g(x)^L[x^t]$$
注意到一旦$t\ge 2003$，$\frac{n!}{(n-t)!}$就一定为0，所以只用考虑2003项即可。

预处理斯特林数，暴力做多项式乘法，时间复杂度为$O(2003F+2003^2logL)$

#include<bits/stdc++.h>
using namespace std;
const int N=50005,mod=2003;
int n,K,L,F,s[N][2500];
int Pow(int x,int y)
{
    int res=1;x%=mod;
    for(;y;y>>=1,x=x*x%mod)
        if(y&1)res=res*x%mod;
    return res;
}
struct Poly
{
    int deg,a[N];
    Poly(){deg=0;memset(a,0,sizeof(a));}
    inline friend Poly mul(const Poly &A,const Poly &B)
    {
        Poly res;res.deg=min(mod,A.deg+B.deg);
        for(int i=0;i<=A.deg;i++)
            for(int j=0;j<=B.deg&&i+j<=mod;j++)
                res.a[i+j]=(res.a[i+j]+A.a[i]*B.a[j])%mod;
        return res;
    }
    inline friend Poly Pow(Poly A,int b)
    {
        Poly res;res.a[0]=1;
        for(;b;b>>=1,A=mul(A,A))
            if(b&1)res=mul(res,A);
        return res;
    }
}G;
void pre()
{
    for(int i=1;i<=F;i++)
    {
        s[i][1]=1;
        for(int j=2;j<=min(i,mod);j++)
            s[i][j]=(s[i-1][j-1]+s[i-1][j]*j)%mod;
    }
    G.deg=F;
    for(int i=1;i<=mod;i++)G.a[i]=s[F][i];
    G=Pow(G,L);
}
int main()
{
    freopen("vodka.in","r",stdin);
    freopen("vodka.out","w",stdout);
    scanf("%d%d%d%d",&n,&K,&L,&F);n%=mod,K%=mod;
    pre();
    int inv=Pow(K,mod-2),num=1,ans=0;
    for(int i=1;i<=L*F&&num;i++)
    {
        num=num*(n-i+1)%mod*inv%mod;
        ans=(ans+num*G.a[i])%mod;
    }
    cout<<ans<<'\n';
    return 0;
}