大家好，我是Weekoder！
今天要讲的内容是矩阵加速！
这时候就有人说了：
\(\tiny{\texttt{Weekoder 这么蒻，怎么会矩阵啊。还给我们讲，真是十恶不赦！}}\)
不不不，容我解释。在经过我的研究后，我发现基本的矩阵运算和矩阵加速都并没有那么难。只要继续往下看，相信你也能学会！
注意：以下内容的学习难度将会用颜色表示，与洛谷题目难度顺序一致，即 \(\color{#FE4C61}\texttt{红}\color{#000000}<\color{#F39C11}\texttt{橙}\color{#000000}<\color{#FFC116}\texttt{黄}\color{#000000}<\color{#52C41A}\texttt{绿}\color{#000000}\)。（并不对标洛谷题目难度，只作为学习难易度参考）
\[\huge\texttt{Part 1}\small\texttt{ Definition} \]\[\color{#FE4C61}\texttt{定义} \]矩阵和二维数组很像，是由 \(m\times n\) 个数排列成 \(m\) 行 \(n\) 列的一张表，由于排列出来的表是一个矩形，故称其为矩阵。矩阵长这个样子：
\[\begin{pmatrix} a_{11} & a_{12} & a_{13} & \cdots & a_{1n} \\ a_{21} & a_{22} & a_{23} & \cdots & a_{2n} \\ a_{31} & a_{32} & a_{33} & \cdots & a_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ a_{m1} & a_{m2} & a_{m3} & \cdots & a_{mn} \end{pmatrix} \]可以看到，矩阵中的每个元素都有着对应的行和列，我们把一个矩阵记作 \(A\)，第 \(i\) 行 \(j\) 列的元素即为 \(a_{ij}\)。更形式化的，写作：
\[A=(a_{ij})\in \mathbb{F^{m\times n}} \]其中 \(\mathbb{F}\) 为数域，一般取为实数域 \(\mathbb{R}\) 或复数域 \(\mathbb{C}\)。（看不懂没事，蒟蒻自行走开QWQ）
\[\huge\texttt{Part 2}\small\texttt{ Special matrices} \]\[\color{#FE4C61}\texttt{特殊矩阵} \]\(\texttt{1.零矩阵}\)
元素全部为 \(0\) 的矩阵称为零矩阵。像这样：
\[\begin{pmatrix} 0 & 0 & \cdots & 0 \\ 0 & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 0 \end{pmatrix} \]零矩阵记作 \(0_{m\times n}\)，就是在 \(0\) 下面加上矩阵的大小 \(m\times n\)。你可以把零矩阵看做数字 \(0\)，任何数乘以 \(0\) 都得 \(0\)。
\(\texttt{2.对角矩阵}\)
只有主对角线上的元素有值，其余元素为 \(0\) 的矩阵称为对角矩阵。
注：主对角线为矩阵中从左上角到右下角的一条对角线。
\[\begin{pmatrix} a_1 & 0 & \cdots & 0 \\ 0 & a_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & a_n \end{pmatrix} \]对角矩阵根据主对角线上的值，记作 \(\text{diag(}\text{a}_1,\text{a}_2,\ldots,\text{a}_n\text{)}\)。
\(\texttt{3.单位矩阵}\)
主对角线上的元素均为 \(1\)，其余元素为 \(0\) 的矩阵称为单位矩阵。
\[\begin{pmatrix} 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{pmatrix} \]单位矩阵记作 \(I\)。
记得分数中的概念分数单位吗？矩阵单位和分数单位的“地位”差不多，代表的都是最基础的，最小的独立个体。你可以把单位矩阵看做数字 \(1\)，任何数乘以 \(1\) 都等于它本身。
最基础的，常见的特殊矩阵就是这些了。当然，还有很多的特殊矩阵，不过我们暂时用不到。
\[\huge\texttt{Part 3}\small\texttt{ Matrix operations} \]\[\color{#F39C11}\texttt{矩阵运算} \]\(\texttt{1.相等}\)
若对于矩阵 \(A,B\)，所有的 \(i,j\) 都有 \(a_{ij}=b_{ij}\) 且矩阵的行和列相等，则称矩阵 \(A,B\) 相等。
其实就是两个矩阵长得一模一样。
\[\begin{pmatrix} a_{11} & a_{12} & a_{13} & \cdots & a_{1n} \\ a_{21} & a_{22} & a_{23} & \cdots & a_{2n} \\ a_{31} & a_{32} & a_{33} & \cdots & a_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ a_{m1} & a_{m2} & a_{m3} & \cdots & a_{mn} \end{pmatrix} = \begin{pmatrix} b_{11} & b_{12} & b_{13} & \cdots & b_{1n} \\ b_{21} & b_{22} & b_{23} & \cdots & b_{2n} \\ b_{31} & b_{32} & b_{33} & \cdots & b_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ b_{m1} & b_{m2} & b_{m3} & \cdots & b_{mn} \end{pmatrix} \]\(\texttt{2.矩阵加（减）法}\)
若要求 \(A,B\) 两个矩阵之和，即 \(C=A+B\)，则对于任意 \(i,j\)，满足 \(c_{ij}=a_{ij}+b_{ij}\)。要求矩阵行列相等。
总结一句话：对应位置相加。
\[\begin{aligned} {} & \begin{pmatrix} a_{11} & a_{12} & a_{13} & \cdots & a_{1n} \\ a_{21} & a_{22} & a_{23} & \cdots & a_{2n} \\ a_{31} & a_{32} & a_{33} & \cdots & a_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ a_{m1} & a_{m2} & a_{m3} & \cdots & a_{mn} \end{pmatrix} + \begin{pmatrix} b_{11} & b_{12} & b_{13} & \cdots & b_{1n} \\ b_{21} & b_{22} & b_{23} & \cdots & b_{2n} \\ b_{31} & b_{32} & b_{33} & \cdots & b_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ b_{m1} & b_{m2} & b_{m3} & \cdots & b_{mn} \end{pmatrix} \\ = & \begin{pmatrix} a_{11}+b_{11} & a_{12}+b_{12} & a_{13}+b_{13} & \cdots & a_{1n}+b_{1n} \\ a_{21}+b_{21} & a_{22}+b_{22} & a_{23}+b_{23} & \cdots & a_{2n}+b_{2n} \\ a_{31}+b_{31} & a_{32}+b_{32} & a_{33}+b_{33} & \cdots & a_{3n}+b_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ a_{m1}+b_{m1} & a_{m2}+b_{m2} & a_{m3}+b_{m3} & \cdots & a_{mn}+b_{mn} \end{pmatrix} \end{aligned} \]矩阵加法满足交换律和结合律：
\[A+B=B+A \]\[(A+B)+C=A+(B+C) \]减法同理，对应位置相减。
\[\begin{aligned} {} & \begin{pmatrix} a_{11} & a_{12} & a_{13} & \cdots & a_{1n} \\ a_{21} & a_{22} & a_{23} & \cdots & a_{2n} \\ a_{31} & a_{32} & a_{33} & \cdots & a_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ a_{m1} & a_{m2} & a_{m3} & \cdots & a_{mn} \end{pmatrix} - \begin{pmatrix} b_{11} & b_{12} & b_{13} & \cdots & b_{1n} \\ b_{21} & b_{22} & b_{23} & \cdots & b_{2n} \\ b_{31} & b_{32} & b_{33} & \cdots & b_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ b_{m1} & b_{m2} & b_{m3} & \cdots & b_{mn} \end{pmatrix} \\ = & \begin{pmatrix} a_{11}-b_{11} & a_{12}-b_{12} & a_{13}-b_{13} & \cdots & a_{1n}-b_{1n} \\ a_{21}-b_{21} & a_{22}-b_{22} & a_{23}-b_{23} & \cdots & a_{2n}-b_{2n} \\ a_{31}-b_{31} & a_{32}-b_{32} & a_{33}-b_{33} & \cdots & a_{3n}-b_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ a_{m1}-b_{m1} & a_{m2}-b_{m2} & a_{m3}-b_{m3} & \cdots & a_{mn}-b_{mn} \end{pmatrix} \end{aligned} \]\(\texttt{3.矩阵数乘}\)
数 \(\lambda\)（一个数字） 乘以矩阵 \(A\)，记作 \(\lambda A\)，即为矩阵数乘运算。若有 \(B=\lambda A\)，则对于任意 \(i,j\) 都满足 \(b_{ij}=\lambda a_{ij}\)。
还是一句话：对应位置相乘。
\[\lambda\begin{aligned} {} & \begin{pmatrix} a_{11} & a_{12} & a_{13} & \cdots & a_{1n} \\ a_{21} & a_{22} & a_{23} & \cdots & a_{2n} \\ a_{31} & a_{32} & a_{33} & \cdots & a_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ a_{m1} & a_{m2} & a_{m3} & \cdots & a_{mn} \end{pmatrix} = \begin{pmatrix} \lambda a_{11} & \lambda a_{12} & \lambda a_{13} & \cdots & \lambda a_{1n} \\ \lambda a_{21} & \lambda a_{22} & \lambda a_{23} & \cdots & \lambda a_{2n} \\ \lambda a_{31} & \lambda a_{32} & \lambda a_{33} & \cdots & \lambda a_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \lambda a_{m1} & \lambda a_{m2} & \lambda a_{m3} & \cdots & \lambda a_{mn} \end{pmatrix} \end{aligned} \]\[\huge\texttt{Part 3.5}\small\texttt{ Matrix multiplication} \]\[\color{#FFC116}\texttt{矩阵乘法} \]虽然矩阵乘法也属于矩阵运算，但难度比前面的都高，而且是今天的重点内容，所以单独放出来讲，故记为 \(\texttt{Part 3.5}\)。（话说你们没有发现难度变成黄了吗）
上例题！（虽然难度是橙）
先看矩阵乘法的定义：若有 \(n\) 行 \(m\) 列矩阵 \(A\) 和 \(m\) 行 \(k\) 列的矩阵 \(B\)（\(A\) 的行与 \(B\) 的列相等），则 \(n\) 行 \(k\) 列的矩阵 \(C=A\times B\) 满足
\[c_{ij}=\sum_{l=1}^m a_{il}\times b_{lj} \]只要枚举 \(i,j\)（范围是 \(n,k\)），并套用公式就能用 \(O(n^3)\) 的时间复杂度解决这个问题。
我知道，这看起来根本不是新手蒟蒻能看懂的。那我就用人话来讲讲矩阵乘法。
矩阵乘法并不是一个一个乘，而是行对应列乘。怎么个乘法呢？我们来看看下面两个矩阵相乘的例子。
\[\begin{pmatrix} 5 & 2 & 3 \\ 7 & 9 & 4 \\ \end{pmatrix} \times \begin{pmatrix} 2 & 6 & 8 & 1 \\ 0 & 9 & 1 & 3 \\ 2 & 4 & 4 & 1 \end{pmatrix} \]第一个矩阵为 \(A\)，第二个矩阵为 \(B\)。
我们先取出 \(A\) 的第一行。像这样：
\[\begin{pmatrix} 5 & 2 & 3 \end{pmatrix} \]再取出 \(B\) 的第一列。像这样：
\[\begin{pmatrix} 2 \\ 0 \\ 2 \end{pmatrix} \]不对，你给我转过来。
\[\begin{pmatrix} 2 & 0 & 2 \end{pmatrix} \]现在终于可以相乘了。逐位相乘得出结果：
\[\begin{pmatrix} 5 & 2 & 3 \end{pmatrix} \times \begin{pmatrix} 2 & 0 & 2 \end{pmatrix} = \begin{pmatrix} 5\times2 & 2\times0 & 3\times2 \end{pmatrix} = \begin{pmatrix} 10 & 0 & 6 \end{pmatrix} \]得出了结果 \(\begin{pmatrix} 10 & 0 & 6 \end{pmatrix}\)。再将每一位相加：
\[\begin{pmatrix} 10 & 0 & 6 \end{pmatrix} \to 10+0+6=16 \]还记得我们之前是怎么取的吗？我们取了 \(A\) 的第一行和 \(B\) 的第一列（注意加粗部分），所以答案就存储在 \(C\) 的第一行第一列。还没搞懂？更通用一点：我们取了 \(A\) 的第 \(x\) 行和 \(B\) 的第 \(y\) 列（注意加粗部分），所以答案就存储在 \(C\) 的第 \(x\) 行第 \(y\) 列。也就是说，当我们想要获取矩阵 \(C\) 的第 \(x\) 行 \(y\) 列的时候，就需要取 \(A\) 的第 \(x\) 行和 \(B\) 的第 \(y\) 列，相乘再相加。由于 \(A\) 的行数与 \(B\) 的列数相等，取出来的数列才可以逐位相乘（不然元素个数不一样）。而取出来的数列长度就是 \(m\)，所以可以用 \(O(m)\) 求和，总时间复杂度 \(O(nmk)=O(n^3)\)。
最后，可以看看代码辅助理解。
#include <bits/stdc++.h>
using namespace std;
const int N = 105;
int n, m, k, a[N][N], b[N][N]; // 用二维数组存矩阵 A,B
int main() {
cin >> n >> m >> k;
for (int i = 1; i <= n; i++)
for (int j = 1; j <= m; j++)
cin >> a[i][j]; // 输入矩阵 A
for (int i = 1; i <= m; i++)
for (int j = 1; j <= k; j++)
cin >> b[i][j]; // 输入矩阵 B
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= k; j++) { // 枚举 C 矩阵 n 行 k 列的每个元素
// 以下部分为模拟刚刚讲的矩阵乘法
int sum = 0; // 求和，sum 即为 C_ij
for (int l = 1; l <= m; l++)
sum += a[i][l] * b[l][j]; // 求和，A 的行和 B 的列，建议模拟一下过程加强理解
cout << sum << " "; // 输出 sum（C_ij）
}
cout << "\n"; // 记得换行！
}
return 0; // 完美的结束
}
这样就能愉快地切掉这道题了。请完成这道题再继续！
矩阵乘法满足以下性质：
结合律：\((AB)C=A(BC)\)
分配律：\((A+B)C=AC+BC\)
矩阵乘法不满足交换律。（这是重点！）
有了矩阵乘法，我们还可以结合上面的特殊矩阵得到一些性质：
\[A\times I=A \]\[A\times 0_{m\times n}=0_{m\times n} \]\[\huge\texttt{Part 4}\small\texttt{ Matrix fast power} \]\[\color{#52C41A}\texttt{矩阵封装 & 矩阵快速幂} \]快到今天的主题了！上例题！
点开题目后的你 be like：
这是啥呀？
我来让题目描述“缩点水”：
给定一个 \(n\) 行 \(n\) 列的矩阵 \(A\)，求 \(A^k\)，即 \(\underbrace{A\times A\times A\times\cdots\times A\times A}_{k\texttt{ 次}}\)。
第一思路：暴力！直接做 \(k\) 次矩阵乘法，时间复杂度 \(O(kn^3)\)。看看数据范围：
\(0\le k\le10^{12}\)
考虑放弃做题。
那我们该怎么优化呢？看到需要计算 \(A^k\)，我突然想到了一个算法：快速幂！但是矩阵快速幂该怎么写呢？答案是：和正常的快速幂一样，矩阵也能使用快速幂，只不过快速幂中的乘法变成了矩阵乘法。但是矩阵乘法太难写，有没有什么办法能让矩阵乘法也像普通的乘法一样，只要写一个 * 乘号就行了呢？
注意：不会快速幂的话可以先简单看看我写的文章。
回到主题，有没有什么办法能只要写一个 * 乘号就能进行矩阵乘法呢？其实我们可以用结构体把矩阵封装起来，再用重载运算符就行了。关于重载运算符，可以参考这些资料。
定义一个矩阵类型的结构体可以写成这样：
struct Matrix {
};
我们需要在里面用一个二维数组存储矩阵。我们还可以写一个结构体初始化函数，只要定义了一个矩阵，就自动清零，免去清零的麻烦。
struct Matrix {
int a[N][N]; // N 为矩阵大小
Matrix() {
memset(a, 0, sizeof a);
}
};
最后，把矩阵乘法写进去。
struct Matrix {
ll a[N][N];
Matrix() {
memset(a, 0, sizeof a);
}
Matrix operator*(const Matrix &x)const {
Matrix res;
for (int i = 1; i <= n; i++)
for (int j = 1; j <= n; j++)
for (int k = 1; k <= n; k++)
res.a[i][j] = (res.a[i][j] % MOD + a[i][k] % MOD * x.a[k][j] % MOD) % MOD;
return res;
}
};
注意，这里一定要写成 a[i][k] * x.a[k][j]，不能写成 x.a[i][k] * a[k][j]，因为矩阵乘法不满足交换律！
这样，结构体封装部分就完成了。
我们要定义两个矩阵：\(a\) 和 \(base\)。\(a\) 是输入的矩阵，\(base\) 是答案矩阵，所以 \(base\) 需要初始化成 \(I\)（单位矩阵），写一个初始化函数 \(\operatorname{init}\)，如下：
void init() {
for (int i = 1; i <= n; i++) base.a[i][i] =1;
}
初始化完以后，就可以执行快速幂了，计算 \(A^k\) 了，让 \(base\) 乘 \(A\)。矩阵快速幂核心代码如下：
void expow(ll b) {
while (b) {
if (b & 1) base = base * a;
a = a * a, b >>= 1;
}
}
有一点需要注意的就是，不能写成 base *= a 等形式，因为重载运算符定义的是 *，没有定义 *=，所以需要将 *= 展开。
最后，就可以输出 \(base\) 了。展示全部代码：
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 105, MOD = 1e9 + 7;
int n;
ll k;
struct Matrix {
ll a[N][N];
Matrix() {
memset(a, 0, sizeof a);
}
Matrix operator*(const Matrix &x)const {
Matrix res;
for (int i = 1; i <= n; i++)
for (int j = 1; j <= n; j++)
for (int k = 1; k <= n; k++)
res.a[i][j] = (res.a[i][j] % MOD + a[i][k] % MOD * x.a[k][j] % MOD) % MOD;
return res;
}
}a, base;
void init() {
for (int i = 1; i <= n; i++) base.a[i][i] =1;
}
void expow(ll b) {
while (b) {
if (b & 1) base = base * a;
a = a * a, b >>= 1;
}
}
int main() {
cin >> n >> k;
for (int i = 1; i <= n; i++)
for (int j = 1; j <= n; j++)
cin >> a.a[i][j];
init();
expow(k);
for (int i = 1; i <= n; putchar('\n'), i++)
for (int j = 1; j <= n; j++)
cout << base.a[i][j] << " ";
return 0;
}
终于到了最后的 \(\color{red}\texttt{BOSS 关卡}\) 了！你们有信心吗？加油！
点击此处进入 \(\color{red}\texttt{BOSS 关卡}\) ......
点开题目 \(\color{red}\texttt{BOSS 关卡}\) 后的你 be like（梅开二度）：
这和矩阵有什么关系吗？？？
我直接一个递推！
- 对于 \(100\%\) 的数据 \(1 \leq T \leq 100\)，\(1 \leq n \leq 2 \times 10^9\)。
\(O(Tn)\) 这 \(2\times10^{11}\) 的复杂度实在无法接受。
（呜呜呜我再也不学 c艹 了）
没关系，先看看思路！
因为发现当 \(x\le3\) 时答案为 \(1\)，所以这是最基础的情况。我们可以构造一个只有一列的矩阵：
\[\begin{pmatrix} a_3 & a_2 & a_1 \end{pmatrix} = \begin{pmatrix} 1 & 1 & 1 \end{pmatrix} \]显然，这三个元素都是 \(1\)。
那么，假设我想要得到 \(a_4\)，该怎么办呢？所以，我们需要进行一种运算，让上面的矩阵变化一下，像这样：
\[\begin{pmatrix} a_3 & a_2 & a_1 \end{pmatrix} \to \begin{pmatrix} a_4 & a_3 & a_2 \end{pmatrix} \]更加通用一点：
\[\begin{pmatrix} a_x & a_{x-1} & a_{x-2} \end{pmatrix} \to \begin{pmatrix} a_{x+1} & a_{x} & a_{x-1} \end{pmatrix} \]可以发现，矩阵中的每个元素的项数都向前推进了 \(1\)。那么，我们大概可以写出伪代码：
如果 \(x\le3\)
输出 \(1\)
否则
执行运算 \(n-3\) 次（重要！）
并输出答案矩阵 \(1\) 行 \(1\) 列
特判（对于特殊情况的判断）和输出应该没什么问题，主要是为什么运算恰好要执行 \(n-3\) 次呢？稍微画个图模拟一下就好了。
还是假设要获取 \(a_4\)，则执行运算 \(4-3=1\) 次。在执行 \(1\) 次运算后，
\[\begin{pmatrix} a_3 & a_2 & a_1 \end{pmatrix} \]变为
\[\begin{pmatrix} a_4 & a_3 & a_2 \end{pmatrix} \]这样就刚好在第 \(1\) 行 \(1\) 列得到 \(a_4\) 啦！
那么，说了这么久，这个神秘的运算是什么呢？当当当当~，他就是我们的——矩阵乘法！
没错，所谓的变换，其实就是乘上了一个特殊的矩阵！那么，这个矩阵长什么样呢？让我们一起来推理吧。
（此处应配上推理の小曲）
我们可以先列一个表格，表格的行代表矩阵 \(\begin{pmatrix}a_3 & a_2 & a_1 \end{pmatrix}\) 的元素，列代表递推时与这些元素相关的元素。像这样：（表格可能在博客里渲染不出来，凑合着看吧，抱歉）
| \(a_x\) | \(a_{x-1}\) | \(a_{x-2}\) | |
|---|---|---|---|
| \(a_{x-1}\) | |||
| \(a_{x-2}\) | |||
| \(a_{x-3}\) |
好了，对于 \(a_x\)，我们该怎么填他那一列呢？我们可以观察到递推式 \(a_x=a_{x-1}+a_{x-3}\)，所以有：
\[a_x=a_{x-1}\times1+a_{x-2}\times0+a_{x-3}\times1 \]观察系数 \(1,0,1\)，把这些系数填入表格中：
| \(a_x\) | \(a_{x-1}\) | \(a_{x-2}\) | |
|---|---|---|---|
| \(a_{x-1}\) | \(1\) | ||
| \(a_{x-2}\) | \(0\) | ||
| \(a_{x-3}\) | \(1\) |
后面的也以此类推：
\[a_{x-1}=a_{x-1}\times1+a_{x-2}\times0+a_{x-3}\times0 \]\[a_{x-2}=a_{x-1}\times0+a_{x-2}\times1+a_{x-3}\times0 \]| \(a_x\) | \(a_{x-1}\) | \(a_{x-2}\) | |
|---|---|---|---|
| \(a_{x-1}\) | \(1\) | \(1\) | \(0\) |
| \(a_{x-2}\) | \(0\) | \(0\) | \(1\) |
| \(a_{x-3}\) | \(1\) | \(0\) | \(0\) |
这样，我们就可以推出这个神秘的矩阵了：
\[\begin{pmatrix} 1 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{pmatrix} \]好了，现在我们终于知道了，一次神秘操作，就是将让 \(\begin{pmatrix}a_3 & a_2 & a_1 \end{pmatrix}\) 这个矩阵乘上\( \begin{pmatrix} 1 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{pmatrix} \)。这时候就有人问了：
一次矩阵乘法的时间复杂度还没有递推快，这根本就没有优化嘛。
等等！我们把这个式子展开：
\[\begin{aligned} {} & \begin{pmatrix} 1 & 1 & 1 \end{pmatrix} \times \begin{pmatrix} 1 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{pmatrix} \times \begin{pmatrix} 1 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{pmatrix} \cdots \times \begin{pmatrix} 1 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{pmatrix} \\ = & \begin{pmatrix} 1 & 1 & 1 \end{pmatrix} \times \begin{pmatrix} 1 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{pmatrix} ^{n-3} \end{aligned} \]不是吧！这居然变成了一个矩阵快速幂？！！
也就是说，我们可以用快速幂计算 \(\begin{pmatrix} 1 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{pmatrix} ^{n-3}\)，并乘上初始矩阵 \(\begin{pmatrix} 1 & 1 & 1 \end{pmatrix}\)。这样，我们成功地把时间复杂度从 \(O(Tn)\) 优化到了 \(O(T\log n)\)！（矩阵快速幂是 \(O(\log n)\)，因为矩阵很小，矩阵乘法只计算 \(9\) 次，是一个很小的常数）
下面奉上代码：（标准的矩阵加速思想）
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MOD = 1e9 + 7;
int T, n;
struct Matrix {
ll a[5][5];
Matrix() {
memset(a, 0, sizeof a);
}
Matrix operator*(const Matrix &x)const { // 矩阵乘法
Matrix res;
for (int i = 1; i <= 3; i++)
for (int j = 1; j <= 3; j++)
for (int k = 1; k <= 3; k++)
res.a[i][j] = (res.a[i][j] % MOD + a[i][k] % MOD * x.a[k][j] % MOD) % MOD;
return res;
}
void mems() {
memset(a, 0, sizeof a);
}
}ans, base;
void init() { // 初始化两个矩阵
ans.mems(), base.mems(); // 记得清空！
ans.a[1][1] = ans.a[1][2] = ans.a[1][3] = 1;
base.a[1][1] = base.a[1][2] = base.a[2][3] = base.a[3][1] = 1;
}
void expow(int b) { // 矩阵快速幂，是在 ans 矩阵的基础上乘的
while (b) {
if (b & 1) ans = ans * base;
base = base * base, b >>= 1;
}
}
int main() {
cin >> T;
while (T --) {
cin >> n;
init(); // 初始化不能忘
if (n <= 3) { // 特判
cout << "1\n";
continue;
}
expow(n - 3); // 计算特殊矩阵的 n - 3 次方，已经乘到了 ans 里
cout << ans.a[1][1] << "\n"; // 输出答案！芜湖！
}
return 0; // 快乐结束
}
就这样，我们完成了矩阵加速递推。
再次声明矩阵快速幂（矩阵加速）时间复杂度：\(O(N^3\log n)\)，其中 \(N\) 为矩阵的行数（列数），\(n\) 为快速幂的规模 \(a^n\)。
小提示：关于 \(base\) 矩阵的构造
就是这个 \(\begin{pmatrix} 1 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{pmatrix}\) 矩阵。
可以这样：我们要推导出 \(\begin{pmatrix}a_x & a_{x-1} & a_{x-2}\end{pmatrix}\)，那么这个矩阵从哪里来？当然是从 \(\begin{pmatrix}a_{x-1} & a_{x-2} & a_{x-3}\end{pmatrix}\) 来。所以，表格才长这样：
| \(a_x\) | \(a_{x-1}\) | \(a_{x-2}\) | |
|---|---|---|---|
| \(a_{x-1}\) | |||
| \(a_{x-2}\) | |||
| \(a_{x-3}\) |
那么，能不能构造一个行列数各不相同的矩阵，而不是一个 \(n\times n\) 的矩阵呢？答案是不可以，因为我们要计算 \(\begin{pmatrix} 1 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{pmatrix}\) 这种矩阵的幂，那如果行和列不相等，相乘的两个矩阵的行列也不相等，就无法进行矩阵乘法。比如这个：
\[\begin{pmatrix} 1 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} \times \begin{pmatrix} 1 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} \]可以看到，左边 \(2\) 行，右边 \(3\) 列，显然不相等，无法进行矩阵乘法。
\[\huge\texttt{Part 6}\small\texttt{ Thank you!} \]\[\color{#FE4C6E}\texttt{你}\color{#F39C11}\texttt{居}\color{#FFC116}\texttt{然}\color{#52C41A}\texttt{看}\color{#3498DB}\texttt{完}\color{#9D3DCF}\texttt{了}\color{#0E1D69}\texttt{！} \]这篇文章花费了我很多时间，希望你喜欢！
对了，你学会了吗？是不是，矩阵也并没有那么难？
这应该是我的【精选】文章中的第一篇，没想到写的是矩阵方面的。
总之，很感谢你的阅读！希望你能从我这学到点东西！