手动博客搬家: 本文发表于20181212 09:37:21, 原地址https://blog.csdn.net/suncongbo/article/details/84962727

呜啊怎么又是数学了啊。。。数学比例\(\frac{16}{33}=0.4848\)

orz yhx-12243神仙

题目链接: https://codeforces.com/contest/947/problem/E

题意:
有一个\([0,n]\)的随机数\(x\)初始为\(i\)的概率为\(p_i\). \(m\)次操作每次从\([0,x]\)中等概率随机选择一个数\(y\)把\(x\)变成\(y\). 对每个\(i=0,1,...,n\)求\(m\)次之后\(x=i\)的概率。

题解:
嗯，是个线代神题。本题解很长，希望大家都能耐心看懂。(当时我看题解也看了三小时多)
一、naïve的dp和矩乘优化
首先，一个\(O(nm)\)的\(dp\)很好想: 设\(f[i][j]\)表示\(i\)轮后\(x=j\)的概率，则轻易列出方程: \(f[i][j]=\sum^{n}_{k=i} \frac{dp[i-1][k]}{k+1}\).
按此\(dp\)时间复杂度\(O(nm)\), 太慢了。考虑优化。
写出转移矩阵: 令\((n+1)\)维列向量为\(dp[i]\)数组, 则
\(\bm f= \begin{bmatrix} 1 & \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & ... & \frac{1}{n+1}\\ 0 & \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & ... & \frac{1}{n+1}\\ 0 & 0 & \frac{1}{3} & \frac{1}{4} & ... & \frac{1}{n+1} \\ & & & & ... \\ 0 & 0 & 0 & 0 & ... & \frac{1}{n+1} \end{bmatrix}\times \bm f_{i-1}\)
令上面的转移矩阵为\(\bm A\). 如果直接利用矩阵乘法优化，时间复杂度\(O(n^3\log m)\), 还是太慢了。
即使是使用Hamilton-Cayley定理，复杂度也难以降到\(O(n^2)\)以下。
那么我们该怎么办呢？于是一个神奇的做法出现了——矩阵对角化。
二、对转移矩阵特征系统的深入了解
在对矩阵进行对角化之前，我们先看一个概念——特征系统（包括特征根和特征向量）。
在这里特征系统相关知识不再阐述，随便一本线代书上就会讲。
那我们观察转移矩阵\(\bm A\), 发现一个事实——\((n+1)\)阶矩阵\(\rm A\)共有\((n+1)\)个特征根，它们分别是\(1,\frac{1}{2},\frac{1}{3},...,\frac{1}{n+1}\). 这是因为\(\rm A\)是上对角矩阵，其特征多项式为\(f(\lambda)=|\lambda \bm I-\bm A|=\prod^{n+1}_{i=1} (\lambda-\frac{1}{i})\).
求出这个还不够，要对矩阵进行对角化，我们还需要求出它的特征向量。（什么你说为什么要求？一会就知道了T_T）
很好，我们对较小的矩阵进行手算，得出结论: 第\(i\)个特征根\(\lambda_i=\frac{1}{i+1}\)的特征向量为\[\bm v_i=\begin{bmatrix} (-1)^i & (-1)^{i+1}{i\choose 1} & (-1)^{i+2}{i\choose 2} & ... & (-1)^{i+n}{i\choose n}\end{bmatrix}^T\]即\(v_{i,j}=(-1)^{i+j}{i\choose j}\). (注意这里\(v_{i,j}\)表示的是第\(i\)个特征(列)向量\(\bm v_i\)的第\(j\)个元素，如果放到特征向量构成的矩阵中应该是第\(j\)行第\(i\)列. 为了避免混乱，下文将用\(v_{i,j}\)表示第\(i\)个特征列向量的第\(j\)个元素，\(V_{i,j}\)表示特征向量构成的矩阵中第\(i\)行第\(j\)列的元素，即\(V_{i,j}=v_{j,i}\).)
即: \[\bm V=\begin{bmatrix} 1 & -1 & 1 & -1 & ... & (-1)^n \\ 0 & 1 & -2 & 3 & ... & (-1)^{n+1}{n\choose 1} \\ 0 & 0 & 1 & -3 & ... & (-1)^{n+2}{n\choose 2} \\ & & & & ... \\ 0 & 0 & 0 & 0 & ... & (-1)^{n+n}{n\choose n} \end{bmatrix}\]
其中\(\bm V\)为每一个\(\bm v\)列向量一起构成的矩阵。一会我们将看到，这个矩阵将起到非常重要的作用。
嗯，我们刚刚通过对较小的矩阵进行手算的方式得到了特征向量\(\bm V\), 现在我们尝试通过理论推导去证明这个公式。
我们先形式化特征向量的表达式: \(v_{i,j}=(-1)^{i+j}{i\choose j}, V_{i,j}=(-1)^{i+j} {j\choose i}\)
根据特征向量的定义，对于特征(列)向量\(\bm v_i\)我们需要证明\(\bm A\bm v_i=\lambda_i\bm v_i\).
我们考虑写出式子: \(\bm A\bm v_i\)是一个\(n\times n\)方阵和\(n\)阶列向量的乘积，考虑乘积的第\(r\)个元素就是矩阵\(\bm A\)的第\(r\)行和\(\bm v_i\)的内积，则\[\sum^{n}_{j=0} A_{r,j}v_{i,j} =\sum^{n}_{j=r} \frac{1}{j+1}(-1)^{i+j}{i\choose j} \\ =\sum^{i}_{j=r} \frac{1}{j+1} (-1)^{i+j}\frac{i!}{j!(i-j)!}\\ =\frac{1}{i+1}\sum^{i}_{j=r}(-1)^{i+j}\frac{(i+1)!}{(j+1)!(i-j)!}\\ =\frac{1}{i+1}\sum^{i}_{j=r}(-1)^{i+j}{i+1\choose j+1}\\ =\frac{1}{i+1}\sum^{i}_{j=r}(-1)^{i+1}{j-i+1\choose j+1}\\ =\frac{1}{i+1}(-1)^{i+1}\sum^{i}_{j=r}{j-i-1\choose j+1}\\ = \frac{1}{i+1}(-1)^{i+1}({i-i-1+1\choose i+1}-{r-i-1\choose r})\\ =\frac{1}{i+1}(-1)^i{r-i-1\choose r}\\ =(-1)^{i+r}{i\choose r}=\lambda_i v_{i,r}\]
嗯，这一步推导值得仔细体会。一共包含九行，其中第5和9行的原理是\({-a\choose b}={a+b-1\choose b}(-1)^b\); 第7行的原理是对杨辉三角任何一斜列求和，等于最右下角元素的下一个减去最左上角元素的左一个：\(\sum^{n}_{j=0} {a+j\choose b+j}={a+n+1\choose b}-{a+n\choose b-1}\).
下面仔细审视如此奇怪的推导的意图: 我们发现第三行把\(\frac{1}{j+1}\)巧妙地转化成了\(\frac{1}{i+1}\), 第五行把\((-1)^{i+j}\)变成了\((-1)^{i+1}\), 然后这样一个组合数求和的如此难以处理的问题被去掉了杂质，变成了一个纯粹的杨辉三角一斜列的求和。这就是要百费周折地各种化成负的再化回来的原因。这一段推导真的有很多精妙的地方，希望自己能够借鉴。好了，恢复正题。
刚才我们成功证明了特征向量\(\bm v_i\)的表达式，然后这有什么用吗？

