Codeforces 947E Perpetual Subtraction (线性代数、矩阵对角化、DP)

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

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

四、最后的优化——从$O(n^2)$到$O(n\log n)$
先考虑如何快速左乘${\mathbf{B}}^m$. 显然因为是对角矩阵，直接$O(n)$做即可。
然后考虑如何快速左乘${\mathbf{\Phi }}^{-1}$: 推一发设$b_k={\sum }_{i=k}^n\left(\genfrac{}{}{0ex}{}{i}{k}\right)a_i={\sum }_{i=k}^n\frac{i!}{k!(i-k)!}{a}_i$, 则$k!b_k={\sum }_{i=k}^n\frac{i!a_i}{(i-k)!}$. 令${\alpha }_k=k!a_k,{\beta }_k=k!b_k,{\gamma }_k=k!$, ${\beta }_k={\sum }_{i=k}^n{\alpha }_i{\gamma }_{i-k}$把数组$\beta$和$\alpha$倒过来可得${\beta }_{n-k}={\sum }_{i=k}^n{\alpha }_{n-i}{\gamma }_{i-k}={\sum }_{i+j=n-k}{\alpha }_i{\gamma }_j$. 看出来了吧？是个卷积。愉快地使用FFT解决，$O(n\log n)$.
快速左乘$\mathbf{\Phi }$同理$k!b_k={\sum }_{i=0}^n\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*A
 ntt(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*C
 for(int i=0; i<_dgr; i++) c[i] = c[i]*b[i]%P;
 //calculate G = PHI*C
 for(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。