矩阵乘法
考虑一个 \(n\times m\)（即 n 行 m 列）的矩阵乘上一个 \(m\times k\) 的矩阵，乘法后得到 \(n\times k\) 的矩阵。
代数的写法就是
\[C_{i,j}=\sum_{t=1}^m A_{it}\cdot B_{tj}
\]
在写的时候，先枚举 \(i,k\)，如果 \(A_{ik}\) 为 \(0\) 的话就直接跳过这一维能够部分优化常数。
\(C_{i,j}\) 由矩阵 \(A\) 的第 \(i\) 行和矩阵 \(B\) 的第 \(j\) 列的乘积和构成，注意不要搞混了。
注意到矩阵只有结合律而没有交换律。
所以 \(A,B,C\) 是三个矩阵，那么 \(ABC\) 可以是 \(A(BC)\) 但不能是 \(BCA\)。
矩阵快速幂优化计算
把矩阵看作一个数，用对于数做快速幂的方法对于矩阵做快速幂，就是矩阵快速幂了，实现上通常会定义一个矩阵类/结构体，封装一个矩阵乘法。
这样就能快速地求出一个矩阵的 \(n\) 次方了。
矩阵快速幂只适用于方阵，即 \(n\times n\) 的矩阵，并且约定 \(A^0=I\)，\(I\) 为单位矩阵。
P1939
这道题目，要求求出一个数列的第 \(x\) 项，且 \(x\) 很大。
\[a_x=
\begin{cases}
1&x\in\{1,2,3\}\\
a_{x-1}+a_{x-3}&x\geq 4
\end{cases}
\]
那么我们考虑使用一个矩阵 \(\begin{bmatrix}a_x\\a_{x-1}\\a_{x-2}\end{bmatrix}\) 来描述一个状态。
发现状态间有
\[\begin{bmatrix}
a_x\\
a_{x-1}\\
a_{x-2}
\end{bmatrix}
=
\begin{bmatrix}
1&0&1\\
1&0&0\\
0&1&0\\
\end{bmatrix}
\cdot
\begin{bmatrix}
a_{x-1}\\
a_{x-2}\\
a_{x-3}\\
\end{bmatrix}
\]
可以依照下图理解一下，对应行和对应列的相同颜色的项相乘，然后和就是红色的那个值了。
\[\begin{bmatrix}
\color{red}{a_x}\\
a_{x-1}\\
a_{x-2}
\end{bmatrix}
=
\begin{bmatrix}
\color{blue}1&\color{green}0&\color{orange}1\\
1&0&0\\
0&1&0\\
\end{bmatrix}
\cdot
\begin{bmatrix}
\color{blue}a_{x-1}\\
\color{green}a_{x-2}\\
\color{orange}a_{x-3}\\
\end{bmatrix}
\\
\begin{bmatrix}
a_x\\
\color{red}a_{x-1}\\
a_{x-2}
\end{bmatrix}
=
\begin{bmatrix}
1&0&1\\
\color{blue}1&\color{green}0&\color{orange}0\\
0&1&0\\
\end{bmatrix}
\cdot
\begin{bmatrix}
\color{blue}a_{x-1}\\
\color{green}a_{x-2}\\
\color{orange}a_{x-3}\\
\end{bmatrix}
\\
\begin{bmatrix}
a_x\\
a_{x-1}\\
\color{red}a_{x-2}
\end{bmatrix}
=
\begin{bmatrix}
1&0&1\\
1&0&0\\
\color{blue}0&\color{green}1&\color{orange}0\\
\end{bmatrix}
\cdot
\begin{bmatrix}
\color{blue}a_{x-1}\\
\color{green}a_{x-2}\\
\color{orange}a_{x-3}\\
\end{bmatrix}
\]
然后我们发现
\[\begin{bmatrix}
a_x\\
a_{x-1}\\
a_{x-2}
\end{bmatrix}
=\begin{bmatrix}
1&0&1\\
1&0&0\\
0&1&0\\
\end{bmatrix}
\cdot
\begin{bmatrix}
a_{x-1}\\
a_{x-2}\\
a_{x-3}\\
\end{bmatrix}
\]
\[\begin{bmatrix}
a_{x-1}\\
a_{x-2}\\
a_{x-3}
\end{bmatrix}
=\begin{bmatrix}
1&0&1\\
1&0&0\\
0&1&0\\
\end{bmatrix}
\cdot
\begin{bmatrix}
a_{x-2}\\
a_{x-3}\\
a_{x-4}\\
\end{bmatrix}
\]
于是考虑递归地将下面的式子带入到上面的式子，就有了
\[\begin{bmatrix}
a_x\\
a_{x-1}\\
a_{x-2}
\end{bmatrix}
=\begin{bmatrix}
1&0&1\\
1&0&0\\
0&1&0\\
\end{bmatrix}^{x-3}
\cdot
\begin{bmatrix}
a_{3}\\
a_{2}\\
a_{1}\\
\end{bmatrix}
\]
发现这个是可以用矩阵快速幂快速计算的，并且比较容易。
我们把做快速幂的这个 \(3\times 3\) 的矩阵称作转移矩阵。
有些时候，在这个 \(A^{x-3}\) 的形式算完之后，因为 $$\begin{bmatrix}a_3\a_2\a_1\end{bmatrix}$$ 这个矩阵的形式比较简单，可以手动把这一步算出来，然后用 \(A\) 中的项直接表示出这次矩阵乘法的结果，代码会好写一些。
比如这道题目，可以计算出 \(A^{x-2}\) 的值，发现那么我们要求的 \(a_x\) 其实对应的就是目标矩阵的第二行那个数字（原图中 \(a_{x-1}\) 的位置），也就是 \(A_{2,1}\)，就是这个转移矩阵的第二行的第一个数字。
关于写法的问题
可以发现其实上面的转移式还有一种写法
\[\begin{bmatrix}
a_x&a_{x-1}&a_{x-2}
\end{bmatrix}
=
\begin{bmatrix}
a_{x-1}&a_{x-2}&a_{x-3}
\end{bmatrix}
\cdot
\begin{bmatrix}
1&1&0\\
0&0&1\\
1&0&0
\end{bmatrix}
\]
这样计算出来的意思是一样的。
（注意看这两种写法的转移矩阵是不一样的，准确来说做了个转置的操作。）
特别注意，在线段树维护矩阵的时候，两种写法的 pushup 操作是顺序不同的，后面会具体进行说明。
矩阵快速幂优化dp
那么我们考虑在大多数的 dp 中，转移的式子会写成这种形式
\(f_x=a_{x-1}f_{x-1}+a_{x-2}f_{x-2}+\dots +a_{x-k}f_{x-k}\)
这里的 \(a_{x-1}\) 这些系数对于每个 \(f_x\) 都是相同的，也就是说这些都是常数。
也可以说这个时候 \(f_x\) 为线性递推数列。
下面这句话不影响理解文章，可以选择跳过。
当然线性递推数列可以通过特征方程的方式求出通项公式，但是通向公式通常会带有根号导致无法直接计算，有时还有不可用根式表达的情况（Abel–Ruffini），还需要手动写出根号类运算支持根号的运算，也并不简单，并且在 mod 意义下（正确性存疑）。
考虑应用刚才的知识，类比方法
\[\begin{bmatrix}
f_x\\
f_{x-1}\\
f_{x-2}\\
\vdots\\
f_{x-k+1}
\end{bmatrix}
=
\begin{bmatrix}
a_{x-1}&a_{x-2}&a_{x-3}&\cdots&a_{x-k+1}&a_{x-k}\\
1&0&0&\cdots&0&0\\
0&1&0&\cdots&0&0\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\
0&0&0&\cdots&1&0
\end{bmatrix}
\cdot
\begin{bmatrix}
f_{x-1}\\
f_{x-2}\\
f_{x-3}\\
\vdots\\
f_{x-k}
\end{bmatrix}
\]
发现这个转移可以用这种形式来表示。
其实就是把原来 dp 转移的系数平推在第一行了。
但是并不是所有的 dp 的转移形式都长成这个样子，还有另外一种形式。
比如说有些 dp 的转移函数长成这个样子
\[dp_{i,j,k}=\sum a_{i-1,x,y}dp_{i-1,x,y}
\]
\(j,k\) 这两位都不大。
那么这个时候就要把 \(j,k\) 压到一行上，然后做矩阵快速幂。
以 \(j\in\{0,1\}，k\in\{0,1\}\) 举例子。
为了方便表示，下面记录 \(dp_{i,x,y}\) 中 \(dp_{i-1,j,k}\) 的贡献为 \(f_{x,y,j,k}\)，即 \(dp_{i,x,y}=\sum f_{x,y,j,k}dp_{i-1,j,k}\)。
\[\begin{bmatrix}
dp_{i,0,0}\\
dp_{i,0,1}\\
dp_{i,1,0}\\
dp_{i,1,1}\\
\end{bmatrix}
=
\begin{bmatrix}
f_{0,0,0,0}&f_{0,0,0,1}&f_{0,0,1,0}&f_{0,0,1,1}\\
f_{0,1,0,0}&f_{0,1,0,1}&f_{0,1,1,0}&f_{0,1,1,1}\\
f_{1,0,0,0}&f_{1,0,0,1}&f_{1,0,1,0}&f_{1,0,1,1}\\
f_{1,1,0,0}&f_{1,1,0,1}&f_{1,1,1,0}&f_{1,1,1,1}\\
\end{bmatrix}
\cdot
\begin{bmatrix}
dp_{i+1,0,0}\\
dp_{i+1,0,1}\\
dp_{i+1,1,0}\\
dp_{i+1,1,1}\\
\end{bmatrix}
\]
这个转移显然是成立的，注意一下和上面的矩阵是有区别的，实现的时候通常写一个函数 id(j,k) 确定矩阵上的位置。
这个题，就是 \(id(j,k)=2j+k，id(j,k)\in\{0,1,2,3\}\)。
也就是 \((0,0)\) 对应这个矩阵的第一个位置，\((0,1)\) 对应这个矩阵的第二个位置，\((1,0)\) 对应这个矩阵的第三个位置...。
例子是 P5303 和 P5371。
注意到有些时候可能会有更多维度，矩阵乘法都压到一行上。
考虑矩阵是一个 \(n\times n\) 的，那么时间复杂度就是 \(n^3\log x\)，其中 \(x\) 是要求 \(dp_x\) 的值，这里的 \(n^3\) 使用上文说过的优化方法（跳过 \(a_{i,j}=0\) 的位置）通常跑不满。
注意到有些时候矩阵可能会很大，比如 P5303 我的做法\(n\)是\(40+\) ， P5371 我的做法 \(n\)是\(20+\)，这个时候可以先写一个暴力的 dp ，然后把暴力 dp 的每个状态用上述方法映射到一个 id ，原来暴力 dp 的转移\(dp_{i,j,k,l}+=a\cdot dp_{i-1,j2,k2,l2}\)，写成\(M_{id(j,k,l),id(j2,k2,l2)}+=a\)或者\(M_{id(j2,k2,l2),id(j,k,l)}+=a\)，取决于你的写法是那种，矩阵乘向量还是向量乘矩阵。
这里\(M\)是转移矩阵，\(a\)是 dp 时的系数。
线段树维护矩阵乘法
注意到上文所说的矩阵快速幂维护的矩阵都有一个特点，就是 \(dp_i\) 和 \(dp_{i-1}\) 的转移系数 \(a_{i,j}\) 是固定的，并不会因为 \(i\) 的变化而变化。
但是考虑下面这个问题。
P12087
给出一个很大的数（\(10^5\) 位量级），每次询问不大于这个数的不包含子串 13 的数有多少个。
到这里是可以使用树形 dp 进行维护的。
\(dp_{i,j,k}\) 表示处理完了前 \(i\) 位的情况，第 \(j\) 位是不是 \(1\)，是否（\(k\)）卡满了上届，后续一共有多少种填法。
转移可以考虑典型的树形 dp 的转移方法
int LIM;
if(lim)LIM=x[i+1]-'0';
else LIM=9;
for(int j=0;j<=LIM;++j){if(f and j==3)continue;(num+=dfs(i+1,(j==1),lim&(j==(x[i+1]-'0'))))%=mod;
}
return dp[i][f][lim]=num;
但是这个题目还有两个变化。
- 查询区间 \([l,r]\) 直接组成一个数字，这个字问题的答案。
- 修改原来的数字的某一位。
发现这一整个转移过程都可以使用上文提到的以 \(j\in\{0,1\}，k\in\{0,1\}\) 举例子的部分用矩阵的方式表示出转移。
注意到转移和 \(a\) 的值有关，因此考虑使用线段树进行维护，每个节点都维护一个转移矩阵，线段树的性质可以支持单点修改矩阵和区间查询矩阵。
注意一个重要的细节
如果写法是
\[\begin{bmatrix}
dp_{i,0,0}\\
dp_{i,0,1}\\
dp_{i,1,0}\\
dp_{i,1,1}\\
\end{bmatrix}
=
\begin{bmatrix}
f_{0,0,0,0}&f_{0,0,0,1}&f_{0,0,1,0}&f_{0,0,1,1}\\
f_{0,1,0,0}&f_{0,1,0,1}&f_{0,1,1,0}&f_{0,1,1,1}\\
f_{1,0,0,0}&f_{1,0,0,1}&f_{1,0,1,0}&f_{1,0,1,1}\\
f_{1,1,0,0}&f_{1,1,0,1}&f_{1,1,1,0}&f_{1,1,1,1}\\
\end{bmatrix}
\cdot
\begin{bmatrix}
dp_{i+1,0,0}\\
dp_{i+1,0,1}\\
dp_{i+1,1,0}\\
dp_{i+1,1,1}\\
\end{bmatrix}
\]
记录 \(x_n\) 为转移矩阵。
\[f_n=x_n f_{n-1}\\
f_{n-1}=x_{n-1} f_{n-2}\\
f_n=x_nx_{n-1}\cdots x_2 f_1
\]
你会发现这里的矩阵是从大下标一直乘到小下标的，而矩阵乘法是不具有交换律的，因此写成 pushup 的时候要记得写右子树矩阵乘左子树矩阵。
不过如果转移矩阵写成
\[\begin{bmatrix}
dp_{i,0,0}&dp_{i,0,1}&dp_{i,1,0}&dp_{i,1,1}
\end{bmatrix}
=
\begin{bmatrix}
dp_{i-1,0,0}&dp_{i-1,0,1}&dp_{i-1,1,0}&dp_{i-1,1,1}
\end{bmatrix}
\cdot
\begin{bmatrix}
f_{0,0,0,0}&f_{0,1,0,0}&f_{1,0,0,0}&f_{1,1,0,0}\\
f_{0,0,0,1}&f_{0,1,0,1}&f_{1,0,0,1}&f_{1,1,0,1}\\
f_{0,0,1,0}&f_{0,1,1,0}&f_{1,0,1,0}&f_{1,1,1,0}\\
f_{0,0,1,1}&f_{0,1,1,1}&f_{1,0,1,1}&f_{1,1,1,1}\\
\end{bmatrix}
\]
就不存在这种问题了，展开递推式子是按照小下标到大下标的顺序的，线段树直接按照 \(ls*rs\) 就是对的，在使用线段树维护矩阵的时候推荐这张写法。
线段树维护链上问题的矩阵乘法
有些时候，树上查询链的问题，将链看作一个数组，也能按照上文的方法进行矩阵优化。
那么就有一个思路，把树做一个树剖，将重链看作一个正常的数列用线段树维护，就能解决链上问题了。
注意链上问题要特别注意顺序，是向上走还是向下走，因为矩阵没有交换律。
例题，P14468，代码细节或许比较多。