喂鸽子
有 \(n\) 只鸽子，每只鸽子吃 \(k\) 粒玉米就能吃饱。
每秒钟，你会等概率随机选择一只鸽子，并投喂一粒玉米。
问所有鸽子都吃饱的期望时间，对 \(998244353\) 取模。
\(1 \leq n \leq 50\)，\(1 \leq k \leq 1000\)。
考虑使用 min-max 容斥，令 \(S = \{1,2,\dots,n\}\)，则答案为：
\[\sum\limits_{\varnothing \subsetneq T \subseteq S} (-1)^{|T|-1} f(T)
\]
其中 \(f(T)\) 代表喂饱 \(T\) 中任意一只鸽子，期望投喂的玉米粒数。
我们考虑如何计算 \(f(T)\)，容易发现所有鸽子本质相同，所以 \(f(T)\) 只与 \(|T|\) 有关。
我们记 \(f_i\) 为喂饱 \(1\) 到 \(i\) 中任意一只鸽子，期望投喂的玉米粒数，则答案为：
\[\sum\limits_{i=1}^{n} (-1)^{i-1} \dbinom{n}{i} f_i
\]
我们现在只需要求出 \(f_i\) 即可。
考虑有一只鸽子在第 \(s\) 秒吃饱，等价于所有鸽子在第 \(0\) 秒到第 \(s-1\) 秒都没有吃饱。
期望转概率，则 \(f_i = \sum\limits_{j \geq 0} g_{i,j}\)，\(g_{i,j}\) 表示前 \(i\) 只鸽子喂了 \(j\) 粒玉米都没有被喂饱的概率。
我们设 \(h_{i,j}\) 表示只对前 \(i\) 只鸽子投喂，前 \(i\) 只鸽子喂了 \(j\) 粒玉米都没有被喂饱的概率。
此时我们将其代入到上面，可得：
\[f_i = \sum\limits_{j \geq 0} \sum\limits_{t=0}^{j} h_{i,t} \dbinom{j}{t} (\dfrac{i}{n})^{t} (\dfrac{n-i}{n})^{j-t}
\]
我们发现 \(h_{i,t} > 0\) 一定需要满足 \(t \leq i \times (k-1)\)，交换求和顺序：
\[f_i = \sum\limits_{t=0}^{i \times (k-1)} h_{i,t} (\dfrac{i}{n})^{t} \sum\limits_{j \geq t} \dbinom{j}{t} (\dfrac{n-i}{n})^{j-t}
\]
\(j \geq t\) 很不好，考虑将 \(j\) 减去 \(t\)。
\[f_i = \sum\limits_{t=0}^{i \times (k-1)} h_{i,t} (\dfrac{i}{n})^{t} \sum\limits_{j} \dbinom{j + t}{t} (\dfrac{n-i}{n})^{j}
\]
考虑 \(\sum\limits_{j} \dbinom{j + t}{t} (\dfrac{n-i}{n})^{j}\) 这一部分，我们令 \(x = \dfrac{n-i}{n}\)，则原式可化为 \(\sum\limits_{j} \dbinom{j + t}{t} x^j\)。
有结论：\(\sum\limits_{j} \dbinom{j + t}{t} x^j = \dfrac{1}{(1-x)^{t+1}}\)。
也就是说，我们得到了：
\[f_i = \sum\limits_{t=0}^{i \times (k-1)} \dfrac{h_{i,t} (\frac{i}{n})^{t}}{(1-\frac{n-i}{n})^{t+1}}
\]
你发现这很蠢，因为 \(1 - \dfrac{n-i}{n} = \dfrac{i}{n}\)，也就是说：
\[f_i = \sum\limits_{t=0}^{i \times (k-1)} h_{i,t} \times \dfrac{n}{i}
\]
即：
\[f_i \times \dfrac{i}{n} = \sum\limits_{t=0}^{i \times (k-1)} h_{i,t}
\]
也就是说，我们求出 \(h_{i,t}\) 即可。
\(h_{i,t}\) 是概率，很难计算，不妨转成方案数。
我们令 \(p_{i,t}\) 表示只对前 \(i\) 只鸽子投喂，前 \(i\) 只鸽子喂了 \(t\) 粒玉米都没有被喂饱的方案数，显然 \(h_{i,t} = \dfrac{p_{i,t}}{i^t}\)。
考虑求出 \(p_{i,t}\)，显然可以递推求解，使用背包的思想：
\[p_{i,t} = \sum\limits_{j=0}^{\min(t,k-1)} \dbinom{t}{j} p_{i-1,t-j}
\]
为了方便处理，我们做如下变换：
\[p_{i,t} = \sum\limits_{j=0}^{min(t,k-1)} \dfrac{t!}{j! \times (t-j)!} \times p_{i-1,t-j}
\]
也就是：
\[\dfrac{p_{i,t}}{t!} = \sum\limits_{j=0}^{min(t,k-1)} \dfrac{1}{j!} \times \dfrac{p_{i-1,t-j}}{(t-j)!}
\]
使用 NTT 转移即可。
复杂度 \(O(n^2 k \log nk)\)，可以通过。
#include<iostream>
#include<cstdio>
using namespace std;
const int mod=998244353,root=15311432,inv_2=499122177;
const int L=23,N=1<<16;
int omega[1<<L],inv_omega[1<<L];
void init(){
omega[0]=1;
for(int i=1;i<(1<<L);i++){
omega[i]=(long long)omega[i-1]*root%mod;
}
inv_omega[0]=1;
for(int i=1;i<(1<<L);i++){
inv_omega[i]=omega[(1<<L)-i];
}
}
int t,buf[N];
void NTT(int* a){
for(int i=t;i>=1;i--){
for(int j=0;j<(1<<t);j+=(1<<i)){
for(int k=0;k<(1<<i-1);k++){
buf[j+k]=a[j+k*2];
buf[j+k+(1<<i-1)]=a[j+k*2+1];
}
for(int k=0;k<(1<<i);k++){
a[j+k]=buf[j+k];
}
}
}
for(int i=1;i<=t;i++){
for(int j=0;j<(1<<t);j+=(1<<i)){
for(int k=0;k<(1<<i-1);k++){
buf[j+k]=(a[j+k]+(long long)omega[k<<L-i]*a[j+k+(1<<i-1)]%mod+mod)%mod;
buf[j+k+(1<<i-1)]=(a[j+k]-(long long)omega[k<<L-i]*a[j+k+(1<<i-1)]%mod+mod)%mod;
}
for(int k=0;k<(1<<i);k++){
a[j+k]=buf[j+k];
}
}
}
}
void INTT(int* a){
for(int i=t;i>=1;i--){
for(int j=0;j<(1<<t);j+=(1<<i)){
for(int k=0;k<(1<<i-1);k++){
buf[j+k]=((long long)(a[j+k]+a[j+k+(1<<i-1)])+mod)%mod*inv_2%mod;
buf[j+k+(1<<i-1)]=((long long)(a[j+k]-a[j+k+(1<<i-1)])+mod)%mod*inv_2%mod*inv_omega[k<<L-i]%mod;
}
for(int k=0;k<(1<<i);k++){
a[j+k]=buf[j+k];
}
}
}
for(int i=1;i<=t;i++){
for(int j=0;j<(1<<t);j+=(1<<i)){
for(int k=0;k<(1<<i-1);k++){
buf[j+k*2]=a[j+k];
buf[j+k*2+1]=a[j+k+(1<<i-1)];
}
for(int k=0;k<(1<<i);k++){
a[j+k]=buf[j+k];
}
}
}
}
long long pow_mod(long long num1,long long num2){
num1=(num1%mod+mod)%mod;
num2=(num2%(mod-1)+(mod-1))%(mod-1);
long long num3=1;
while(num2){
if(num2&1){
num3=num3*num1%mod;
}
num1=num1*num1%mod;
num2>>=1;
}
return num3;
}
long long fact[50010],invfact[50010];
long long p[60][50010],h[60][50010],f[60];
int time_fact[N],fact_ans[N],fact_final[N];
int main(){
init();
int n,k;
scanf("%d %d",&n,&k);
fact[0]=invfact[0]=1;
for(int i=1;i<=n*k;i++){
fact[i]=fact[i-1]*i%mod;
invfact[i]=invfact[i-1]*pow_mod(i,-1)%mod;
}
fact_ans[0]=1;
p[0][0]=1;
h[0][0]=1;
for(int i=0;i<=k-1;i++){
time_fact[i]=invfact[i];
}
t=16;
NTT(time_fact);
long long ans=0;
for(int i=1;i<=n;i++){
NTT(fact_ans);
for(int j=0;j<N;j++){
fact_final[j]=(long long)time_fact[j]*fact_ans[j]%mod;
fact_ans[j]=fact_final[j];
}
INTT(fact_ans);
INTT(fact_final);
for(int j=0;j<=i*(k-1);j++){
p[i][j]=(long long)fact_final[j]*fact[j]%mod;
h[i][j]=p[i][j]*pow_mod(i,-j)%mod;
f[i]=(f[i]+h[i][j])%mod;
}
f[i]*=n*pow_mod(i,-1)%mod;
f[i]%=mod;
long long pre=fact[n]*invfact[i]%mod*invfact[n-i]%mod*f[i]%mod;
if(i%2==1){
ans+=pre;
ans%=mod;
}
else{
ans-=pre;
ans=(ans%mod+mod)%mod;
}
}
printf("%lld",ans);
return 0;
}
然而上面这份代码只能在 qoj 通过，uoj 过不去。
我们考虑三次 NTT 还是太多余了，很容易优化到一次：
#include<iostream>
#include<cstdio>
using namespace std;
const int mod=998244353,root=15311432,inv_2=499122177;
const int L=23,N=1<<16;
int omega[1<<L],inv_omega[1<<L];
void init(){
omega[0]=1;
for(int i=1;i<(1<<L);i++){
omega[i]=(long long)omega[i-1]*root%mod;
}
inv_omega[0]=1;
for(int i=1;i<(1<<L);i++){
inv_omega[i]=omega[(1<<L)-i];
}
}
int t,buf[N];
void NTT(int* a){
for(int i=t;i>=1;i--){
for(int j=0;j<(1<<t);j+=(1<<i)){
for(int k=0;k<(1<<i-1);k++){
buf[j+k]=a[j+k*2];
buf[j+k+(1<<i-1)]=a[j+k*2+1];
}
for(int k=0;k<(1<<i);k++){
a[j+k]=buf[j+k];
}
}
}
for(int i=1;i<=t;i++){
for(int j=0;j<(1<<t);j+=(1<<i)){
for(int k=0;k<(1<<i-1);k++){
buf[j+k]=(a[j+k]+(long long)omega[k<<L-i]*a[j+k+(1<<i-1)]%mod+mod)%mod;
buf[j+k+(1<<i-1)]=(a[j+k]-(long long)omega[k<<L-i]*a[j+k+(1<<i-1)]%mod+mod)%mod;
}
for(int k=0;k<(1<<i);k++){
a[j+k]=buf[j+k];
}
}
}
}
void INTT(int* a){
for(int i=t;i>=1;i--){
for(int j=0;j<(1<<t);j+=(1<<i)){
for(int k=0;k<(1<<i-1);k++){
buf[j+k]=((long long)(a[j+k]+a[j+k+(1<<i-1)])+mod)%mod*inv_2%mod;
buf[j+k+(1<<i-1)]=((long long)(a[j+k]-a[j+k+(1<<i-1)])+mod)%mod*inv_2%mod*inv_omega[k<<L-i]%mod;
}
for(int k=0;k<(1<<i);k++){
a[j+k]=buf[j+k];
}
}
}
for(int i=1;i<=t;i++){
for(int j=0;j<(1<<t);j+=(1<<i)){
for(int k=0;k<(1<<i-1);k++){
buf[j+k*2]=a[j+k];
buf[j+k*2+1]=a[j+k+(1<<i-1)];
}
for(int k=0;k<(1<<i);k++){
a[j+k]=buf[j+k];
}
}
}
}
int pow_mod(int num1,int num2){
num1=(num1%mod+mod)%mod;
num2=(num2%(mod-1)+(mod-1))%(mod-1);
int num3=1;
while(num2){
if(num2&1){
num3=(long long)num3*num1%mod;
}
num1=(long long)num1*num1%mod;
num2>>=1;
}
return num3;
}
int fact[50010],invfact[50010];
int p[60][50010],h[60][50010],f[60];
int time_fact[N],fact_ans[N],fact_final[N],tmp_final[N];
int main(){
init();
int n,k;
scanf("%d %d",&n,&k);
fact[0]=invfact[0]=1;
for(int i=1;i<=n*k;i++){
fact[i]=(long long)fact[i-1]*i%mod;
invfact[i]=(long long)invfact[i-1]*pow_mod(i,-1)%mod;
}
fact_final[0]=1;
p[0][0]=1;
h[0][0]=1;
for(int i=0;i<=k-1;i++){
time_fact[i]=invfact[i];
}
while((1<<t)<=n*(k-1)){
t++;
}
NTT(time_fact);
NTT(fact_final);
int ans=0;
for(int i=1;i<=n;i++){
for(int j=0;j<N;j++){
tmp_final[j]=fact_final[j]=(long long)time_fact[j]*fact_final[j]%mod;
}
INTT(tmp_final);
for(int j=0;j<=i*(k-1);j++){
p[i][j]=(long long)tmp_final[j]*fact[j]%mod;
h[i][j]=(long long)p[i][j]*pow_mod(i,-j)%mod;
f[i]=(f[i]+h[i][j])%mod;
}
f[i]=(long long)n*pow_mod(i,-1)%mod*f[i]%mod;
int pre=(long long)fact[n]*invfact[i]%mod*invfact[n-i]%mod*f[i]%mod;
if(i%2==1){
ans+=pre;
ans%=mod;
}
else{
ans-=pre;
ans=(ans%mod+mod)%mod;
}
}
printf("%d",ans);
return 0;
}
浙公网安备 33010602011771号