三、矩阵对角化加速运算——从\(O(n^3\log m)\)到\(O(n^2)\)
这里是本题的重点。
首先，我们考虑现在要做的事情：快速计算该矩阵的\(m\)次幂。假设这个矩阵是一个对角矩阵\(\rm diag(x_0,x_1,...,x_n)\)，则它的\(m\)次方非常容易计算: \(\rm diag(x_0^m,x_1^m,...,x_n^m)\).
我们现在希望快速计算矩阵\(\bm f\)的\(m\)次幂，于是我们可以考虑把它转化为对角矩阵的幂运算。
我们找到一个可逆矩阵\(\bm \Phi\)满足\(\bm\Phi^{-1}\bm A\bm\Phi=\bm B\)其中\(\bm B\)为对角矩阵。有一个定理告诉我们，\(n\times n\)矩阵\(\bm A\)可对角化当且仅当 \(\bm A\)有\(n\)个线性无关的特征向量。
然后我们考虑如何计算\(\bm A^m\): \(\bm \Phi^{-1}\bm A\bm \Phi=\bm B\)移项可得\(\bm A=\bm\Phi\bm B\bm\Phi^{-1}\), 则\(\bm A^m=(\bm \Phi\bm B\bm \Phi^{-1})^m=\bm \Phi\bm B(\bm \Phi^{-1}\bm \Phi\bm A)^{m-1}\bm\Phi^{-1}=\bm\Phi\bm B^m\bm\Phi^{-1}\), 因此只要求出\(\bm\Phi\bm B^m\bm\Phi^{-1}\)即可。
下面我们考虑如何求\(\bm\Phi\): 有一个定理告诉我们，\(\bm\Phi=\begin{bmatrix}\bm v_1 & \bm v_2 & \bm v_3 & ... & \bm v_n \end{bmatrix}\), 即为\(n\)个列特征向量构成的矩阵。我们可以验证一下这一点: \(\bm A\bm \Phi=\bm A \begin{bmatrix} \bm v_1 & \bm v_2 & \bm v_3 & ... & \bm v_n \end{bmatrix}=\begin{bmatrix} \bm A\bm v_1 & \bm A\bm v_2 & \bm A\bm v_3 & ... & \bm A\bm v_n \end{bmatrix}=\begin{bmatrix}\lambda_1 \bm v_1 & \lambda_2 \bm v_2 & \lambda_3 \bm v_3 & ... & \lambda_n \bm v_n\end{bmatrix}=\bm \Phi \rm diag(\lambda_1 ,\lambda_2 , \lambda_3 , ... , \lambda_n)\), 从而\(\bm\Phi^{-1}\bm A\bm\Phi=\rm diag(\lambda_1,\lambda_2,\lambda_3,...,\lambda_n)\).
于是，在这道题中，我们构造出桥接矩阵\(\bm \Phi\)和它的逆矩阵，然后进行计算。
根据上面的推导，\(\bm \Phi_{i,j}=(-1)^{i+j} {j\choose i}\), 那如何求出\(\bm \Phi^{-1}\)呢？我们发现其实\(\bm \Phi^{-1}\)就是\(\bm\Phi\)的每一项取绝对值的结果, \(\bm\Phi^{-1}_{i,j}={j\choose i}.\)\[\bm\Phi^{-1}=\begin{bmatrix} 1 & 1 & 1 & 1 & ... & {n\choose 0} \\ 0 & 1 & 2 & 3 & ... & {n\choose1} & \\ 0 & 0 & 1 & 3 & ... & {n\choose 2} \\ & && & ... \\ 0 & 0 & 0 & 0 & ... &{n\choose n}\end{bmatrix}\] 我们可以用二项式反演推出这一点。\[\sum^{n}_{k=0}\bm \Phi^{-1}_{i,k}\bm\Phi_{k,j}=\sum^{n}_{k=0}\bm\Phi^{-1}_{i,k}(-1)^{k+j}{j\choose k}\\ =\sum^n_{k=0}\bm\Phi^{-1}_{i,k}(-1)^{j-k}{j\choose k}=[i=j]\\ \bm\Phi^{-1}_{i,j}=\sum^{n}_{k=0}[i=k]{j\choose k}={j\choose k}\]
最后一步用到了二项式反演。
就这样，我们求出了该矩阵的桥接矩阵\(\bm\Phi\)及其逆矩阵\(\bm\Phi^{-1}\), 这样我们的问题转化为了将一个矩阵\(\bm f_0\)依次左乘\(\bm\Phi^{-1}, \bm B^m, \bm\Phi\)得到\(\bm f_m\).暴力做是\(O(n^2)\)的。

