「Luogu P4707」重返现世

Problem

Description

n 种元素,每个单位时间会随机取得一种元素,取得第 i 种元素的概率为 \frac{p_{i}}{m}
求取得 k 种不同元素的期望时间。

答案对 998244353 取模。

Constraints

1\le n\le 10^{3},\max\left(1,n-10\right)\le k\le n,1\le m\le 10^{4},\sum_{i=1}^{n}p_{i}=m

Solution

Solution 1

orz Sdchr.

Analysis

应用扩展 \min - \max 容斥可得答案为:

\begin{aligned} &\sum_{T\subseteq\left\{1,2,...,n\right\}}(-1)^{|T|-n+k-1}\binom{|T|-1}{n-k}\frac{m}{p_{T}}\\ =&m\sum_{i=n-k+1}^{n}(-1)^{i-n+k-1}\binom{i-1}{n-k}\sum_{s=1}^{m}\frac{f_{n,i,s}}{s}\\ =&m\sum_{s=1}^{m}\frac{1}{s}\sum_{i=n-k+1}^{n}(-1)^{i-n+k-1}\binom{i-1}{n-k}f_{n,i,s} \end{aligned}

其中:

p_{T}=\sum_{i\in T}p_{i}

f_{i,j,s}=\sum_{T\in\left\{1,2,...,i\right\},|T|=j}[p_{T}=s]=f_{i-1,j,s}+f_{i-1,j-1,s-p_{i}}

f 进行 DP 显然是 O(nmk) 的,于是考虑将 \frac{1}{s} 的系数作为一个整体进行 DP:

g_{i,j,s}=\sum_{sz=j}^{i}(-1)^{sz-j}\binom{sz-1}{j-1}f_{i,sz,s}

即在仅考虑前 i 个元素的前提下概率和为 s 的集合的带系数贡献之和。

发现当 sz<j 时, \displaystyle\binom{sz-1}{j-1}=0 ,当 i<sz 时, f_{i,sz,s}=0 ,故上式亦可写成:

g_{i,j,s}=\sum_{sz=1}^{n}(-1)^{sz-j}\binom{sz-1}{j-1}f_{i,sz,s}

考虑 g 的转移方程,有:

g_{i,j,s}=g_{i-1,j,s}+g_{i-1,j-1,s-p_{i}}+x

\begin{aligned} x &= g_{i,j,s}-g_{i-1,j,s}-g_{i-1,j-1,s-p_{i}} \\ &= \sum_{sz=1}^{n}(-1)^{sz-j}\binom{sz-1}{j-1}\left(f_{i,sz,s}-f_{i-1,sz,s}\right)-\sum_{sz=1}^{n}(-1)^{sz-j+1}\binom{sz-1}{j-2}f_{i-1,sz,s-p_{i}} \\ &= \sum_{sz=1}^{n}(-1)^{sz-j}\binom{sz-1}{j-1}f_{i-1,sz-1,s-p_{i}}-\sum_{sz=1}^{n}(-1)^{sz-j+1}\binom{sz-1}{j-2}f_{i-1,sz,s-p_{i}} \\ &= \sum_{sz=1}^{n}(-1)^{sz-j+1}\binom{sz}{j-1}f_{i-1,sz,s-p_{i}}-\sum_{sz=1}^{n}(-1)^{sz-j+1}\binom{sz-1}{j-2}f_{i-1,sz,s-p_{i}} \\ &= \sum_{sz=1}^{n}(-1)^{sz-j+1}\left(\binom{sz}{j-1}-\binom{sz-1}{j-2}\right)f_{i-1,sz,s-p_{i}} \\ &= \sum_{sz=1}^{n}(-1)^{sz-j+1}\binom{sz-1}{j-1}f_{i-1,sz,s-p_{i}} \\ &= -g_{i-1,j,s-p_{i}} \end{aligned}

故:

