矩阵快速幂与线性递推

前言

02-04 快速幂与逆元 里已经介绍过矩阵快速幂的基本原理：把递推式 $a_n = 2a_{n-1} + a_{n-2}$ 写成矩阵形式，然后用快速幂在 $O(\log N)$ 内求出第 $N$ 项。那篇的例题 5 给出了一个 $2 \times 2$ 矩阵的完整实现。

这篇是矩阵快速幂的进阶篇。我们会处理更高阶的递推（Tribonacci 三阶递推、$k$ 阶递推），以及更复杂的转移矩阵构造方法。核心思想不变：任何常系数线性递推都可以写成矩阵乘法。

问题的本质

从递推到矩阵

斐波那契递推 $f_n = f_{n-1} + f_{n-2}$ 可以写成：

$\begin{pmatrix} f_n \\ f_{n-1} \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} f_{n-1} \\ f_{n-2} \end{pmatrix}$

展开验证：$f_n = 1 \cdot f_{n-1} + 1 \cdot f_{n-2}$，$f_{n-1} = 1 \cdot f_{n-1} + 0 \cdot f_{n-2}$。✓

递推 $n-1$ 次就得到：

$\begin{pmatrix} f_n \\ f_{n-1} \end{pmatrix} = M^{n-2} \begin{pmatrix} f_2 \\ f_1 \end{pmatrix}$

其中 $M = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}$。$M^{n-2}$ 用快速幂 $O(\log N)$ 算出。

高阶递推的推广

三阶递推 $a_n = a_{n-1} + a_{n-2} + a_{n-3}$（Tribonacci）的转移矩阵是 $3 \times 3$：

$\begin{pmatrix} a_n \\ a_{n-1} \\ a_{n-2} \end{pmatrix} = \begin{pmatrix} 1 & 1 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{pmatrix} \begin{pmatrix} a_{n-1} \\ a_{n-2} \\ a_{n-3} \end{pmatrix}$

规律：转移矩阵的第一行就是递推系数（$1, 1, 1$），其余行是单位向量偏移（把第 $i$ 行第 $i-1$ 列设为 1，用来”向下移动”状态）。

一般地，$k$ 阶常系数线性递推 $a_n = c_1 a_{n-1} + c_2 a_{n-2} + \cdots + c_k a_{n-k}$ 对应 $k \times k$ 转移矩阵：

$M = \begin{pmatrix} c_1 & c_2 & \cdots & c_k \\ 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & & \ddots & \vdots \\ 0 & 0 & \cdots & 1 & 0 \end{pmatrix}$

复杂度

$K \times K$ 矩阵乘法 $O(K^3)$，快速幂 $O(\log N)$ 次乘法，总时间 $O(K^3 \log N)$。$K=3, N=10^{18}$ 约 $27 \times 60 \approx 1600$ 次乘法，瞬间跑完。

理论 + 代码

通用矩阵快速幂模板

#include <cstdio>
#include <cstring>
using namespace std;

const long long MOD = 1000000007;

struct Mat {
    static const int MAXK = 10;
    long long a[MAXK][MAXK];
    int n; // 矩阵大小 n × n
    Mat(int n = 2) : n(n) { memset(a, 0, sizeof(a)); }
};

Mat mul(const Mat& A, const Mat& B) {
    Mat C(A.n);
    for (int i = 0; i < A.n; i++)
        for (int k = 0; k < A.n; k++)
            for (int j = 0; j < A.n; j++)
                C.a[i][j] = (C.a[i][j] + A.a[i][k] * B.a[k][j]) % MOD; // ① 三重循环
    return C;
}

Mat power(Mat A, long long e) {
    Mat R(A.n);
    for (int i = 0; i < A.n; i++) R.a[i][i] = 1;  // ② 单位矩阵
    while (e > 0) {
        if (e & 1) R = mul(R, A);     // ③ 累乘
        A = mul(A, A);                 // ④ 自乘
        e >>= 1;
    }
    return R;
}

逐行解析：

• ① 矩阵乘法三重循环。注意 k 放中间可以减少 cache miss。每次乘法后取模。
• ② 单位矩阵 $I$：$I_{ii} = 1$，其余 0。$I \times A = A$，是矩阵乘法的”1”。
• ③④ 和普通快速幂结构完全一样——③ 累乘到结果，④ 底数自乘。

斐波那契的转移矩阵

$K=2$ 的特例，在 02-04 已经实现过。这里用通用模板：

