【UOJ449】【集训队作业2018】喂鸽子（Min-Max容斥+NTT优化DP）

点此看题面

• 有\(n\)只鸽子，每秒随机喂一只，求所有鸽子都被喂了至少\(k\)次的期望时间。
• \(n\le50,k\le10^3\)

---

据说正解是一个很神奇的算法，但想想太过神仙没去深究。

不过感觉我这种应该算是比较套路的朴素做法？还是挺容易就能推出来的。

\(Min-Max\)容斥

要求所有鸽子都被喂了至少\(k\)次显然很棘手，立刻想到\(Min-Max\)容斥。

考虑它的公式：

\[E(max(S))=\sum_{T\subseteq S}(-1)^{|T|-1}E(min(T)) \]

然后发现此题中鸽子之间没有差异，所以相同大小的不同鸽子集合其实都是等价的，可以改写成：

\[E(max(n))=\sum_{i=1}^n(-1)^{i-1}C_n^i E(min(i)) \]

那么现在就是要求\(i\)只鸽子中出现一只被喂了\(k\)次的期望时间。

简单期望\(DP\)+无脑\(NTT\)优化

考虑到一个经典公式：

\[E(X)=\sum_{x\ge 0}P(X>x) \]

所以我们只要求出\(i\)只鸽子总共被喂了\(j\)次还没有出现一只被喂了\(k\)次的概率。

注意，我们这里的\(j\)次指的是作用在这\(i\)只鸽子上的有效喂食次数，也就是我们这里算出的期望最后还要除以一次喂食作用在这\(i\)只鸽子上的概率（\(\frac in\)）换算到全局。

众所周知概率等于合法方案数除以总方案数，总方案数显然是\(i^j\)，而合法方案数可以设一个\(f_{i,j}\)表示\(i\)只鸽子总共被喂了\(j\)次还没有出现一只被喂了\(k\)次的方案数来计算。

每次可以直接从\(f_{i-1}\)转移过来：

\[f_{i,j+l}=\sum f_{i-1,j}\times C_{j+l}^l \]

套路地暴拆组合数：

\[\frac{f_{i,j+l}}{(j+l)!}=\sum\frac{f_{i-1,j}}{j!}\times \frac1{l!} \]

显然这是一个卷积的形式，因此我们只要维护好\(f\)的\(EGF\)形式，每次与阶乘逆元的生成函数卷一卷就完成了转移。

然后就结束了。

代码：\(O(nk\log nk)\)

#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define N 50
#define K 1000
#define X 998244353
#define C(x,y) (1LL*Fac[x]*IFac[y]%X*IFac[(x)-(y)]%X)
using namespace std;
int n,k,f[N*K+5];I int QP(RI x,RI y) {RI t=1;W(y) y&1&&(t=1LL*t*x%X),x=1LL*x*x%X,y>>=1;return t;}
int Fac[N*K+5],IFac[N*K+5];I void InitFac(CI S)//预处理阶乘和阶乘逆元
{
	RI i;for(Fac[0]=i=1;i<=S;++i) Fac[i]=1LL*Fac[i-1]*i%X;
	for(IFac[i=S]=QP(Fac[S],X-2);i;--i) IFac[i-1]=1LL*IFac[i]*i%X;
}
namespace Poly//多项式模板
{
	#define PR 3
	int P,L,A[N*K<<2],B[N*K<<2],R[N*K<<2];I void NTT(int* s,CI op)
	{
		RI i,j,k,x,y,U,S;for(i=0;i^P;++i) i<R[i]&&(swap(s[i],s[R[i]]),0);
		for(i=1;i^P;i<<=1) for(U=QP(QP(PR,op),(X-1)/(i<<1)),j=0;j^P;j+=i<<1) for(S=1,
			k=0;k^i;S=1LL*S*U%X,++k) s[j+k]=((x=s[j+k])+(y=1LL*S*s[i+j+k]%X))%X,s[i+j+k]=(x-y+X)%X;
	}
	I void Mul(CI n,int* a,CI m,int* b)//把a,b卷起来
	{
		RI i;P=1,L=0;W(P<=n+m) P<<=1,++L;
		for(i=0;i^P;++i) A[i]=i<=n?a[i]:0,B[i]=i<=m?b[i]:0,R[i]=(R[i>>1]>>1)|((i&1)<<L-1);
		for(NTT(A,1),NTT(B,1),i=0;i^P;++i) A[i]=1LL*A[i]*B[i]%X;
		RI t=QP(P,X-2);for(NTT(A,X-2),i=0;i<=n+m;++i) a[i]=1LL*A[i]*t%X;
	}
}
int main()
{
	RI i,j,l,t,s=0,p,Inv;for(scanf("%d%d",&n,&k),InitFac(n*k),f[0]=i=1;i<=n;++i)
	{
		Poly::Mul((i-1)*(k-1),f,k-1,IFac);//卷上阶乘逆元的生成函数得到当前f的EGF
		for(p=1,Inv=QP(i,X-2),t=j=0;j<=i*(k-1);p=1LL*p*Inv%X,++j) t=(t+1LL*f[j]*Fac[j]%X*p)%X;//共喂了j次还没有一只被喂k次的概率
		t=1LL*t*n%X*Inv%X,s=(s+(i&1?1LL:X-1LL)*C(n,i)%X*t)%X;//除以i/n换算到全局，Min-Max容斥更新答案
	}return printf("%d\n",s),0;
}