- 有(n)只鸽子,每秒随机喂一只,求所有鸽子都被喂了至少(k)次的期望时间。
- (nle50,kle10^3)
据说正解是一个很神奇的算法,但想想太过神仙没去深究。
不过感觉我这种应该算是比较套路的朴素做法?还是挺容易就能推出来的。
(Min-Max)容斥
要求所有鸽子都被喂了至少(k)次显然很棘手,立刻想到(Min-Max)容斥。
考虑它的公式:
[E(max(S))=sum_{Tsubseteq 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_{xge 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} imes C_{j+l}^l
]
套路地暴拆组合数:
[frac{f_{i,j+l}}{(j+l)!}=sumfrac{f_{i-1,j}}{j!} imes frac1{l!}
]
显然这是一个卷积的形式,因此我们只要维护好(f)的(EGF)形式,每次与阶乘逆元的生成函数卷一卷就完成了转移。
然后就结束了。
代码:(O(nklog 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
",s),0;
}