ARC136F - PGF 学习笔记初步
lsj2009
·
·
题解
Description
给定一个 h\times w 的 01 矩阵和非负整数序列 \{b_n\}，接下来每秒会在 h\times w 个格子中均匀随机地选取一个将其取反。
问期望第几秒会第一次满足：对于任意的 i 有第 i 行恰有 b_i 个 1。
## Solution
闲话：
我终于会 PGF 了！！！！
这真的只是铜牌题吗？？？我只能说，too hard for me。
由于笔者是 PGF 初学者，所以本篇文章**会讲得非常详细**以来加深笔者自己地理解和帮助其他初学者。
下面进入正文：
### 前置知识：
1. 基本组合数学知识（如二项式定理等）。
2. 基本微积分知识（如求导，泰勒展开等）。
阅读此文你**并不需要有关于概率生成函数（PGF）的知识**，因为本文可能会从头到尾地介绍一遍。
你甚至可能**不需要普通生成函数（OGF）和指数生成函数（EGF）的知识**（虽然理论上是这样，但是我觉得想要做这题还是得了解生成函数）。
### 符号约定：
- $F(x)$ 表示序列 $f$ 的普通生成函数（OGF）：$F(x)=\sum\limits_{i\ge 0} f_ix^i$。
- $\hat{F}(x)$ 表示序列 $f$ 的指数生成函数（EGF）：$\hat{F}(x)=\sum\limits_{i\ge 0} \frac{f_i}{i!}x^i=\sum\limits_{i\ge 0} \frac{[x^i]F(x)}{i!}x^i$。
- $[x^i]F(x)$ 表示多项式 $F(x)$ 其第 $i$ 次项的系数。
- $F*G$ 表示的是两个多项式卷积得到的多项式，具体的，若 $H=F*G$，则我们有：
- $[x^n]H(x)=\sum\limits_{m=0}^n [x^m]F(x)\times [x^{n-m}]G(x)$。
- 这个东西也可以写作 $H(x)=F(x)G(x)$。
- $F'(x)$ 表示的是对多项式 $F(x)$ 求导后得到的多项式。
- $\Pr(A)$ 表示的是 $A$ 事件成立的概率。
- $\operatorname{E}(X)$ 表示的是随机变量 $X$ 的期望取值。
### Pre - 概率生成函数（PGF）
如果我们现在有一个**取值是非负整数**的随机变量 $X$，则我们记其概率生成函数（PGF）为：
$$
F_X(x)=\sum\limits_{i\ge 0} \Pr(X=i)x^i
$$
下面简记为 $F(x)$。
或者从另一个角度理解：即序列 $\{\Pr(X=1),\Pr(X=2),\cdots\}$ 的普通生成函数（OGF）。
更一般的，我们的 PGF 可能是关于**一系列事件**而非一个随机变量的（特别的，我们可以将 $X=1,2,\cdots$ 看作是**一系列不同的事件**）。
在这题中，我们的随机变量 $X$ 也就是**停时时间**（即第一次符合条件的时间）。
则我们的答案（**期望停时时间**）就是：
$$
\operatorname{E}(X)=\sum\limits_{i\ge 0} \Pr(X=i)\cdot i
$$
我们发现：如果对 $F(x)$ 求导，我们就可以得到：
$$
F'(x)=\sum\limits_{i\ge 0} \Pr(X=i)\cdot i\cdot x^{i-1}
$$
如果我们此时带入 $x=1$，就可以惊讶地发现：
$$
\begin{aligned}
F'(1) & = \sum\limits_{i\ge 0} \Pr(X=i)\cdot i\cdot 1^{i-1}\\
& =\sum\limits_{i\ge 0} \Pr(X=i)\cdot i\\
& =\operatorname{E}(X)
\end{aligned}
$$
也就是说，$F'(1)$ 就是我们想要的答案！
### Pre - 解决「第一次」问题
我们令 $F(x)$ 为答案的 PGF。
我们注意到本题中想要的是「**第一次**符合条件的时间」，这是很烦的，我们考虑这样子处理：
我们记 $S,T$ 分别为题目中的初始和最终局面，事件 $A(U\to V,i)$ 表示的是从初始是 $U$ 局面经过**恰好** $i$ 时刻变成了 $V$ 局面，则我们列出关于 $A(S\to T,*)$ 的 PGF $G(x)$ 和关于 $A(T\to T,*)$ 的 PGF $H(x)$：
$$
G(x)=\sum\limits_{i\ge 0} \Pr(A(S\to T,i))x^i\\
H(x)=\sum\limits_{i\ge 0} \Pr(A(T\to T,i))x^i
$$
则我们有：$G=F*H$。
- 证明：我们可以认为是，当我们在 $i$ 时刻**第一次**从初始状态到达了最终状态，然后**又经过**了 $j$ 时刻的一系列操作**从最终状态又回到了最终状态**；又因为这两个事件是独立的，则我们有 $[x^i]F(x)\times [x^j]H(x)\to [x^{i+j}]G(x)$，即 $G=F*H$。
注意到我们关心的是 $F'(1)$ 的值，则根据 $G=F*H$，我们有 $F=\frac{G}{H}$，根据商的求导法则，我们有：
$$
F'(1)=\frac{G'(1)H(1)-H'(1)G(1)}{H^2(1)}
$$
则我们现在将问题转换为了求 $G(1),G'(1),H(1),H'(1)$ 的值，而这样子就消除了「第一次」的限制。
下面我们默认只讨论 $G(1),G'(1)$ 的求解，对于 $H(1),H'(1)$ 则是几乎一模一样的，在实现上只需要改一下进行 dp 时传进去的参数即可。
### Part 1 - 一个 dp
注意到我们一个格子**翻了奇数次相当于翻了一次，翻了偶数次则相当于没翻**。
我们不妨假定：每个格子**最多只会被翻转一次**。考虑计算从 $S\to T$ 的方案总数。
注意到我们只关心每一行 $1$ 的个数，所以考虑按行 dp。
记录 $f_{i,j}$ 表示考虑前 $i$ 行翻转恰好 $j$ 次使得满足条件的的方案总数。
我们记第 $i$ 行在初始状态有 $a_i$ 个 $1$，而最终状态就是题目给定的有 $b_i$ 个 $1$。
则考虑枚举将多少个 $1$ change 成了 $0$。
如果有 $k$ 个 $0\to 1$，则我们应有 $b_i-(a_i-k)$ 个 $1\to 0$，则得到转移方程：
$$
f_{i-1,j}\binom{a_i}{k}\binom{w-a_i}{b_i-(a_i-k)}\to f_{i,j+k+(b_i-(a_i-k))}
$$
记 $n=h\times w$ 为格子总数量，则这个 dp 是 $\mathcal{O}(n^2)$ 的，可以接受。
我们设 $\{c_n\}$ 为最终的答案，即 $c_i=f_{h,i}$。
### Part 2 - 构造生成函数
现在我们对于每一个格子给出其被翻转偶数 or 奇数次的 PGF：
$$
P_0(x)=\sum\limits_{i\equiv0\pmod{2}} \frac{x^i}{n^i}\\
P_1(x)=\sum\limits_{i\equiv1\pmod{2}} \frac{x^i}{n^i}
$$
这时候我们看起来只需要枚举实际上有多少个格子被翻转（即被翻转奇数次），然后将对应的生成函数卷起来即可！
但是实际上并不是，因为**不同格子之间的翻转是区分顺序的**，所以我们要使用 EGF 进行区分（这是因为两个 EGF 的卷积表述成最终的序列是会乘上一个组合数的系数），即得到：
$$
\hat{G}(x)=\sum\limits_{i\ge 0} c_i\hat{P_1}^i(x)\hat{P_0}^{n-i}(x)
$$
尝试推导 $\hat{P_0}$ 和 $\hat{P_1}$。
我们给出：
$$
P(x)=\sum\limits_{i\ge 0} \frac{x^i}{n^i}\\
Q(x)=\sum\limits_{i\ge 0} (-1)^i\frac{x^i}{n^i}=\sum\limits_{i\ge 0} \frac{x^i}{(-n)^i}
$$
可以得到：
$$
\begin{aligned}
\hat{P}(x) & = \sum\limits_{i\ge 0} \frac{x^i}{i!n^i}\\
& = \sum\limits_{i\ge 0} \frac{(\frac{x}{n})^i}{i!}\\
& = e^{\frac{x}{n}}
\end{aligned}
$$
最后一步是考虑 $e^x$ 的泰勒展开。
同理我们可以得到的是：$\hat{Q}(x)=e^{-\frac{x}{n}}$。
我们注意到有：
$$
P_0(x)=\frac{P(x)+Q(x)}{2}
$$
则可以得到的是：
$$
\hat{P_0}(x)=\frac{\hat{P}(x)+\hat{Q}(x)}{2}
$$
也就是对等式两边同时取 EGF，因为我们注意到 $[x^n]$ 的相对关系没有发生变化，所以等号仍然成立。
即
$$
\hat{P_0}(x)=\frac{e^{\frac{x}{n}}+e^{-\frac{x}{n}}}{2}
$$
同理有：
$$
\hat{P_1}(x)=\frac{e^{\frac{x}{n}}-e^{-\frac{x}{n}}}{2}
$$
则我们得到了：
$$
\hat{G}(x)=\sum\limits_{i\ge 0} c_i\left(\frac{e^{\frac{x}{n}}-e^{-\frac{x}{n}}}{2}\right)^i\left(\frac{e^{\frac{x}{n}}+e^{-\frac{x}{n}}}{2}\right)^{n-i}
$$
### Part 2 - 一个柿子的推
考虑化简上面 EGF：
$$
\begin{aligned}
\hat{G}(x) & = \sum\limits_{i\ge 0} c_i\left(\frac{e^{\frac{x}{n}}-e^{-\frac{x}{n}}}{2}\right)^i\left(\frac{e^{\frac{x}{n}}+e^{-\frac{x}{n}}}{2}\right)^{n-i}\\
& = \frac{1}{2^n}\sum\limits_{i=1}^n c_i\left(e^{\frac{x}{n}}-e^{-\frac{x}{n}}\right)^i\left(e^{\frac{x}{n}}+e^{-\frac{x}{n}}\right)^{n-i}\\
& = \frac{1}{2^n}\sum\limits_{i=1}^n c_i\sum\limits_{j=0}^i \binom{i}{j}(e^{\frac{x}{n}})^j(-e^{-\frac{x}{n}})^{i-j}\sum\limits_{k=0}^{n-i} \binom{n-i}{k}(e^{\frac{x}{n}})^k(e^{-\frac{x}{n}})^{n-i-k}\\
& = \frac{1}{2^n}\sum\limits_{i=1}^n c_i\sum\limits_{j=0}^i \sum\limits_{k=0}^{n-i}\binom{i}{j}\binom{n-i}{k}(e^{\frac{x}{n}})^{j+k}(e^{-\frac{x}{n}})^{n-(j+k)}(-1)^{i-j}\\
& = \frac{1}{2^n}\sum\limits_{i=1}^n c_i\sum\limits_{j=0}^i \sum\limits_{k=0}^{n-i}\binom{i}{j}\binom{n-i}{k}(e^{\frac{x}{n}})^{2(j+k)-n}(-1)^{i-j}\\
& = \frac{1}{2^n}\sum\limits_{k=0}^{n}e^{(\frac{2k-n}{n})x}\sum\limits_{i=1}^n c_i\sum\limits_{j=0}^i \binom{i}{j}\binom{n-i}{k-j}(-1)^{i-j}\\
& = \frac{1}{2^n}\sum\limits_{k=0}^{n}e^{(\frac{2k-n}{n})x}\sum\limits_{i=1}^n c_i[x^k](x-1)^{i}(x+1)^{n-i}\\
& = \frac{1}{2^n}\sum\limits_{k=0}^{n}e^{(\frac{2k-n}{n})x}[x^k]\sum\limits_{i=1}^n c_i(x-1)^{i}(x+1)^{n-i}\\
\end{aligned}
$$
其中第三个等号处用二项式定理将其暴力展开，倒数第三个等号处改为枚举 $j+k$，倒数第二个等号处将 $\sum\limits_{j=0}^i \binom{i}{j}\binom{n-i}{k-j}(-1)^{i-j}$ 改为其 OGF 写法（考虑二项式定理将其展开即可获证）。
### Part 3 - 计算答案
我们先暴力计算出 $\sum\limits_{i=1}^n c_i(x-1)^{i}(x+1)^{n-i}$ 的多项式 $C(x)$。
则柿子化为
$$
\hat{G}(x)=\frac{1}{2^n}\sum\limits_{k=0}^{n}e^{(\frac{2k-n}{n})x}[x^k]C(x)
$$
这时候我们去掉 EGF，还原回 OGF：
$$
G(x)=\frac{1}{2^n}\sum\limits_{k=0}^{n}\frac{[x^k]C(x)}{1-(\frac{2k-n}{n})x}
$$
这一部分即为先将 $e^{(\frac{2k-n}{n})x}$ 泰勒展开得到：
$$
\frac{1}{0!}+\frac{x^{\frac{2k-n}{n}}}{1!}+\frac{(x^{\frac{2k-n}{n}})^2}{2!}+\frac{(x^{\frac{2k-n}{n}})^3}{3!}+\cdots
$$
由于 EGF->OGF 会将 $[x^i]$ 去 $\times i!$，所以其 OGF 为：
$$
1+x^{\frac{2k-n}{n}}+(x^{\frac{2k-n}{n}})^2+(x^{\frac{2k-n}{n}})^3+\cdots
$$
然后再化为封闭形式：
$$
\frac{1}{1-(\frac{2k-n}{n})x}
$$
此时我们代入 $x=1$ 即可得到 $G(1)$ 的值，然而有一个问题是当 $k=n$ 时我们会得到 $\frac{1}{0}$，这无法处理。
我们考虑令 $A(x)=(1-x)G(x),B(x)=(1-x)H(x)$，此时我们仍有 $F=\frac{A}{B}$，所以不影响答案。
而这时当 $k=n$ 时我们外层的 $1-x$ 和内层的 $1-x$ 相互抵消，所以对 $A(1)$ 恰好贡献为 $1\times [x^n]C(x)$。
而对于 $k\ne n$ 时，我们记 $t=\frac{2k-n}{n},v=[x^k]C(x)$，则所得多项式为 $\frac{v(x-1)}{1-tx}$。
- 这时 $x-1=0$，不对 $A(1)$ 产生贡献。
- 此时导函数为 $\frac{(t-1)v}{(1-tx)^2}$，对 $A'(1)$ 贡献为 $\frac{v}{t-1}$。
注意别忘记除掉 $2^n$。
同理我们可以计算出 $B(1)$ 和 $B'(1)$ 的值，则可以得到：
$$
F'(1)=\frac{A'(1)B(1)-B'(1)A(1)}{B^2(1)}
$$
这就是答案。
暴力实现复杂度是 $\mathcal{O}(n^2)$ 的，已经可以通过；使用一些多项式知识可以做到 $\mathcal{O}(n\operatorname{polylog}(n))$。
下面给出 $\mathcal{O}(n^2)$ 的实现。
## Code
```cpp
#include<bits/stdc++.h>
//#pragma GCC optimize(3,"Ofast","inline")
//#define int long long
#define i128 __int128
#define ll long long
#define ull unsigned long long
#define uint unsigned int
#define ld double
#define PII pair<int,int>
#define INF 0x3f3f3f3f
#define INFLL 0x3f3f3f3f3f3f3f3f
#define chkmax(a,b) a=max(a,b)
#define chkmin(a,b) a=min(a,b)
#define rep(k,l,r) for(int k=l;k<=r;++k)
#define per(k,r,l) for(int k=r;k>=l;--k)
#define cl(f,x) memset(f,x,sizeof(f))
#define pcnt(x) __builtin_popcount(x)
#define lg(x) (31-__builtin_clz(x))
using namespace std;
void file_IO() {
// system("fc .out .ans");
freopen(".in","r",stdin);
freopen(".out","w",stdout);
}
bool M1;
template<int p>
struct mint {
int x;
mint() {
x=0;
}
mint(int _x) {
x=_x;
}
friend mint operator + (mint a,mint b) {
return a.x+b.x>=p? a.x+b.x-p:a.x+b.x;
}
friend mint operator - (mint a,mint b) {
return a.x<b.x? a.x-b.x+p:a.x-b.x;
}
friend mint operator * (mint a,mint b) {
return 1ll*a.x*b.x%p;
}
friend mint operator ^ (mint a,ll b) {
mint res=1,base=a;
while(b) {
if(b&1)
res*=base;
base*=base; b>>=1;
}
return res;
}
friend mint operator ~ (mint a) {
return a^(p-2);
}
friend mint operator / (mint a,mint b) {
return a*(~b);
}
friend mint & operator += (mint& a,mint b) {
return a=a+b;
}
friend mint & operator -= (mint& a,mint b) {
return a=a-b;
}
friend mint & operator *= (mint& a,mint b) {
return a=a*b;
}
friend mint & operator /= (mint& a,mint b) {
return a=a/b;
}
friend mint operator ++ (mint& a) {
return a+=1;
}
friend mint operator -- (mint& a) {
return a-=1;
}
};
const int MOD=998244353;
#define mint mint<MOD>
const int N=5e3+5;
mint jc[N],inv_jc[N];
void init(int n=5000) {
jc[0]=1;
rep(i,1,n)
jc[i]=jc[i-1]*i;
inv_jc[n]=~jc[n];
per(i,n-1,0)
inv_jc[i]=inv_jc[i+1]*(i+1);
}
mint C(int n,int m) {
if(m<0||n<m)
return 0;
return jc[n]*inv_jc[n-m]*inv_jc[m];
}
int n,h,w;
mint tmp[N][N],t[N];
void calc(int a[],int b[],mint c[]) {
rep(i,0,h) {
rep(j,0,n)
tmp[i][j]=0;
}
tmp[0][0]=1;
rep(i,1,h) {
rep(j,0,w)
t[j]=0;
rep(j,0,a[i]) {
int d=b[i]-a[i]+j;
if(j+d>=0)
t[j+d]+=C(a[i],j)*C(w-a[i],d);
}
rep(j,0,w) {
rep(k,0,n-j)
tmp[i][j+k]+=tmp[i-1][k]*t[j];
}
}
rep(i,0,n)
c[i]=tmp[h][i];
}
mint p[N];
void calc(mint c[],mint f[]) {
rep(i,0,n)
p[i]=C(n,i);
rep(i,0,n) {
if(i) {
rep(j,1,n)
p[j]-=p[j-1];
per(j,n,1)
p[j]=p[j-1]-p[j];
p[0]=MOD-p[0];
}
rep(j,0,n)
f[j]+=p[j]*c[i];
}
}
int a[N],b[N];
char s[N];
mint c1[N],c2[N],f[N],g[N];
void solve() {
scanf("%d%d",&h,&w);
rep(i,1,h) {
scanf("%s",s+1);
int cnt=0;
rep(j,1,w)
cnt+=s[j]=='1';
a[i]=cnt;
}
rep(i,1,h)
scanf("%d",&b[i]);
n=h*w;
calc(a,b,c1);
calc(c1,f);
calc(b,b,c2);
calc(c2,g);
mint sf=0,sg=0,sdf=0,sdg=0;
rep(i,0,n-1) {
mint val=(mint(2*i)-mint(n))/mint(n);
sdf+=f[i]/(val-mint(1));
sdg+=g[i]/(val-mint(1));
}
sf+=f[n];
sg+=g[n];
mint inv=(~mint(2))^n;
sf*=inv; sg*=inv;
sdf*=inv; sdg*=inv;
printf("%d\n",((sdf*sg-sf*sdg)/(sg*sg)).x);
}
bool M2;
signed main() {
//file_IO();
int testcase=1;
init();
//scanf("%d",&testcase);
while(testcase--)
solve();
fprintf(stderr,"used time = %ldms\n",1000*clock()/CLOCKS_PER_SEC);
fprintf(stderr,"used memory = %lldMB\n",(&M2-&M1)/1024/1024);
return 0;
}
```