g_{i,j,s}=g_{i-1,j,s}+g_{i-1,j-1,s-p_{i}}-g_{i-1,j,s-p_{i}}

边界条件为 g_{i,0,0}=1

时间复杂度 O(nm(n - k))

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
#include<cstdio>
#include<cstring>
#include<algorithm>

using i64=long long;
using uchar=unsigned char;

constexpr int maxk(10);
constexpr int maxm(10000);
constexpr int p(998244353);

template<typename _Tp>
inline void swap(_Tp&x,_Tp&y)
{_Tp z=x;x=y;y=z;}

namespace IOManager{
constexpr int FILESZ(131072);
char buf[FILESZ];
const char*ibuf=buf,*tbuf=buf;

struct IOManager{
inline char gc()
{return(ibuf==tbuf)&&(tbuf=(ibuf=buf)+fread(buf,1,FILESZ,stdin),ibuf==tbuf)?EOF:*ibuf++;}

template<typename _Tp>
inline operator _Tp(){
_Tp s=0u;char c=gc();
for(;c<48;c=gc());
for(;c>47;c=gc())
s=(_Tp)(s*10u+c-48u);
return s;
}
};
}IOManager::IOManager io;

int f[2][maxk+2][maxm+1],
inv[maxm+1];

int main(){
const auto pls=[](const auto&x,const auto&y){return x+y<p?x+y:(x+y-p);};
const auto sub=[](const auto&x,const auto&y){return x<y?x-y+p:(x-y);};

const int n=io,k=n-(int)io+1,m=io;

inv[1]=1;
for(int i=2;i<=m;++i)
inv[i]=(i64)inv[p%i]*(p-p/i)%p;

auto cur=f[0],las=f[1];
las[0][0]=1;
for(int i=1,sp=0;i<=n;++i,swap(cur,las)){
memset(cur,0,(sizeof f)>>1);
cur[0][0]=1;

const int cp=io;sp+=cp;
for(int j=1;j<=k;++j){
std::copy(las[j],las[j]+cp,cur[j]);
for(int s=cp;s<=sp;++s)
cur[j][s]=sub(pls(las[j][s],las[j-1][s-cp]),las[j][s-cp]);
}
}

i64 ans=0;
for(int s=1;s<=m;++s)
ans+=(i64)las[k][s]*inv[s]%p;
printf("%d",ans%p*m%p);

return 0;
}

Solution 2

orz alan_cty.

Analysis

容易得到第 i 种元素的 EGF 为 \hat{f}_{i}(z)=e^{\frac{p_{i}}{m}z}-1
枚举最后取得的元素 v ,设 g_{v}(z) 为取得 v 前序列的 EGF,易有:

g_{v}(z)=\sum_{T\in\left\{1,2,...,n\right\}-\left\{v\right\},|T|=k-1}\prod_{i\in T}\hat{f}_{i}(z)

则答案即为:

\sum_{v=1}^{n}\frac{p_{v}}{m}\sum_{0\le i}\left(i+1\right)i!\left[z^{i}\right]g_{v}(z)

考虑将 g_{v}(z) 写为 \sum_{0\le s}a_{v,s}e^{\frac{s}{m}z} 的形式,易知 e^{\frac{s}{m}z} 对答案的贡献为:

\sum_{0\le i}\left(i+1\right)i!\left[z^{i}\right]e^{\frac{s}{m}z}=\sum_{0\le i}\left(i+1\right)\left(\frac{s}{m}\right)^{i}=\left(\frac{m}{m-s}\right)^{2}

对于 a_{v,s} ,可先 DP 求出 \displaystyle{\prod_{i=1}^{n}\hat{f}_{i}(z)} 中各项的系数,在枚举 v 时将 v 的贡献减去即可。

则答案为:

