ARC136F - PGF 学习笔记初步

Description

给定一个 \(h\times w\) 的 \(01\) 矩阵和非负整数序列 \(\{b_n\}\)，接下来每秒会在 \(h\times w\) 个格子中均匀随机地选取一个将其取反。

问期望第几秒会第一次满足：对于任意的 \(i\) 有第 \(i\) 行恰有 \(b_i\) 个 \(1\)。

\(1\le h,w\le 50\)。

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 - 解决「第一次」问题

我们注意到本题中想要的是「第一次符合条件的时间」，这是很烦的，我们考虑这样子处理：

我们记 \(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_1}^{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=1}^n 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)^{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}})^{2(j+k)-n}(-1)^{n-i-k}\\ & = \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)^{n-i-(k-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)^{n-i-(k-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}\) 泰勒展开得到：

\[1+x^{\frac{2k-n}{n}}+(x^{\frac{2k-n}{n}})^2+(x^{\frac{2k-n}{n}})^3+\cdots \]

然后再化为 OGF 的封闭形式：

\[\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{2t-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}\)。

注意别忘记除掉 \(\frac{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

#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;
}