JZOJ 5680 绝对伏特加

Description

进行$n$轮随机数生成，每轮的数值在[$1$,$K$]中随机出现。记$a_i$表示$n$轮投掷中，数值$i$出现的次数，求$(\Pi_{i=1}^La_i)^F$的期望。答案对$2003$取模。

Solution

$0≤n,K≤10^9$ ,$0<F≤10^3$,$0<L∗F≤50000$,$L≤K$

Solution

先把每个数被选中的概率都看成1，最后答案乘上$1\over {K^n}$即可。
考虑$1$~$L$的贡献的生成函数$F(x)$，显然他们的生成函数都一样，考虑$EGF$

$$F(x)=\sum_{i\geq 0}\frac{i^Fx^i}{i!}$$

考虑用斯特林数把$i^F$拆开，得

$$F(x)=\sum_{i\geq 0}{x^i\over i!}\sum_{j=1}^F\left\{\begin{array}{c}{F}\\{j}\end{array}\right\}i^{\underline j}$$

$$F(x)=\sum_{j=1}^F\left\{\begin{array}{c}{F}\\{j}\end{array}\right\}\sum_{i\geq j}\frac{x^i}{(i-j)!}$$

$$F(x)=\sum_{j=1}^F\left\{\begin{array}{c}{F}\\{j}\end{array}\right\}\sum_{i\geq 0}{x^{i+j}\over i!}$$

$$F(x)=\sum_{j=1}^F\left\{\begin{array}{c}{F}\\{j}\end{array}\right\}x^j\sum_{i\geq 0}{x^{i}\over i!}$$

$$F(x)=\sum_{j=1}^F\left\{\begin{array}{c}{F}\\{j}\end{array}\right\}x^je^x$$

$$F(x)=e^x\sum_{j=1}^F\left\{\begin{array}{c}{F}\\{j}\end{array}\right\}x^j$$

于是乎我们变得出了$1$~$L$每一个数贡献的生成函数$F(x)$，则$1$~$L$所有数贡献的生成函数即为$F^L(x)$。

接下来再考虑剩下的$K-L$个数的贡献的生成函数，由于剩下的数不会对答案有贡献，只会让方案数有所增加，因此可以很容易得出剩下的数的生成函数$G(x)$

$$G(x)=\sum_{i\geq 0}\frac{(K-L)^ix^i}{i!}=\sum_{i\geq 0}\frac{[(K-L)x]^i}{i!}=e^{K-L}$$

所以答案就是让我们求$[x^n]$$n!F^L(x)G(x)$ $mod \ 2003$。

$$n!F^L(x)G(x)=n!e^K(\sum_{j=1}^F\left\{\begin{array}{c}{F}\\{j}\end{array}\right\}x^j)^L$$

考虑$[x^i]n!e^K\ mod\ 2003$的值
$[x^i]n!e^K=n^{\underline {n-i}}K^i$
观察上式，当$i<n-2003$时$[x^i]n!e^K$一定为$2003$的倍数，换言之，就是$n!e^K$中只有$n-2003$至第$n$项是有用的，同时也就意味着$(\sum_{j=1}^F\left\{\begin{array}{c}{F}\\{j}\end{array}\right\}x^j)^L$只有前$2003$项是有用的，套上个快速幂暴力求即可。

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

Code

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>

#define fo(i,j,l) for(int i=j;i<=l;++i)
#define fd(i,j,l) for(int i=j;i>=l;--i)

using namespace std;
typedef long long ll;
const ll N=2300,mo=2003;

int S[2][N];
int f[N],g[N];
ll n,K,L,F;

inline ll ksm(int o,ll t)
{
    int y=1;t=t%(mo-1);
    for(;t;t>>=1,o=o*o%mo)
    if(t&1)y=y*o%mo;
    return y;
}

inline int min(int a,int b)
{return a<b?a:b;} 

inline void zj()
{
    fd(i,mo,0){
        int op=0;
        fo(l,0,i)op=(op+g[l]*g[i-l])%mo;
        g[i]=op;
    }
}

inline void lj()
{
    fd(i,mo,0){
        int op=0;
        fo(l,0,i)op=(op+f[l]*g[i-l])%mo;
        f[i]=op;
    }
}

int main()
{
    cin>>n>>K>>L>>F;
    int v=0,u=0;
    S[0][0]=1;
    fo(i,1,F){
        v=u^1;
        S[v][0]=0;
        fo(l,1,min(i,mo))S[v][l]=(S[u][l-1]+S[u][l]*l)%mo;
        if(i<=mo)S[v][i]=1;
        u=v;
    }
    f[0]=1; fo(i,0,2003)g[i]=S[u][i];
    int r=ksm(K%mo,n),ny=ksm(K%mo,mo-2);
    for(;L;L>>=1,zj())if(L&1)lj();
    int op=1,ans=0;
    fo(i,0,mo){
        ans=(ans+f[i]*r%mo*op)%mo;
        r=r*ny%mo;
        op=op*((n-i)%mo)%mo;
    }
    int yy=ksm(K%mo,mo-2);
    ans=ans*ksm(yy,n)%mo;
    cout<<ans;
}