四、最后的优化——从\(O(n^2)\)到\(O(n\log n)\)
先考虑如何快速左乘\(\bm B^m\). 显然因为是对角矩阵，直接\(O(n)\)做即可。
然后考虑如何快速左乘\(\bm \Phi^{-1}\): 推一发设\(b_k=\sum^{n}_{i=k} {i\choose k}a_i=\sum^{n}_{i=k}\frac{i!}{k!(i-k)!}a_i\), 则\(k!b_k=\sum^{n}_{i=k}\frac{i!a_i}{(i-k)!}\). 令\(\alpha_k=k!a_k, \beta_k=k!b_k, \gamma_k=k!\), \(\beta_k=\sum^{n}_{i=k}\alpha_i\gamma_{i-k}\)把数组\(\beta\)和\(\alpha\)倒过来可得\(\beta_{n-k}=\sum^{n}_{i=k}\alpha_{n-i}\gamma_{i-k}=\sum_{i+j=n-k}\alpha_i\gamma_j\). 看出来了吧？是个卷积。愉快地使用FFT解决，\(O(n\log n)\).
快速左乘\(\bm \Phi\)同理\(k!b_k=\sum^{n}_{i=0}\frac{(-1)^{i-k}}{(i-k)!}(i!a_i)\)。总时间复杂度\(O(n\log n)\).
问题至此解决！

代码实现

(说了这么多其实就是fft一下)