\sum_{v=1}^{n}\frac{p_{v}}{m}\sum_{s=1}^{m}a_{v,s}\left(\frac{m}{m-s}\right)^{2}=m\sum_{v=1}^{n}p_{v}\sum_{s=1}^{m}\frac{a_{v,s}}{\left(m-s\right)^{2}}

时间复杂度 O\left(nm\left(n-k\right)\right)

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
#include<cstdio>
#include<cstring>
#include<algorithm>

using i64=long long;
using uchar=unsigned char;

constexpr int maxn(1000);
constexpr int maxk(10);
constexpr int maxm(10000);
constexpr int p(998244353);

namespace IOManager{
constexpr int FILESZ(131072);
char buf[FILESZ];
const char*ibuf=buf,*tbuf=buf;

struct IOManager{
inline char gc()
{return(ibuf==tbuf)&&(tbuf=(ibuf=buf)+fread(buf,1,FILESZ,stdin),ibuf==tbuf)?EOF:*ibuf++;}

template<typename _Tp>
inline operator _Tp(){
_Tp s=0u;char c=gc();
for(;c<48;c=gc());
for(;c>47;c=gc())
s=(_Tp)(s*10u+c-48u);
return s;
}
};
}IOManager::IOManager io;

int a[maxk+1][maxm+1],
pr[maxn+1],
inv[maxm+1];

int main(){
const auto op=[](const auto&v){return v?p-v:0;};
const auto pls=[](const auto&x,const auto&y){return x+y<p?x+y:(x+y-p);};
const auto sub=[](const auto&x,const auto&y){return x<y?x-y+p:(x-y);};

const int n=io,k=n-(int)io,m=io;

inv[1]=1;
for(int i=2;i<=m;++i)
inv[i]=(i64)inv[p%i]*(p-p/i)%p;
for(int i=2;i<=m;++i)
inv[i]=(i64)inv[i]*inv[i]%p;

a[0][0]=1;
for(int i=1;i<=n;++i){
pr[i]=io;
const int cp=pr[i];

for(int sz=k;sz;--sz){
int*const&cur=a[sz],*const&las=a[sz-1];
for(int s=m;cp<=s;--s)
cur[s]=sub(pls(cur[s-cp],las[s]),cur[s]);
for(int s=0;s!=cp;++s)
cur[s]=sub(las[s],cur[s]);
}
int*const&cur=a[0];
for(int s=m;cp<=s;--s)
cur[s]=sub(cur[s-cp],cur[s]);
for(int s=0;s!=cp;++s)
cur[s]=op(cur[s]);
}

int*const&use=a[k];
i64 ans=0;
for(int i=1;i<=n;++i){
const int cp=pr[i];

int*const&ocur=a[0];
for(int s=0;s!=cp;++s)
ocur[s]=op(ocur[s]);
for(int s=cp;s<=m;++s)
ocur[s]=sub(ocur[s-cp],ocur[s]);
for(int sz=1;sz<=k;++sz){
int*const&cur=a[sz],*const&las=a[sz-1];
for(int s=0;s!=cp;++s)
cur[s]=sub(las[s],cur[s]);
for(int s=cp;s<=m;++s)
cur[s]=sub(pls(cur[s-cp],las[s]),cur[s]);
}

i64 t=0;
for(int s=0;s<=m-cp;++s)
t+=(i64)use[s]*inv[m-s]%p;
ans+=t%p*cp%p;

for(int sz=k;sz;--sz){
int*const&cur=a[sz],*const&las=a[sz-1];
for(int s=m;cp<=s;--s)
cur[s]=sub(pls(cur[s-cp],las[s]),cur[s]);
for(int s=0;s!=cp;++s)
cur[s]=sub(las[s],cur[s]);
}
int*const&cur=a[0];
for(int s=m;cp<=s;--s)
cur[s]=sub(cur[s-cp],cur[s]);
for(int s=0;s!=cp;++s)
cur[s]=op(cur[s]);
}

printf("%d",ans%p*m%p);

return 0;
}