// 斐波那契: f[n] = f[n-1] + f[n-2]
Mat M(2);
M.a[0][0] = 1; M.a[0][1] = 1;  // 第一行 = [1, 1]
M.a[1][0] = 1; M.a[1][1] = 0;  // 第二行 = [1, 0]

Mat R = power(M, N - 2);        // M^{N-2}
// f[N] = R.a[0][0] * f[2] + R.a[0][1] * f[1]
long long ans = (R.a[0][0] + R.a[0][1]) % MOD; // f[2]=1, f[1]=1

Tribonacci 的转移矩阵

$K=3$，递推 $a_n = a_{n-1} + a_{n-2} + a_{n-3}$：

Mat M(3);
M.a[0][0] = 1; M.a[0][1] = 1; M.a[0][2] = 1; // 第一行 = [1, 1, 1]
M.a[1][0] = 1; M.a[1][1] = 0; M.a[1][2] = 0; // 第二行 = [1, 0, 0]
M.a[2][0] = 0; M.a[2][1] = 1; M.a[2][2] = 0; // 第三行 = [0, 1, 0]

构造规则：第一行是递推系数 $(c_1, c_2, c_3)$。第 $i$ 行 ($i \ge 2$) 的第 $i-1$ 列为 1，其余为 0——这行的效果是把 $a_{n-i+1}$ “传递”给下一轮的状态向量。

例题

例题 1：M&A 054 — Fibonacci Hard（mod $10^9$）

题目：斐波那契数列 $f_1 = 1, f_2 = 1, f_n = f_{n-1} + f_{n-2}$。求 $f_N \bmod 10^9$。

数据范围：$3 \le N \le 10^{18}$

—— AtCoder M&A 054

分析：$2 \times 2$ 矩阵快速幂。注意模数是 $10^9$ 不是 $10^9+7$。

输入样例：10，输出：55

代码：

#include <cstdio>
#include <cstring>
using namespace std;

const long long MOD = 1000000000; // 注意：10^9 不是 10^9+7

struct Mat {
    long long a[2][2];
};

Mat mul(Mat A, Mat B) {
    Mat C; memset(C.a, 0, sizeof(C.a));
    for (int i = 0; i < 2; i++)
        for (int k = 0; k < 2; k++)
            for (int j = 0; j < 2; j++)
                C.a[i][j] = (C.a[i][j] + A.a[i][k] * B.a[k][j]) % MOD;
    return C;
}

Mat power(Mat A, long long e) {
    Mat R; memset(R.a, 0, sizeof(R.a));
    R.a[0][0] = R.a[1][1] = 1;                        // ① 单位矩阵
    while (e > 0) {
        if (e & 1) R = mul(R, A);
        A = mul(A, A);
        e >>= 1;
    }
    return R;
}

int main() {
    long long N;
    scanf("%lld", &N);

    Mat M; memset(M.a, 0, sizeof(M.a));
    M.a[0][0] = 1; M.a[0][1] = 1;                     // ② 转移矩阵
    M.a[1][0] = 1; M.a[1][1] = 0;

    Mat R = power(M, N - 2);                           // ③ M^{N-2}
    long long ans = (R.a[0][0] + R.a[0][1]) % MOD;    // ④ f[2]*M^{N-2}[0][0] + f[1]*M^{N-2}[0][1]
    printf("%lld\n", ans);
}

逐行解析：

• ① $I = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}$。快速幂的初始值。
• ② $M = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}$。
• ③ 从 $f_2$ 推到 $f_N$ 需要乘 $N-2$ 次。
• ④ $f_N = M^{N-2}[0][0] \cdot f_2 + M^{N-2}[0][1] \cdot f_1 = R[0][0] \times 1 + R[0][1] \times 1$。

验证：$N=10$，$M^8 = \begin{pmatrix} 34 & 21 \\ 21 & 13 \end{pmatrix}$。$f_{10} = 34 + 21 = 55$。✓

---

例题 2：M&A 056 — Recurrence Formula 2（Tribonacci）

题目：$a_1 = 1, a_2 = 1, a_3 = 2$，$a_n = a_{n-1} + a_{n-2} + a_{n-3}$。求 $a_N \bmod 10^9+7$。

数据范围：$4 \le N \le 10^{18}$

—— AtCoder M&A 056

分析：$3 \times 3$ 矩阵快速幂。递推系数 $(1, 1, 1)$。

输入样例：10，输出：149

代码：

#include <cstdio>
#include <cstring>
using namespace std;

const long long MOD = 1000000007;

struct Mat {
    long long a[3][3];
};