//Wrong Coding:
//FFT to (dgr<<1)
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#define llong long long
#define modinc(x) {if(x>=P) x-=P;}
using namespace std;
const int N = 1<<19;
const int LGN = 19;
const int P = 998244353;
const int G = 3;
llong fact[N+3];
llong finv[N+3];
llong a[N+3];
llong b[N+3];
llong c[N+3];
llong d[N+3];
llong e[N+3];
llong f[N+3];
llong g[N+3];
llong fftid[N+3];
llong tmp1[N+3],tmp2[N+3];
int n; llong m;
llong quickpow(llong x,llong y)
{llong cur = x,ret = 1ll;for(int i=0; y; i++){if(y&(1ll<<i)) {ret = ret*cur%P; y-=(1ll<<i);}cur = cur*cur%P;}return ret;
}
llong mulinv(llong x) {return x<=N ? finv[x]*fact[x-1]%P : quickpow(x,P-2);}
void init_fftid(int dgr)
{int len = 0; for(int i=0; i<=LGN; i++) if(dgr==(1<<i)) {len = i; break;}fftid[0] = 0ll;for(int i=1; i<dgr; i++) fftid[i] = ((fftid[i>>1])>>1)|((i&1)<<(len-1));
}
int getdgr(int x)
{int ret = 1; while(ret<x) ret<<=1;return ret;
}
void ntt(int dgr,int coe,llong poly[],llong ret[])
{init_fftid(dgr);for(int i=0; i<dgr; i++) ret[i] = poly[i];for(int i=0; i<dgr; i++) if(i<fftid[i]) swap(ret[i],ret[fftid[i]]);for(int i=1; i<=(dgr>>1); i<<=1){llong tmp = quickpow(G,(P-1)/(i<<1));if(coe==-1) tmp = mulinv(tmp);for(int j=0; j<dgr; j+=(i<<1)){llong expn = 1ll;for(int k=0; k<i; k++){llong x = ret[k+j],y = ret[k+i+j]*expn%P;ret[k+j] = x+y; modinc(ret[k+j]);ret[k+i+j] = x-y+P; modinc(ret[k+i+j]);expn = expn*tmp%P;}}}if(coe==-1){llong tmp = mulinv(dgr);for(int i=0; i<dgr; i++) ret[i] = ret[i]*tmp%P;}
}
int main()
{fact[0] = 1ll;for(int i=1; i<=N; i++) fact[i] = fact[i-1]*i%P;finv[N] = quickpow(fact[N],P-2);for(int i=N-1; i>=0; i--) finv[i] = finv[i+1]*(i+1)%P;scanf("%d%I64d",&n,&m);for(int i=0; i<=n; i++) scanf("%I64d",&a[i]);for(int i=0; i<=n; i++){b[i] = mulinv(i+1);b[i] = quickpow(b[i],m);}int _dgr = n+1,dgr = getdgr(_dgr);for(int i=0; i<_dgr; i++) d[i] = fact[i]*a[i]%P;for(int i=0; i<_dgr-1-i; i++) swap(d[i],d[_dgr-1-i]);for(int i=0; i<_dgr; i++) e[i] = finv[i];for(int i=0; i<_dgr; i++) f[i] = (i&1)==0 ? finv[i] : (P-finv[i])%P;//calculate C = INVPHI*Antt(dgr<<1,1,d,tmp1); ntt(dgr<<1,1,e,tmp2);for(int i=0; i<(dgr<<1); i++) tmp1[i] = tmp1[i]*tmp2[i]%P;ntt(dgr<<1,-1,tmp1,c);for(int i=_dgr; i<(dgr<<1); i++) c[i] = 0ll;for(int i=0; i<_dgr-1-i; i++) swap(c[i],c[_dgr-1-i]);for(int i=0; i<dgr; i++) c[i] = c[i]*finv[i]%P;//calculate C = B*Cfor(int i=0; i<_dgr; i++) c[i] = c[i]*b[i]%P;//calculate G = PHI*Cfor(int i=0; i<_dgr; i++) d[i] = fact[i]*c[i]%P;for(int i=0; i<_dgr-1-i; i++) swap(d[i],d[_dgr-1-i]);ntt(dgr<<1,1,d,tmp1); ntt(dgr<<1,1,f,tmp2);for(int i=0; i<(dgr<<1); i++) tmp1[i] = tmp1[i]*tmp2[i]%P;ntt(dgr<<1,-1,tmp1,g);for(int i=0; i<_dgr-1-i; i++) swap(g[i],g[_dgr-1-i]);for(int i=0; i<dgr; i++) g[i] = g[i]*finv[i]%P;for(int i=0; i<=n; i++) printf("%I64d ",g[i]);return 0;
}

后事
——卧槽这E怎么这么难……
——其实这场还有F。