第 8.4 章：矩阵快速幂

8.4.0 从一道题出发：数到宇宙尽头

🔢 题目：Fibonacci，但 N 很大 求斐波那契数列第 N 项 mod (10⁹+7)。定义 f(0) = 0，f(1) = 1，f(n) = f(n−1) + f(n−2)。
数据规模： 0 ≤ N ≤ 10¹⁸。时间限制：1 秒。

先试试你会的方法

思路 1：递推。 for 循环从 2 推到 N。N = 10¹⁸ 时，即使每秒做 10⁹ 次加法，也需要约 32 年。淘汰。

思路 2：记忆化递归 / 数学公式（通项公式）。 前者同样是 O(N)；后者（Binet 公式）依赖 √5 的浮点精度，N 稍大就在模意义下彻底失真——模运算世界里没有 √5（本章范围内）。

思路 3：矩阵。 把「递推一步」写成一个 2×2 矩阵的乘法，递推 N 步 = 矩阵自乘 N−1 次。矩阵幂……不就是幂吗？——用第 8.1 章的快速幂思想，O(log N) 次矩阵乘法搞定。10¹⁸ ≈ 2⁶⁰，只需约 60 轮，毫秒级出答案。

本章路线：

学习目标：

• 实现带取模的矩阵乘法
• 用快速幂计算矩阵的 n 次幂，O(k³ log n)
• 把线性递推转化为矩阵形式（以 Fibonacci 为例）
• 判断一道题的递推能否用矩阵加速

8.4.1 矩阵乘法：定义与实现

定义

矩阵 A 是 n×m 的表，B 是 m×s 的表。乘积 C = A×B 是 n×s 的表，其中：

只有当 A 的列数等于 B 的行数时才能相乘。以 2×2 为例：

⎡ 1  2 ⎤   ⎡ 3  1 ⎤   ⎡ 1×3+2×5   1×1+2×4 ⎤   ⎡ 13  9 ⎤
⎣ 4  5 ⎦ × ⎣ 5  4 ⎦ = ⎣ 4×3+5×5   4×1+5×4 ⎦ = ⎣ 37 24 ⎦

实现

// 解答：矩阵乘法（带取模）
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef vector<vector<ll>> Matrix;

const ll MOD = 1e9 + 7;

// A: n×m, B: m×s → C: n×s
Matrix mat_mul(const Matrix& A, const Matrix& B) {
    int n = A.size(), m = B.size(), s = B[0].size();
    Matrix C(n, vector<ll>(s, 0));
    for (int i = 0; i < n; i++)
        for (int k = 0; k < m; k++)
            for (int j = 0; j < s; j++)
                C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % MOD;
    return C;
}

int main() {
    Matrix A = {{1, 2}, {4, 5}};
    Matrix B = {{3, 1}, {5, 4}};
    Matrix C = mat_mul(A, B);
    // 输出 13 9 / 37 24
    for (auto& row : C) {
        for (ll v : row) cout << v << " ";
        cout << "\n";
    }
    return 0;
}

三个细节：

1. 循环顺序 i-k-j（而非 i-j-k）：k 在中 layer，B[k][j] 沿 j 连续访问，缓存友好，常数更小。
2. A[i][k] * B[k][j] 必须是 ll * ll——两个模后的数（< 10⁹）相乘约 10¹⁸，int 直接爆炸。
3. 每次累加后立即 % MOD。每项乘积 < 10¹⁸，而 C[i][j] 始终 < MOD，所以「加上一项再取模」的中间值 < 10¹⁸ + MOD，远在 long long 范围内；若累加多项才取模就会溢出。

注意：矩阵乘法不满足交换律（A×B ≠ B×A 一般成立），但满足结合律：(A×B)×C = A×(B×C)。结合律正是快速幂的入场券。

8.4.2 矩阵快速幂：O(k³ log n)

单位矩阵：矩阵世界的「1」

数的快速幂中 result 初始化为 1。矩阵世界里，「乘了等于没乘」的东西是单位矩阵 I——对角线为 1、其余为 0：

⎡ 1 0 ⎤
⎣ 0 1 ⎦

任意矩阵 M 满足 M × I = I × M = M。

算法：与数的快速幂一模一样

因为矩阵乘法满足结合律，M^n = M^(n 的二进制拆分) 的推导原封不动照搬：

// 解答：矩阵快速幂 — 计算 M^n（n 可达 10^18）
Matrix mat_pow(Matrix M, ll n) {
    int sz = M.size();
    Matrix result(sz, vector<ll>(sz, 0));
    for (int i = 0; i < sz; i++) result[i][i] = 1;   // 初始化为单位矩阵！

    while (n > 0) {
        if (n & 1)
            result = mat_mul(result, M);
        M = mat_mul(M, M);   // 底数平方：M^1, M^2, M^4, M^8, ...
        n >>= 1;
    }
    return result;
}