Mat mul(Mat A, Mat B) {
    Mat C; memset(C.a, 0, sizeof(C.a));
    for (int i = 0; i < 3; i++)
        for (int k = 0; k < 3; k++)
            for (int j = 0; j < 3; j++)
                C.a[i][j] = (C.a[i][j] + A.a[i][k] * B.a[k][j]) % MOD;
    return C;
}

Mat power(Mat A, long long e) {
    Mat R; memset(R.a, 0, sizeof(R.a));
    R.a[0][0] = R.a[1][1] = R.a[2][2] = 1;            // ① 3×3 单位矩阵
    while (e > 0) {
        if (e & 1) R = mul(R, A);
        A = mul(A, A);
        e >>= 1;
    }
    return R;
}

int main() {
    long long N;
    scanf("%lld", &N);

    // ② 转移矩阵
    Mat M; memset(M.a, 0, sizeof(M.a));
    M.a[0][0] = 1; M.a[0][1] = 1; M.a[0][2] = 1; // a_n = a_{n-1} + a_{n-2} + a_{n-3}
    M.a[1][0] = 1; M.a[1][1] = 0; M.a[1][2] = 0; // a_{n-1} = a_{n-1}
    M.a[2][0] = 0; M.a[2][1] = 1; M.a[2][2] = 0; // a_{n-2} = a_{n-2}

    Mat R = power(M, N - 3);                           // ③ M^{N-3}

    // ④ a_N = R[0][0]*a_3 + R[0][1]*a_2 + R[0][2]*a_1
    long long ans = (R.a[0][0] * 2 + R.a[0][1] * 1 + R.a[0][2] * 1) % MOD;
    printf("%lld\n", ans);
}

逐行解析：

• ② 转移矩阵 $3 \times 3$。第一行 $(1,1,1)$ 是递推系数。
• ③ 从 $a_3$ 推到 $a_N$ 需要乘 $N-3$ 次。
• ④ 初始状态向量 $(a_3, a_2, a_1) = (2, 1, 1)$。$a_N = R[0][0] \times 2 + R[0][1] \times 1 + R[0][2] \times 1$。

验证：$a = (1, 1, 2, 4, 7, 13, 24, 44, 81, 149, \ldots)$。$a_{10} = 149$。✓

---

例题 3：M&A 055 — Recurrence Formula 1（含系数递推）

题目：$a_1 = 1, a_2 = 1$，$a_n = 2a_{n-1} + a_{n-2}$。求 $a_N \bmod 10^9+7$。

数据范围：$3 \le N \le 10^{18}$

—— AtCoder M&A 055

分析：递推系数 $(2, 1)$——不是 $(1, 1)$ 了。转移矩阵变成 $\begin{pmatrix} 2 & 1 \\ 1 & 0 \end{pmatrix}$。这个题在 02-04 例题 5 已经完整实现过。

代码（与 02-04 例题 5 相同）：

#include <cstdio>
using namespace std;
const long long MOD = 1000000007;

struct Mat {
    long long a[2][2];
    Mat() { a[0][0]=a[0][1]=a[1][0]=a[1][1]=0; }
};

Mat mul(Mat A, Mat B) {
    Mat C;
    for (int i = 0; i < 2; i++)
        for (int k = 0; k < 2; k++)
            for (int j = 0; j < 2; j++)
                C.a[i][j] = (C.a[i][j] + A.a[i][k] * B.a[k][j]) % MOD;
    return C;
}

Mat power(Mat A, long long e) {
    Mat ans;
    ans.a[0][0] = ans.a[1][1] = 1;                     // ① 单位矩阵
    while (e > 0) {
        if (e & 1) ans = mul(ans, A);
        A = mul(A, A);
        e >>= 1;
    }
    return ans;
}

int main() {
    long long N;
    scanf("%lld", &N);
    Mat M;
    M.a[0][0] = 2; M.a[0][1] = 1;                     // ② [2, 1] 不是 [1, 1]
    M.a[1][0] = 1; M.a[1][1] = 0;
    Mat R = power(M, N - 2);                            // ③ M^{N-2}
    long long ans = (R.a[0][0] + R.a[0][1]) % MOD;     // ④ a_2=1, a_1=1
    printf("%lld\n", ans);
}

逐行解析：

• ② 关键区别：$a_n = 2a_{n-1} + 1 \cdot a_{n-2}$，所以转移矩阵第一行是 $(2, 1)$ 而非 $(1, 1)$。
• ③④ 和斐波那契完全一样的流程——只要改转移矩阵。

