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

Posted 2021-02-05 suncongbo

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

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

四、最后的优化——从(O(n^2))到(O(nlog n))
先考虑如何快速左乘(m B^m). 显然因为是对角矩阵，直接(O(n))做即可。
然后考虑如何快速左乘(m Phi^{-1}): 推一发设(b_k=sum^{n}_{i=k} {ichoose 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, eta_k=k!b_k, gamma_k=k!), (eta_k=sum^{n}_{i=k}alpha_igamma_{i-k})把数组(eta)和(alpha)倒过来可得(eta_{n-k}=sum^{n}_{i=k}alpha_{n-i}gamma_{i-k}=sum_{i+j=n-k}alpha_igamma_j). 看出来了吧？是个卷积。愉快地使用FFT解决，(O(nlog n)).
快速左乘(m Phi)同理(k!b_k=sum^{n}_{i=0}frac{(-1)^{i-k}}{(i-k)!}(i!a_i))。总时间复杂度(O(nlog 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。