下面的内容不用看了，直接看 「算法笔记」线性代数 中的“矩阵加速递推”部分吧 qwq
一、矩阵定义
在数学中，矩阵是一个按照长方阵列排列的复数或实数集合。一个 n×m 的矩阵一般这样表示：
(A_{n,m}=egin{bmatrix}a_{1,1}&a_{1,2}&cdots&a_{1,m}\a_{2,1}&a_{2,2}&cdots&a_{2,m}\vdots&vdots&ddots&vdots\a_{n,1}&a_{n,2}&cdots&a_{n,m}end{bmatrix})
所谓方阵，就是一个 n×n 的矩阵，也称作 n 阶矩阵。由于行数和列数相等，所以方阵的幂是有意义的。
方阵的幂是指，A 是一个方阵，将 A 连乘 n 次，即：C=An。若不是方阵则无法进行乘幂运算。
二、矩阵的乘法运算
设 A,B 是两个矩阵，令 C=A×B。
如果 A 是 n×p 的矩阵，B是 m×p 的矩阵，A 和 B 的乘积 C 是一个 n×m 的矩阵。
(C_{i,j}=A_{i,1} imes B_{1,j} + A_{i,2} imes B_{2,j} + A_{i,3} imes B_{3,j} + cdots +A_{i,p} imes B_{p,j} = sum_{k=1}^p A_{i,k} imes B_{k,j})
即乘积 C 的第 i 行第 j 列的元素 Ci,j 等于矩阵 A 的第 i 行的元素与矩阵 B 的第 j 列对应元素乘积之和。
(egin{bmatrix}1&4\2&5\3&6end{bmatrix} imes egin{bmatrix}1&2&3\4&5&6end{bmatrix}=egin{bmatrix}1 imes 1+4 imes4&1 imes2+4 imes5&1 imes3+4 imes6\2 imes1+5 imes4&2 imes2+5 imes5&2 imes3+5 imes6\3 imes1+6 imes4&3 imes2+6 imes5&3 imes4+6 imes6end{bmatrix}=egin{bmatrix}17&22&27\22&29&36\27&36&45end{bmatrix})
题目链接 一个 n×m 的矩阵 A 与一个 m×p 的矩阵 B 相乘得到一个 n×p 的矩阵 C。给定矩阵 A 和矩阵 B，求矩阵 C。这里放一个代码片段。完整代码
memset(c,0,sizeof(c)); for(int i=0;i<n;i++) for(int j=0;j<p;j++) for(int k=0;k<m;k++) c[i][j]=c[i][j]+a[i][k]*b[k][j];
三、矩阵乘法结合律
矩阵乘法不满足交换律，但是满足结合律。
若 (A×B×C=D)，设 (A×B=T)，则
(D_{i,j}=sum_{x=1}^q T_{i,x} imes C_{x,j})
(D_{i,j}=sum_{x=1}^q (sum_{y=1}^p A_{i,y} imes B_{y,x}) imes C_{x,j})
(D_{i,j}=sum_{x=1}^q sum_{y=1}^p A_{i,y} imes B_{y,x} imes C_{x,j})
设 (A×(B×C)=E)，同理可得，
(E_{i,j}=sum_{x=1}^p sum_{y=1}^q A_{i,x} imes B_{x,y} imes C_{y,j})
故 (D=E)。
四、矩乘优化线性递推
1. POJ 3070 斐波那契数列第n项
题目链接 在斐波那契数列中，(F_0=1,F_1=1,F_n=F_{n-1}+F_{n-2})（(n>1)）。现给定整数 (n,m)，(0≤n≤2 imes 10^9,m=10000)，求 (F_n mod m)。
设 (f(n)) 是一个 1×2 的矩阵，(f(n)=egin{bmatrix}F_n&F_{n+1 }end{bmatrix})。我们希望能根据 (f(n-1)=egin{bmatrix}F_{n-1}&F_nend{bmatrix}) 计算出 (f(n))。把 (f(n-1)) 第二列上的数作为 (f(n)) 第一列上的数，并把 (f(n-1)) 第 1 列和第 2 列的数都累加到 (f(n)) 的第二列上。
(F_n=0 imes F_{n-1}+1 imes F_n)，(F_{n+1}=1 imes F_{n-1}+1 imes F_n)。
所以我们令转移矩阵 A 为 (egin{bmatrix}0&1\1&1 end{bmatrix})。(f(n)=f(n-1) imes A)。
(egin{bmatrix}F_{n-1}&F_nend{bmatrix} imes egin{bmatrix}0&1\1&1 end{bmatrix}=egin{bmatrix}F_n&F_{n-1}+F_nend{bmatrix}=egin{bmatrix}F_n&F_{n+1 }end{bmatrix}=egin{bmatrix}F_0&F_1end{bmatrix} imes egin{bmatrix}0&1\1&1 end{bmatrix}^n)
则初值为 (f(0)=egin{bmatrix}F_0&F_1end{bmatrix}=egin{bmatrix}0&1end{bmatrix})，目标为 (f(n)=f(0) imes A^n)。
由于矩阵乘法满足结合律，所以我们可以通过与整数快速幂相同的方法计算 (f(0) imes A^n)，得到矩阵的第 1 列上的数就是 (F_n)。取模只要直接在执行矩阵乘法时，对各个整数的加法、乘法运算加入取模操作。
#include<cstdio> #include<cstring> #define int long long using namespace std; const int N=2,mod=10000; int n,f[N][N]; void mul(int x[N][N],int y[N][N]){ int c[N][N]; memset(c,0,sizeof(c)); for(int i=0;i<N;i++) for(int j=0;j<N;j++) for(int k=0;k<N;k++) c[i][j]=(c[i][j]+x[i][k]*y[k][j])%mod; memcpy(x,c,sizeof(c)); } signed main(){ while(~scanf("%lld",&n)&&n!=-1){ int a[N][N]={{0,1},{1,1}}; memset(f,0,sizeof(f)); for(int i=0;i<N;i++) f[i][i]=1; for(;n;n>>=1,mul(a,a)) if(n&1) mul(f,a); printf("%lld ",f[0][1]); } return 0; }
另一种版本的代码。
(egin{bmatrix}F_{n-1}&F_nend{bmatrix} imes egin{bmatrix}0&1\1&1 end{bmatrix}=egin{bmatrix}F_n&F_{n-1}+F_nend{bmatrix}=egin{bmatrix}F_n&F_{n+1 }end{bmatrix}=egin{bmatrix}F_1&F_2end{bmatrix} imes egin{bmatrix}0&1\1&1 end{bmatrix}^{n+1})
初值为 (egin{bmatrix}F_1&F_2end{bmatrix})。mul 就是求 x*y 得到的矩阵。最后得到矩阵的第一个数就是答案。
顺便补充一下，单位矩阵就是主对角线上的元素都为1，其余元素全为 0 的 n 阶矩阵，称为 n 阶单位矩阵，记为 In 或 En，通常用 I 或 E 表示。在矩阵乘法中，单位矩阵起着特殊的作用，如同数的乘法中的1。
#include<cstdio> #include<cstring> #define int long long using namespace std; const int N=2,mod=10000; int n,f[N][N]; void mul(int x[N][N],int y[N][N]){ int c[N][N]; memset(c,0,sizeof(c)); for(int i=0;i<N;i++) for(int j=0;j<N;j++) for(int k=0;k<N;k++) c[i][j]=(c[i][j]+x[i][k]*y[k][j])%mod; memcpy(x,c,sizeof(c)); } signed main(){ while(~scanf("%lld",&n)&&n!=-1){ int a[N][N]={{0,1},{1,1}}; memset(f,0,sizeof(f)),n++; for(int i=0;i<N;i++) f[i][i]=1; for(;n;n>>=1,mul(a,a)) if(n&1) mul(f,a); printf("%lld ",f[0][0]); } return 0; }
2. Luogu P5110 块速递推
题目链接 给定一个数列 a 满足 (a_0=0,a_1=1,a_n=233a_{n-1}+666a_{n-2})，求 (a_n mod 10^9+7)，多组询问。(n≤10^9,T≤5 imes10^7)。
先考虑朴素的做法，可以直接用矩阵乘法和矩阵快速幂求解。矩乘的式子如下：
(egin{bmatrix}a_n&a_{n-1 }end{bmatrix}=egin{bmatrix}233 imes a_{n-1}+666 imes a_{n-2}\1 imes a_{n-1}+0 imes a_{n-2} end{bmatrix}=egin{bmatrix} 233&666\1&0end{bmatrix} imes egin{bmatrix}a_{n-1}\a_{n-2}end{bmatrix}=egin{bmatrix} 233&666\1&0end{bmatrix}^{n-1} imes egin{bmatrix}a_1\a_0end{bmatrix})
时间复杂度为 O(23T log n)，不能通过此题。
记转移矩阵 (egin{bmatrix} 233&666\1&0end{bmatrix}) 为 A。
要快速求出 Ax，用类似分块的方法，可以令 (t=lfloor sqrt n
floor)，那么可以将任何 x 表示为 ut+v 的形式，其中 (u,v≤O(sqrt n))。(u=lfloor frac{x}{t}
floor,v=x mod t)。
(A^x=(A^t)^u imes A^v)，预处理 ((A^t)^i) 和 (A^i)，就可以 O(1) 求得答案。
3. Luogu P5107 能量采集
题目链接 n 个点 m 条边的无向图，每个点有初始权值 ai。每一个时刻，节点 i 上的权值会平均分成 di+1 份，其中 di 份通过与 i 相邻的边加到与 i 相邻的节点上，剩下 1 份留在 i。求 t 时刻后每个点的权值。不保证无重边自环。
首先在每个点上l连一个自环，那么现在就不存在留在原地的那份权值了。
记 vi,j 为第 i 时刻，节点 j 的权值；记 numi,j 为 (i,j) 之间连的边的数量。
尝试构造 n×n 的矩阵 A，使得
(egin{bmatrix}v_{i-1,1}&v_{i-1,2} cdots v_{i-1,n}end{bmatrix} imes A=egin{bmatrix} v_{i,1}&v_{i,2} cdots v_{i,n} end{bmatrix})
我们有 (v_{i,j}=sum_{x=1}^n v_{i-1,x} imes frac{num_{x,j}}{d_x})。 所以令 (A_{i,j}=frac{num_{i,j}}{d_i}) 即可。
五、一些其他应用
1. 计算本质不同子序列数量
记 ai,j 表示以 i 开头，以 j 结尾的，不作为当前序列子序列的，但是去掉最后一位能作为当前序列子序列的数量。
2. 套线段树/其他数据结构
解决一些诸如动态dp的问题。
由于在证明矩阵乘法的结合律时只用到了加法和乘法的分配律 (a+b)c=ac+bc，所以可以对矩阵乘法进行重定义来解决更多问题。
如我们有 max(a,b)+c=max(a+c,b+c)，于是可以定义 A×B=C 为 (C_{i,j}=max_{k=1}^p A_{i,k}+B_{k,j})。