验证：$a = (1, 1, 3, 7, 17, 41, 99, 239, 577, 1393)$。$a_{10} = 1393$。✓

---

例题 4：M&A 049 — Fibonacci Easy（mod $10^9+7$，$O(N)$ 递推）

题目：斐波那契数列 $f_1 = 1, f_2 = 1$，求 $f_N \bmod 10^9+7$。

数据范围：$1 \le N \le 10^6$

—— AtCoder M&A 049

分析：$N \le 10^6$，$O(N)$ 递推就够了，不需要矩阵快速幂。教材 5.3 B28 Fibonacci Easy 就是这么做的。但如果你已经写了矩阵快速幂模板，它也能轻松处理——两种方法都正确，选择取决于数据范围。

$O(N)$ 代码：

#include <cstdio>
using namespace std;
const int MOD = 1000000007;

int main() {
    int N;
    scanf("%d", &N);
    if (N <= 2) { printf("1\n"); return 0; }
    int a = 1, b = 1;                                  // ① f_1, f_2
    for (int i = 3; i <= N; i++) {
        int c = (a + b) % MOD;                         // ② f_i = f_{i-1} + f_{i-2}
        a = b;                                          // ③ 滚动：f_{i-2} ← f_{i-1}
        b = c;                                          // ④ f_{i-1} ← f_i
    }
    printf("%d\n", b);
}

逐行解析：

• ① 初始值 $f_1 = 1, f_2 = 1$。
• ②③④ 滚动数组——只保留前两项，空间 $O(1)$。
• 复杂度 $O(N)$。$N = 10^6$ 时约 $10^6$ 次加法，毫秒级。

什么时候用矩阵快速幂？ 当 $N \ge 10^9$ 或递推阶数 $K$ 较小时。$N \le 10^6$ 时 $O(N)$ 递推更简单更快。

参考文献

教材讲解 — 競技プログラミングの鉄則 第 5 章

• 5.3 B28 Fibonacci Easy（$O(N)$ 递推 解说）
• 5.4 B29 Power Hard（快速幂 解说）

基础练习 — アルゴリズムと数学 演習問題集

• 049 Fibonacci Easy（$O(N)$ 递推 mod $10^9+7$）【例题】
• 054 Fibonacci Hard（矩阵快速幂 mod $10^9$）【例题】
• 055 Recurrence Formula 1（含系数 $2 \times 2$ 矩阵快速幂）【例题】
• 056 Recurrence Formula 2（Tribonacci $3 \times 3$ 矩阵快速幂）【例题】
• 062 Teleporter（倍增/矩阵快速幂 $K \le 10^{18}$）

---

系列索引

第零章 基础工具

• 00-01 复杂度与暴力枚举
• 00-02 排序
• 00-03 累积和与差分

第一章 搜索技术

• 01-01 二分搜索
• 01-02 双指针与尺取法

第二章 数学基础

• 02-01 整数与整除
• 02-02 素数与筛法
• 02-03 同余与模运算
• 02-04 快速幂与逆元
• 02-05 组合计数
• 02-06 容斥原理
• 02-07 欧拉函数
• 02-08 博弈论基础

第三章 数据结构

• 03-01 栈、队列与单调栈
• 03-02 堆与优先队列
• 03-03 并查集
• 03-04 树状数组
• 03-05 线段树
• 03-06 懒标记线段树
• 03-07 Sparse Table 与倍增
• 03-08 字符串哈希

第四章 图论

• 04-01 图的遍历
• 04-02 最短路—Dijkstra 与 01-BFS
• 04-03 最短路—Bellman-Ford 与 Floyd
• 04-04 拓扑排序
• 04-05 最小生成树
• 04-06 强连通分量与 2-SAT
• 04-07 二分图与网络流
• 04-08 树上问题

第五章 动态规划

• 05-01 DP入门—状态与转移
• 05-02 背包问题族
• 05-03 LIS、LCS与编辑距离
• 05-04 区间DP
• 05-05 状态压缩DP
• 05-06 树形DP与数位DP
• 05-07 矩阵快速幂与线性递推

第六章 贪心

• 06-01 贪心策略与交换论证
• 06-02 区间调度与优化

第七章 字符串

• 07-01 字符串哈希详述
• 07-02 后缀数组

第八章 进阶

• 08-01 卷积与FFT
• 08-02 期望与概率
• 08-03 Zeta / Möbius 变换
• 08-04 计算几何基础
• 08-05 采样与随机化