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

最新推荐文章于 2022-06-15 12:02:18 发布

suncongbo

最新推荐文章于 2022-06-15 12:02:18 发布

阅读量612

点赞数

分类专栏： codeforces 多项式线性代数

本文链接：https://blog.csdn.net/suncongbo/article/details/84962727

版权

多项式同时被 3 个专栏收录

7 篇文章 0 订阅

订阅专栏

codeforces

2 篇文章 0 订阅

订阅专栏

线性代数

1 篇文章 0 订阅

订阅专栏

呜啊怎么又是数学了啊。。。数学比例 $\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 (n m)$ 的 $d p$ 很好想: 设 $f [i] [j]$ 表示 $i$ 轮后 $x = j$ 的概率，则轻易列出方程: $f[i][j]=\sum^{n}_{k=i} \frac{dp[i-1][k]}{k+1}$ .
按此 $d p$ 时间复杂度 $O (n m)$ , 太慢了。考虑优化。
写出转移矩阵: 令 $(n + 1)$ 维列向量为 $d p [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*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。

suncongbo

关注

0
点赞
踩
1

收藏

觉得还不错? 一键收藏
0
评论
Codeforces 947E Perpetual Subtraction (线性代数、矩阵对角化、DP)

呜啊怎么又是数学了啊。。。数学比例1633=0.4848\frac{16}{33}=0.48483316=0.4848题目链接: https://codeforces.com/contest/947/problem/E题意:有一个[0,n][0,n][0,n]的随机数xxx初始为iii的概率为pip_ipi. mmm次操作每次从[0,x][0,x][0,x]中等概率随机选择一个数yyy把...
复制链接

扫一扫