复杂度： 每次 mat_mul 是 O(k³)（k×k 矩阵），循环 O(log n) 轮，总计 O(k³ log n)。k = 2 时约 60×16 = 不足一千次乘法，k = 100 时约 60×10⁶ = 六千万次，仍在时限内。

8.4.3 完整实战：Fibonacci 的矩阵形态

推导转移矩阵

目标：把「从 (n−1, n−2) 推进到 (n, n−1)」写成一次矩阵乘法。把相邻两项摞成状态向量，找一个矩阵 M 使得：

⎛ f(n)   ⎞   ⎛ ?  ? ⎞ ⎛ f(n-1) ⎞
⎝ f(n-1) ⎠ = ⎝ ?  ? ⎠ ⎝ f(n-2) ⎠

逐项凑系数：f(n) = 1·f(n−1) + 1·f(n−2)，f(n−1) = 1·f(n−1) + 0·f(n−2)。于是：

$$\begin{pmatrix} f(n) \ f(n-1) \end{pmatrix} = \begin{pmatrix} 1 & 1 \ 1 & 0 \end{pmatrix} \times \begin{pmatrix} f(n-1) \ f(n-2) \end{pmatrix}$$

递推 n−1 步，就是 M 自乘 n−1 次：

$$\begin{pmatrix} f(n) \ f(n-1) \end{pmatrix} = \begin{pmatrix} 1 & 1 \ 1 & 0 \end{pmatrix}^{n-1} \times \begin{pmatrix} f(1) \ f(0) \end{pmatrix}$$

完整代码

// 解答：Fibonacci 第 N 项 mod 1e9+7 — O(k^3 log n)
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef vector<vector<ll>> Matrix;

const ll MOD = 1e9 + 7;

Matrix mat_mul(const Matrix& A, const Matrix& B) {
    int n = A.size(), m = B.size(), s = B[0].size();
    Matrix C(n, vector<ll>(s, 0));
    for (int i = 0; i < n; i++)
        for (int k = 0; k < m; k++)
            for (int j = 0; j < s; j++)
                C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % MOD;
    return C;
}

Matrix mat_pow(Matrix M, ll n) {
    int sz = M.size();
    Matrix result(sz, vector<ll>(sz, 0));
    for (int i = 0; i < sz; i++) result[i][i] = 1;
    while (n > 0) {
        if (n & 1) result = mat_mul(result, M);
        M = mat_mul(M, M);
        n >>= 1;
    }
    return result;
}

ll fibonacci(ll n) {
    if (n <= 1) return n;
    Matrix trans = {{1, 1}, {1, 0}};
    Matrix R = mat_pow(trans, n - 1);
    // [f(n), f(n-1)]^T = R × [f(1), f(0)]^T = R × [1, 0]^T
    // 第一个分量 = R[0][0]×1 + R[0][1]×0 = R[0][0]
    return R[0][0];
}

int main() {
    for (int i = 0; i <= 10; i++)
        cout << "f(" << i << ") = " << fibonacci(i) << "\n";
    cout << "f(10^18) = " << fibonacci(1000000000000000000LL) << "\n";
    return 0;
}

手算验证：为什么 f(10) = 55

M² = [[1,1],[1,0]] × [[1,1],[1,0]]：

M² = ⎡ 1×1+1×1  1×1+1×0 ⎤   ⎡ 2 1 ⎤
     ⎣ 1×1+0×1  1×1+0×0 ⎦ = ⎣ 1 1 ⎦

M⁴ = M² × M² = ⎡ 2×2+1×1  2×1+1×1 ⎤   ⎡ 5 3 ⎤
               ⎣ 1×2+1×1  1×1+1×1 ⎦ = ⎣ 3 2 ⎦

M⁸ = M⁴ × M⁴ = ⎡ 5×5+3×3  5×3+3×2 ⎤   ⎡ 34 21 ⎤
               ⎣ 3×5+2×3  3×3+2×2 ⎦ = ⎣ 21 13 ⎦

f(10) 需要 M⁹ = M⁸ × M：

M⁹ = M⁸ × M = ⎡ 34×1+21×1  34×1+21×0 ⎤   ⎡ 55 34 ⎤
              ⎣ 21×1+13×1  21×1+13×0 ⎦ = ⎣ 34 21 ⎦

f(10) = M⁹[0][0] × f(1) = 55 × 1 = 55 ✓。事实上能看出 Mⁿ 的规律：左上角恰好是 f(n+1)，右上角是 f(n)——矩阵每一格都在讲 Fibonacci 的故事。

性能对比： 递推 O(N) 对 N = 10¹⁸ 要算 32 年；矩阵快速幂约 60 轮 × 2 次矩阵乘法 × 8 次标量乘法 ≈ 1000 次乘法，不到一微秒。

8.4.4 适用判断：什么递推能矩阵加速？

判据：线性齐次递推（可带常数项）

一个 k 阶递推能写成一阶矩阵形式，当且仅当新状态是旧状态的线性组合（各项次数为 1，可以乘系数、加常数）：

