大家好，我是Weekoder！
今天要講的內容是矩陣加速！
這時候就有人說了：
\(\tiny{\texttt{Weekoder 這麼蒻，怎麼會矩陣啊。還給我們講，真是十惡不赦！}}\)
不不不，容我解釋。在經過我的研究後，我發現基本的矩陣運算和矩陣加速都並沒有那麼難。只要繼續往下看，相信你也能學會！
注意：以下內容的學習難度將會用顏色表示，與洛谷題目難度順序一致，即 \(\color{#FE4C61}\texttt{紅}\color{#000000}<\color{#F39C11}\texttt{橙}\color{#000000}<\color{#FFC116}\texttt{黃}\color{#000000}<\color{#52C41A}\texttt{綠}\color{#000000}\)。（並不對標洛谷題目難度，只作為學習難易度參考）
矩陣和二維陣列很像，是由 \(m\times n\) 個數排列成 \(m\) 行 \(n\) 列的一張表，由於排列出來的表是一個矩形，故稱其為矩陣。矩陣長這個樣子：
可以看到，矩陣中的每個元素都有著對應的行和列，我們把一個矩陣記作 \(A\)，第 \(i\) 行 \(j\) 列的元素即為 \(a_{ij}\)。更形式化的，寫作：
其中 \(\mathbb{F}\) 為數域，一般取為實數域 \(\mathbb{R}\) 或複數域 \(\mathbb{C}\)。（看不懂沒事，蒟蒻自行走開QWQ）
\(\texttt{1.零矩陣}\)
元素全部為 \(0\) 的矩陣稱為零矩陣。像這樣：
零矩陣記作 \(0_{m\times n}\)，就是在 \(0\) 下面加上矩陣的大小 \(m\times n\)。你可以把零矩陣看做數字 \(0\)，任何數乘以 \(0\) 都得 \(0\)。
\(\texttt{2.對角矩陣}\)
只有主對角線上的元素有值，其餘元素為 \(0\) 的矩陣稱為對角矩陣。
注：主對角線為矩陣中從左上角到右下角的一條對角線。
對角矩陣根據主對角線上的值，記作 \(\text{diag(}\text{a}_1,\text{a}_2,\ldots,\text{a}_n\text{)}\)。
\(\texttt{3.單位矩陣}\)
主對角線上的元素均為 \(1\)，其餘元素為 \(0\) 的矩陣稱為單位矩陣。
單位矩陣記作 \(I\)。
記得分數中的概念分數單位嗎？矩陣單位和分數單位的“地位”差不多，代表的都是最基礎的，最小的獨立個體。你可以把單位矩陣看做數字 \(1\)，任何數乘以 \(1\) 都等於它本身。
最基礎的，常見的特殊矩陣就是這些了。當然，還有很多的特殊矩陣，不過我們暫時用不到。
\(\texttt{1.相等}\)
若對於矩陣 \(A,B\)，所有的 \(i,j\) 都有 \(a_{ij}=b_{ij}\) 且矩陣的行和列相等，則稱矩陣 \(A,B\) 相等。
其實就是兩個矩陣長得一模一樣。
\(\texttt{2.矩陣加（減）法}\)
若要求 \(A,B\) 兩個矩陣之和，即 \(C=A+B\)，則對於任意 \(i,j\)，滿足 \(c_{ij}=a_{ij}+b_{ij}\)。要求矩陣行列相等。
總結一句話：對應位置相加。
矩陣加法滿足交換律和結合律：
減法同理，對應位置相減。
\(\texttt{3.矩陣數乘}\)
數 \(\lambda\)（一個數字） 乘以矩陣 \(A\)，記作 \(\lambda A\)，即為矩陣數乘運算。若有 \(B=\lambda A\)，則對於任意 \(i,j\) 都滿足 \(b_{ij}=\lambda a_{ij}\)。
還是一句話：對應位置相乘。
雖然矩陣乘法也屬於矩陣運算，但難度比前面的都高，而且是今天的重點內容，所以單獨放出來講，故記為 \(\texttt{Part 3.5}\)。（話說你們沒有發現難度變成黃了嗎）
上例題！（雖然難度是橙）
先看矩陣乘法的定義：若有 \(n\) 行 \(m\) 列矩陣 \(A\) 和 \(m\) 行 \(k\) 列的矩陣 \(B\)（\(A\) 的行與 \(B\) 的列相等），則 \(n\) 行 \(k\) 列的矩陣 \(C=A\times B\) 滿足
只要列舉 \(i,j\)（範圍是 \(n,k\)），並套用公式就能用 \(O(n^3)\) 的時間複雜度解決這個問題。
我知道，這看起來根本不是新手蒟蒻能看懂的。那我就用人話來講講矩陣乘法。
矩陣乘法並不是一個一個乘，而是行對應列乘。怎麼個乘法呢？我們來看看下面兩個矩陣相乘的例子。
第一個矩陣為 \(A\)，第二個矩陣為 \(B\)。
我們先取出 \(A\) 的第一行。像這樣：
再取出 \(B\) 的第一列。像這樣：
不對，你給我轉過來。
現在終於可以相乘了。逐位相乘得出結果：
得出了結果 \(\begin{pmatrix} 10 & 0 & 6 \end{pmatrix}\)。再將每一位相加：
還記得我們之前是怎麼取的嗎？我們取了 \(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\)
矩陣乘法不滿足交換律。（這是重點！）
有了矩陣乘法，我們還可以結合上面的特殊矩陣得到一些性質：
快到今天的主題了！上例題！
點開題目後的你 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\)，所以這是最基礎的情況。我們可以構造一個只有一列的矩陣：
顯然，這三個元素都是 \(1\)。
那麼，假設我想要得到 \(a_4\)，該怎麼辦呢？所以，我們需要進行一種運算，讓上面的矩陣變化一下，像這樣：
更加通用一點：
可以發現，矩陣中的每個元素的項數都向前推進了 \(1\)。那麼，我們大概可以寫出虛擬碼：
如果 \(x\le3\)
輸出 \(1\)
否則
執行運算 \(n-3\) 次（重要！）
並輸出答案矩陣 \(1\) 行 \(1\) 列
特判（對於特殊情況的判斷）和輸出應該沒什麼問題，主要是為什麼運算恰好要執行 \(n-3\) 次呢？稍微畫個圖模擬一下就好了。
還是假設要獲取 \(a_4\)，則執行運算 \(4-3=1\) 次。在執行 \(1\) 次運算後，
變為
這樣就剛好在第 \(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}\)，所以有：
觀察係數 \(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\) | \(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}a_3 & a_2 & a_1 \end{pmatrix}\) 這個矩陣乘上\( \begin{pmatrix} 1 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{pmatrix} \)。這時候就有人問了：
一次矩陣乘法的時間複雜度還沒有遞推快，這根本就沒有最佳化嘛。
等等！我們把這個式子展開：
不是吧！這居然變成了一個矩陣快速冪？！！
也就是說，我們可以用快速冪計算 \(\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}\) 這種矩陣的冪，那如果行和列不相等，相乘的兩個矩陣的行列也不相等，就無法進行矩陣乘法。比如這個：
可以看到，左邊 \(2\) 行，右邊 \(3\) 列，顯然不相等，無法進行矩陣乘法。
這篇文章花費了我很多時間，希望你喜歡！
對了，你學會了嗎？是不是，矩陣也並沒有那麼難？
這應該是我的【精選】文章中的第一篇，沒想到寫的是矩陣方面的。
總之，很感謝你的閱讀！希望你能從我這學到點東西！