0.前言
矩阵加速是一种在 OI 中常用的 trick ，经常用于优化部分递推式简单且固定的 dp 式子。
另外，确保你在阅读之前对线性代数有一定的了解。
1.前置知识
矩阵是一个二维数表，可以进行加减乘等运算，并且有一些特殊的性质。
矩阵加法很简单，直接将对应位置加起来即可，此处不再赘述（而且接下来也用不到）。
矩阵乘法比较复杂，这里直接给出公式（设两个相乘的矩阵为 (A,B) ，结果矩阵为 (C) ， (p) 为矩阵 (B) 的行数）： (C_{i,j}=sum_{k=1}^{p}A_{i,k} imes B_{k,j}) 。
不难看出，两个矩阵能够相乘的条件是第一个矩阵的列数等于第二个矩阵的行数，但由于矩阵加速用到的绝大部分都是方阵（行数等于列数），所以我们接下来先只研究方阵。
下面给出了两个二阶方阵相乘的例子：
我们知道实数乘法里有一个单位元 (1) ，在矩阵乘法中也有这样的单位元。我们规定一个单位矩阵 (E) ，使得 (E_{i,j}=[i=j]) ，也就是说只有主对角线上的数是 (1) 。
容易发现单位矩阵一定是一个方阵。下面给出了一个三阶单位矩阵：
这时我们不仅要想，矩阵的乘法定义出来了，那它有没有和实数乘法类似的性质呢？
如果有心的读者将这刚刚那两个矩阵调换过来相乘，结果是：
竟然和上面的不一样！所以我们就可以得出一条结论：矩阵相乘不一定满足交换律。
虽然矩阵乘法没有交换律，但它是拥有结合律的（读者自证不难），于是我们可以方便地计算矩阵的幂 (A^n) ：只需要用类似快速幂的方式，将整数乘法换成矩阵乘法即可。
由此，我们可以开始尝试用矩阵来加速解决一些问题了。
2.例题
矩阵加速的题一般是先设计出一个递推状态，然后将其用矩阵实现。
洛谷P1939 【模板】矩阵加速
递推公式已经给出，下面我们考虑如何将其变为矩阵。
发现这个递推公式用到了前三项的值，我们可以把相邻的三项写成一个向量，每次向前转移一个数。
下面我们考虑转移矩阵应该如何设计，也就是说，在左边的矩阵中填上数使得等式成立：
我们可以考虑右边矩阵里每个数从何而来，结合矩阵乘法的定义，易得转移矩阵为：
举个例子：第一行为 ([1 0 1]) ，根据矩阵乘法的定义，它代表着 (a_n=a_{n-1} imes1+a_{n-2} imes0+a_{n-3} imes1) ，也就是题目中给出的式子，其余同理。
请读者暂停一下，仔细思考每个数由何得出，这对接下来的学习至关重要。
理解了这个式子，解法也就呼之欲出：我们只需要求出来这个转移矩阵的 (n-3) 次方，然后将其乘上初始矩阵即可（这里可以直接手动模拟）。也就是：
这样这个题就做完了。
矩阵加速还经常被用于求两点之间长度为一个值的路径条数，此时需对图的邻接矩阵做快速幂，具体原因及方法下面将会提及。
洛谷P4159 [SCOI2009]迷路
这道题让我们求带权有向图中长度为 (t) 的路径条数，我们先考虑无权图中怎么做。
我们可以设计出一个 dp 式子：用 (f[i][j][t]) 代表从 (i) 到 (j) 长度为 (t) 的条数，那么 (f[i][j][t]=sum_{k=1}^n f[i][k][t-1] imes f[k][j][1]) 。
仔细看的话会发现这个式子很像我们的矩阵乘法，而且它的初始状态实际上就是邻接矩阵。所以此时我们可以把邻接矩阵进行快速幂，结果矩阵中的数就代表路径条数。
回到这个题，这个题的路径是带权的，所以不能直接进行矩阵快速幂（读者自证不难），但我们发现邻接矩阵是以字符串的形式给出，也就是边的长度不超过⑨。
这启发我们把边拆开。很容易想到这样一个思路：将整个图分为⑨层，每层的点先向下一层对应点连边，接下来如果有一条边 ((u,v,w)) ，我们就在 (u_1) 和 (v_w) 之间连边。
此时，如果我们想要从 (u) 走到 (v) ，就相当于从 (u_w) 先向下走 (w-1) 步，再向上走到 (v_w) 。可以证明，在这个⑨层的图里跑，答案是不会变的。
那么此时就可以写出代码了：
const int mod=2009;
int n,t;
char mp[110][110];
struct Matrix {
int a[110][110];
}A,ans;
il void Init(Matrix &A){//相当于构造函数
for(rg int i=1;i<=n*9;i++){
for(rg int j=1;j<=n*9;j++){
A.a[i][j]=0;
}
}
}
il Matrix operator * (Matrix A,Matrix B){//矩阵乘法
Matrix C;
Init(C);
for(rg int i=1;i<=n*9;i++){//前两层循环枚举矩阵A的每个数
for(rg int j=1;j<=n*9;j++){
for(rg int k=1;k<=n*9;k++){//第三层循环枚举矩阵B当前列的数
C.a[i][j]=(C.a[i][j]+A.a[i][k]*B.a[k][j]%mod)%mod;
}
}
}//总复杂度为O(n^3)，所以需要在n比较小的时候用
return C;
}
il Matrix pw(Matrix A,int b){//矩阵快速幂
Matrix ans;
Init(ans);
for(rg int i=1;i<=n*9;i++){//赋为单位矩阵
ans.a[i][i]=1;
}
while(b){//和普通快速幂一样
if(b&1)ans=ans*A;
A=A*A,b>>=1;
}
return ans;
}
il int Idx(int i,int j){//给每层的点都编上号
return (i-1)*9+j;
}
int main(){
Read(n),Read(t);
Init(A);
for(rg int i=1;i<=n;i++){
scanf("%s",mp[i]+1);
for(rg int j=2;j<=9;j++)A.a[Idx(i,j)][Idx(i,j-1)]=1;//每个点先向下连边
for(rg int j=1;j<=n;j++){
if(mp[i][j]=='0')continue;
A.a[Idx(i,1)][Idx(j,mp[i][j]-48)]=1;
}
}
cout<<pw(A,t).a[Idx(1,1)][Idx(n,1)]<<endl;//最终输出的是1到n的答案
return 0;
}
3.总结
矩阵加速是一种好用且重要的优化方法，是提高/省选水平的选手必须掌握的。但需要注意的是，它也有一定的局限性，例如矩阵的大小不能过大、转移必须是固定的形式等等。
4.习题
洛谷P3390 【模板】矩阵快速幂
洛谷P1306 斐波那契公约数
洛谷P2151 [SDOI2009]HH去散步
（选做）洛谷P4719 【模板】"动态 DP"&动态树分治