✓ f(n) = a₁·f(n-1) + a₂·f(n-2) + ... + aₖ·f(n-k)          — 纯线性
✓ f(n) = 2·f(n-1) + 3·f(n-2) + 7                          — 常数项：加一行恒为 1 的状态
✓ f(n) = f(n-1) + n²                                       — n 的多项式：把 n², n, 1 也放进状态向量
✗ f(n) = f(n-1) × f(n-2)                                   — 乘法非线性，矩阵无能为力
✗ f(n) = f(n-1)² + 1                                       — 平方非线性

常数项的处理： f(n) = 2f(n−1) + 7，状态向量扩成 [f(n), 1]：

⎛ f(n) ⎞   ⎛ 2 7 ⎞ ⎛ f(n-1) ⎞
⎝  1   ⎠ = ⎝ 0 1 ⎠ ⎝   1    ⎠

含 n 的多项式项： f(n) = f(n−1) + n，状态向量扩成 [f(n), n, 1]，转移矩阵里放好 n ← n+1 的递进。状态每扩一维，k 加 1，复杂度 O(k³ log n) 里 k³ 增长——扩到 k ≈ 100 依然是毫秒级。

典型信号

• 「求第 N 项，N ≤ 10¹⁸，答案取模」——N 大到只能 log 级算法，递推又是线性的 → 矩阵快速幂
• 「走 N 步后的方案数」且转移固定 → 转移矩阵的 N 次幂（图上 N 步路径计数，k = 顶点数）
• 数据结构题里「重复执行 N 次某个线性变换」→ 同理

与 DP 的关系

若 N ≤ 10⁷ 左右，直接 O(N) 递推更简单——矩阵快速幂只在 N 巨大 + 递推线性时才有出场必要。这是典型的 Platinum 套路：题目数据范围就是算法提示。

⚠️ 常见错误

1. result 没初始化为单位矩阵。 初始化成全 0 矩阵的话，result × M 永远是全 0——快速幂结果直接报废。矩阵的「1」是 I，不是全 1，也不是全 0。
2. 乘法顺序写反。 矩阵乘法不满足交换律！result = mat_mul(M, result) 与 mat_mul(result, M) 一般结果不同。按模板顺序写，别「顺手」交换。
3. k³ 里的溢出。 A[i][k] * B[k][j] 两个模后（< 10⁹）的数相乘约 10¹⁸，必须全程 long long。若模数本身 ≥ 2×10⁹（罕见），还要上 __int128。
4. 状态向量维度不够。 常数项忘了加一维、n 的多项式项忘了带上，转移矩阵就「装不下」递推。写出递推后先数清楚有几个自由变量。
5. 忘记矩阵乘法不满足交换律却用了快速幂。 快速幂只需要结合律（成立），但推理「M^(a+b) = M^a × M^b」时顺序不能乱：M^a × M^b = M^(a+b) ✓（同底数才可合并）。
6. n = 0 或 n = 1 的边界。 Fibonacci 里 if (n <= 1) return n; 必须在矩阵计算之前——mat_pow(M, -1) 是未定义行为。

📋 章节小结

前置/后置： 快速幂的数的版本见第 8.1 章：数论基础；图上 N 步路径计数（转移矩阵幂的另一种读法）可对照第 5 章 Floyd–Warshall 的矩阵视角。

🏋️ 练习题

🟢 简单

8.4-E1. 矩阵幂验证
给定 2×2 矩阵 M 和整数 n，输出 M^n 的各元素 mod 10⁹+7。
数据规模： 0 ≤ n ≤ 10¹⁸，M 的元素 ∈ [0, 10⁹]。

8.4-E2. Tribonacci
求满足 t(0)=0, t(1)=1, t(2)=1, t(n)=t(n−1)+t(n−2)+t(n−3) 的数列第 N 项 mod 10⁹+7。
数据规模： 0 ≤ N ≤ 10¹⁸。

🟡 中等

8.4-M1. 爬楼梯（带常数项）
每次可以上 1 阶或 2 阶，但第 0 阶到第 1 阶之间有一段电梯直达（相当于从 f(0) 走到 f(1) 有 2 种方式）。递推为 f(n) = 2·f(n−1) + f(n−2)，f(0)=1, f(1)=2。求第 N 阶的方案数 mod 10⁹+7。
数据规模： 0 ≤ N ≤ 10¹⁸。

🔴 困难

8.4-H1. 图上的 N 步路径计数 (USACO Platinum 风格)
给定有向图（K 个顶点，E 条边）和起点 s、终点 t，求从 s 走恰好 N 步到 t 的路径数 mod 10⁹+7（每步沿一条有向边走）。
数据规模： 1 ≤ K ≤ 100，0 ≤ E ≤ K²，1 ≤ N ≤ 10¹⁸。

第 8.4 章结束 — 另请参阅：数论基础 | 数学部分总览