先说我为什么想写这个。最近在牛客和力扣上连续看到好几道题斐波那契数列把 n 给到 1e18或者让你求“图上走 K 步回到原点的方案数”。这种题如果直接用循环去推一秒都跑不完但只要换成 C 的矩阵快速幂就能把 O(n) 压到 O(k³ log n)k 通常是个不超过 10 的小数字。这篇文章就是我这些年刷题、写竞赛代码、做算法课设时关于矩阵快速幂的一份个人总结适合刚学完递归和循环、准备进阶算法或者正在准备面试手撕代码的 C 开发者。我不打算写一本教科书而是尽量把我踩过的坑、验证过的小技巧都摆出来让你看完能直接动手写。1. 先看一个能把人逼疯的场景斐波那契的 n 到了 1e181.1 普通递推的极限在哪里斐波那契数列大概是很多人第一次接触递推时写的东西long long fib(int n) { if (n 1) return n; long long a 0, b 1, c; for (int i 2; i n; i) { c a b; a b; b c; } return b; }这个小函数在 n 10⁶ 左右还能跑n 10⁸ 就已经明显卡顿n 10¹² 基本等于让你等上大半天。原因很简单它需要把从 1 到 n 的每一项都算出来时间是 O(n)。而很多题目刻意把 n 拉到 10¹⁸就是为了逼你换思路。那能不能像快速幂算 a^b 那样把 2¹⁸ 次循环变成几十次乘法呢能但前提是先把“递推”改写成“矩阵连乘”的形式。这个改写就是整个矩阵快速幂的起点。1.2 什么样的问题值得用矩阵快速幂根据我自己的经验遇到下面几类情况第一反应就该往矩阵快速幂上靠。第一类是线性递推序列。典型代表是斐波那契以及任何形如 F(n) c₁F(n-1) c₂F(n-2) ... 的常系数线性递推。这类题目在算法竞赛里出现频率极高面试中也常拿来做“你会不会优化”的考察点。第二类是图上的路径计数。如果用邻接矩阵 A 表示一张有向图那么 A^k 的第 i 行第 j 列元素恰好表示从节点 i 出发走恰好 k 步到达节点 j 的路径条数。这个性质看起来像魔法实际上就是矩阵乘法的定义在作祟后面我会展开讲。第三类是状态转移型动态规划。有些 DP 的状态转移方程是线性的比如“第 n 秒的状态 常数矩阵 × 第 n-1 秒的状态”典型如马尔可夫链。这种问题同样可以套矩阵快速幂把“一秒一秒地转移”变成“一次性转移 K 步”。不过有一点必须提前说清楚矩阵快速幂只适合“线性”递推。如果递推式里有 F(n-1) × F(n-2) 这种非线性项矩阵是表达不出来的就别硬套了。2. 拆开看原理矩阵乘法 幂的二进制分解2.1 矩阵乘法是怎么乘的很多人被“矩阵”两个字吓住其实你在算法题里用到的矩阵绝大多数是 2×2、3×3 这种小型矩阵。矩阵乘法的定义对于一个大小为 a×b 的矩阵 A 和另一个大小为 b×c 的矩阵 B结果 C 是 a×c 矩阵C[i][j] Σₖ A[i][k] × B[k][j]这个式子的直观含义就是结果矩阵里第 i 行第 j 列的数是 A 的第 i 行和 B 的第 j 列逐位相乘再求和。两个矩阵能相乘的前提是 A 的列数等于 B 的行数。在快速幂场景里我们一般只和方阵打交道也就是行列数相等所以不需要太担心维度不匹配的问题。2.2 斐波那契数列怎么变成矩阵乘法以斐波那契为例。递推式是F(n) F(n-1) F(n-2)如果只看这一项你可能会困惑这跟矩阵有什么关系关键一步是把状态“打包”成向量。设状态向量为S(n) [ F(n) , F(n-1) ]那么 S(n) 的两个分量分别等于F(n) 1 × F(n-1) 1 × F(n-2) F(n-1) 1 × F(n-1) 0 × F(n-2)写成矩阵形式就是[ F(n) ] [ 1 1 ] × [ F(n-1) ] [ F(n-1) ] [ 1 0 ] [ F(n-2) ]也就是说如果我们记转移矩阵 M [[1,1],[1,0]]那么 S(n) M × S(n-1)。继续往下推S(2) M × S(1) S(3) M × S(2) M² × S(1) ... S(n) M^(n-1) × S(1)其中 S(1) [F(1), F(0)]如果按 F(0)0、F(1)1 的惯例就是 [1, 0]。这样一来求 F(n) 就变成了求 M 的 n-1 次幂再和初始状态向量相乘。你可能会问把一个循环问题变成矩阵幂这不是更麻烦了吗关键就在于“求幂”这件事可以非常快。2.3 快速幂的本质把指数拆成二进制快速幂算法本身是个很经典的思想。以 a^13 为例13 的二进制是 1101也就是a^13 a^(841) a^8 × a^4 × a^1我们可以从 a 开始不断对底数自乘得到 a¹, a², a⁴, a⁸, ...然后根据二进制位决定要不要乘进结果里。这样原本需要 13 次乘法现在只需要 4 次自乘加 3 次结果相乘时间复杂度 O(log n)。这个思想完全可以直接搬到矩阵上。矩阵乘法虽然比数字乘法复杂但它同样满足结合律即 (A × B) × C A × (B × C)。满足结合律就意味着我们可以用和数字快速幂一模一样的框架来处理矩阵幂。唯一要注意的是矩阵乘法不满足交换律所以“结果矩阵 × 底数矩阵”的顺序必须固定不能随手换。3. C 实现从结构体到运算符重载3.1 先写一个矩阵结构体我在竞赛里最常用的方式是定义一个结构体内部用一个二维数组存储数据。考虑到矩阵快速幂通常只需要小矩阵直接开固定大小数组是最省心、性能也最好的方案#include iostream #include cstring using namespace std; typedef long long ll; const ll MOD 1000000007LL; const int MAXN 105; struct Matrix { int n, m; ll a[MAXN][MAXN]; Matrix(int n 0, int m 0) : n(n), m(m) { memset(a, 0, sizeof(a)); } };这里我保留了 n、m 表示行列数虽然方阵场景里 n 恒等于 m但保留两个字段能让代码更通用。构造函数里用 memset 清零非常重要避免后续累加时读到未初始化的垃圾值。3.2 重载乘法运算符为了让快速幂代码读起来像普通数字那样自然我会重载 operator*。实现就是三层的矩阵乘法循环Matrix operator*(const Matrix x, const Matrix y) { Matrix res(x.n, y.m); for (int i 0; i x.n; i) { for (int k 0; k x.m; k) { if (x.a[i][k] 0) continue; for (int j 0; j y.m; j) { res.a[i][j] (res.a[i][j] x.a[i][k] * y.a[k][j]) % MOD; } } } return res; }这里我特意调整了循环顺序外层是 i中层是 k内层是 j。从矩阵乘法的定义式看标准的写法应该是 i、j、k但 i、k、j 的顺序可以利用计算机缓存的局部性原理稍微提高一点速度。中间加了个 if (x.a[i][k] 0) continue对于稀疏矩阵效果明显虽然斐波那契矩阵不稀疏但养成这个习惯没坏处。3.3 快速幂主函数矩阵快速幂和普通快速幂的框架完全一致Matrix mat_pow(Matrix base, ll exp) { int n base.n; Matrix res(n, n); for (int i 0; i n; i) res.a[i][i] 1; // 单位矩阵 while (exp 0) { if (exp 1) { res res * base; } base base * base; exp 1; } return res; }res 初始化为单位矩阵是极其关键的一步。单位矩阵相当于数字里的 1任何矩阵乘以单位矩阵都等于它自己。如果初始化为全 0 矩阵那结果也会是全 0这是初学矩阵快速幂时最容易犯的错。想想数字快速幂里为什么 res 初始化为 1这里就完全对应。3.4 完整示例求第 n 个斐波那契数有了上面的结构体、乘法和幂函数求斐波那契就非常简单了int main() { ll n; cin n; if (n 0) { cout 0 endl; return 0; } Matrix M(2, 2); M.a[0][0] 1; M.a[0][1] 1; M.a[1][0] 1; M.a[1][1] 0; Matrix Mn mat_pow(M, n - 1); // 初始状态向量 S(1) [F(1), F(0)] [1, 0] // F(n) Mn.a[0][0] * 1 Mn.a[0][1] * 0 cout Mn.a[0][0] % MOD endl; return 0; }这里有个小细节值得注意我们求的是 M^(n-1)不是 M^n。因为状态向量是从 S(1) 开始的如果从 S(0) 开始定义就变成 M^n 了。这个边界非常容易踩坑务必根据自己定义的 F(0) 和 F(1) 来确定指数。用几个小数据验证一下。按 F(0)0, F(1)1n斐波那契值Mn.a[0][0] 计算结果111211322433555105555至少在小数据上是完全吻合的。4. 第一个落地案例邻接矩阵的 K 步可达路径数4.1 邻接矩阵相乘为什么能数路径斐波那契是矩阵快速幂最常见的教学例子但我觉得最能体现“矩阵乘法之美”的场景还得是图上的路径计数。假设一张有向图的邻接矩阵是 A。A[i][j] 1 表示 i 到 j 有一条边A[i][j] 0 表示没有。那么 A² 的每个元素A²[i][j] Σₖ A[i][k] × A[k][j]这个式子的含义非常直白i 到 j 走两步的路径必须先走到某个中间节点 k再从 k 走到 j。A[i][k] 表示第一步存在A[k][j] 表示第二步存在两者相乘为 1 就说明存在一条经过 k 的两步路径。把所有可能的 k 加起来就是 i 到 j 走两步的路径总数。再用 A³ A² × A 推下去就能得到递推式走 K 步的路径数等于走 K-1 步的路径数再乘一次邻接矩阵。写得更紧凑一点A^k[i][j] i 到 j 恰好走 k 步的路径总数这就是矩阵快速幂在图论中的经典应用。4.2 一个三节点的例子构造一个最简单的有向图节点 1 到 2节点 2 到 3节点 1 到 3节点 3 到 1。邻接矩阵为A [ 0 1 1 ] [ 0 0 1 ] [ 1 0 0 ]手算 A²A²[1][3] A[1][1]A[1][3] A[1][2]A[2][3] A[1][3]A[3][3] 0 1×1 0 1这个 1 对应的路径是 1→2→3。再看 A²[3][2]A²[3][2] A[3][1]A[1][2] A[3][2]A[2][2] A[3][3]A[3][2] 1×1 0 0 1对应路径是 3→1→2。通过这种方式两步路径的数量全部可以被 A² 精确表达。4.3 代码骨架代码其实和第 3 节几乎一样只是把矩阵 M 换成图的邻接矩阵int main() { int n; ll k; cin n k; Matrix A(n, n); for (int i 0; i n; i) { for (int j 0; j n; j) { cin A.a[i][j]; } } Matrix Ak mat_pow(A, k); for (int i 0; i n; i) { for (int j 0; j n; j) { cout Ak.a[i][j] % MOD ; } cout endl; } return 0; }如果题目问的是“是否存在走 K 步到达的路径”那就在矩阵里用布尔值 0/1乘法时把加法换成逻辑或、乘法换成逻辑与。不过更常见的是问方案数这时候记得取模因为不取模的路径数量在 K 较大时会爆炸式增长long long 根本扛不住。4.4 这一类问题的变体除了最简单的路径计数邻接矩阵快速幂还能处理很多变体。比如给每条边一个权值问“走 K 步且总权值最小是多少”这时候把加法换成 min乘法换成加法就成了矩阵的广义乘法。又比如装备了“恰好 K 步”这个限制后很多原本要用 BFS 多层状态的题目也能用矩阵压缩状态转移。我个人觉得这类题的核心并不在于矩阵快速幂本身而在于你能否想到“用矩阵去描述一步转移”。5. 再往上走一步任意常系数线性递推的状态矩阵构造5.1 从“两个依赖项”到“多个依赖项”斐波那契只有 F(n-1) 和 F(n-2) 两项所以状态向量是 2 维。但如果递推式变成F(n) 2F(n-1) 3F(n-3) 5怎么办方法仍然是把状态向量拉长把所有需要用到的“历史项”和常数项都装进去。既然递推式依赖 F(n-1) 和 F(n-3)那就需要保存 F(n-2)否则无法从上一状态滚到下一状态。这就要一个 4 维状态向量S(n) [ F(n), F(n-1), F(n-2), 1 ]最后一个分量是常数 1专门用来吸收递推式末尾的 5。状态转移要满足S(n) M × S(n-1)逐个分量对齐F(n) 2F(n-1) 0F(n-2) 3F(n-3) 5×1 F(n-1) 1F(n-1) 0F(n-2) 0F(n-3) 0×1 F(n-2) 0F(n-1) 1F(n-2) 0F(n-3) 0×1 1 0F(n-1) 0F(n-2) 0F(n-3) 1×1所以转移矩阵是 4×4M [ 2 0 3 5 ] [ 1 0 0 0 ] [ 0 1 0 0 ] [ 0 0 0 1 ]注意里面那个 5 放在第一行最后一列这正是常数项的作用位置。去掉常数项时把最后一列全设为 0 即可。这个“拉长状态向量来吸收常数和缺失项”的技巧是构造矩阵最核心的一步。5.2 多个递推式互相耦合有时候题目会给不止一个递推关系比如G(n) G(n-1) 2H(n-1) H(n) 3G(n-1) 4H(n-1)这种两个递推式互相依赖的情况用数组模拟也能做但一旦 n 到了 1e18就必须矩阵化了。状态向量取 S(n) [G(n), H(n)]转移矩阵是M [ 1 2 ] [ 3 4 ]于是 S(n) M^(n-1) × S(1)。本质上任何一个“当前状态完全由上一时刻状态线性映射而来”的系统都可以写成矩阵形式。你要做的就是把这种线性映射关系找出来填到矩阵里。5.3 构造矩阵的经验法则根据我自己的实操经验构造状态矩阵时有几条可以照着做先确定递推式的阶数 k也就是最多依赖到 F(n-k) 的哪一项。状态向量里除了 k 个历史项如果需要常数项就再加一个分量 1。转移矩阵的第一行严格按递推式系数填写缺少的项填 0常数项填递推式尾部的常数。第一行之外的行作用是把向量整体“向下平移”第 i 行只在第 i-1 列写 1其余全是 0保证 S(n-1) 的旧分量能在 S(n) 中保留。常数分量对应的行最后一位写 1确保常数项一直保持为 1。6. 实测中必须注意的坑溢出、边界和性能6.1 单位矩阵初始化最常见的翻车点我第一次写矩阵快速幂时结果全是 0排查了半天发现 res 初始成了全 0 矩阵。这个错误太经典了以至于我现在看到别人代码里结果全零第一反应就是去看单位矩阵。单位矩阵的定义很简单对角线为 1其余为 0。写的时候用 for(int i0;in;i) res.a[i][i]1 即可千万别直接 memset(res, 0, sizeof(res)) 就完事。6.2 乘法溢出的层级问题取模运算里有个隐蔽的坑。假设 MOD 1e97两个小于 MOD 的数字相乘结果是 1e18 量级已经非常接近 long long 的极限 9.22e18勉强能放下。但如果再乘一次比如 x.a[i][k] * y.a[k][j] % MOD 这一步没取模而是直接累加到 res.a[i][j] 上累加出去的数可能超过 long long 承受范围。我常用的写法是res.a[i][j] (res.a[i][j] x.a[i][k] * y.a[k][j] % MOD) % MOD;先让乘法结果取一次模再和已有结果相加最后再取一次模。这样每一步的操作数都控制在 MOD 级别比较安全。如果题目给的 MOD 特别大接近 1e18乘法本身就可能溢出那就得用 __int128 过渡或者写一个快速乘法函数。普通场景下1e9 这个量级配合上面的写法足够稳定。6.3 快速幂指数的边界到底是谁不少人在斐波那契问题里纠结到底传 n 还是 n-1。这里有一个笨但可靠的验证方法先用暴力循环算出前几个 F(n)然后用矩阵算同样的 n比对一下。我个人偏好把初始状态定义为 S(0) [F(0), F(1)]然后用 S(n) M^n × S(0)这样指数就是 n遇到 n0 时返回单位矩阵即可边界更直观。但无论采用哪种定义请务必在代码里显式处理 n0 的情况。很多人只顾着算 M^(n-1)忘了 n0 时指数是 -1会直接崩溃。6.4 固定数组和 vector 的性能对比矩阵快速幂的复杂度是 O(k³ log n)其中 k 是矩阵边长。k 2 或 3 时用 vectorvector 和固定数组差别不大但当 k 到 50 甚至 100三层循环里的内存访问就非常关键了。竞赛中我几乎只用固定数组因为在栈上连续分布的数据对缓存更友好vector 的多次间接寻址会带来不小的常数开销。如果事先不知道最大 k可以用 const int MAXN 105 这种固定上限然后给结构体里放一个 n 字段记录实际的矩阵维度乘法时只遍历 0 到 n-1。这样既保有固定数组的性能又不至于写死每个矩阵都是 100×100 导致浪费。6.5 验证模板的一个高效方法矩阵快速幂代码虽然短但矩阵乘法、指数循环、取模位置任何一个细节出错结果都可能面目全非。所以我会在正式提交前专门写一个暴力的朴素递推函数让 n 取 0 到 20 的小范围逐项对拍。这一步看起来笨但能省下大量调试时间。尤其是当你换了模数或改过矩阵维度之后对拍能立刻抓出问题。提示在写矩阵快速幂时如果结果需要取模我建议矩阵内部的乘法也全部走模运算而不是在最后统一取模。因为中间结果可能膨胀得非常快等到最后再处理就晚了。7. 一点更进阶的优化思路减少矩阵乘法次数7.1 二进制幂里能不能省几次乘法快速幂的常规写法是while(exp) { if(exp 1) res res * base; base base * base; exp 1; }。这个写法里即使当前二进制位是 0base 的自乘也一定执行。有人会试图通过预计算所有 base 的 2^k 次幂再组合相乘从理论上减少一些乘法次数。但实际效果非常有限因为省去的只是少量结果矩阵乘法而 base 自乘是必须做的。在竞赛里常规写法已经足够快。7.2 利用广义矩阵乘法处理 min-plus 问题前面提过把矩阵乘法的加法换成 min、乘法换成加法就可以解决“最短 K 步路径”之类的优化问题。这其实是把矩阵快速幂思想迁移到了半环结构上。只要运算满足结合律快速幂框架就能继续用。这个扩展在刷题时偶尔会遇到了解原理能帮你多一条解题思路。7.3 什么时候不要用矩阵快速幂矩阵快速幂的性能瓶颈是 O(k³ log n)这个 k 严格取决于状态向量的维度。一旦 k 超过 100即使 log n 只有 60总共也要执行 100³ × 60 ≈ 6000 万次乘法运算考虑到取模的开销已经在超时边缘了。如果递推式的阶数很大通常会改用线性递推的优化算法比如 Kitamasa 算法它能做到 O(k² log n) 甚至 O(k log k log n)。不过刷题和面试中题目设计者一般不会让 k 大到这个程度所以我不会优先推荐上更复杂的算法但至少要知道矩阵快速幂不是万能的。8. 几个可以直接抄的小模板和个人习惯8.1 一个精简的通用模板平时我用得最多的模板其实就是把第 3 部分的结构体、乘法和快速幂整合到一起。为了防止每次重写时遗漏细节我会把它保存在本地代码片段里。模板的关键部分已经在这篇文章里给出过我再补充一个针对“多点查询”的优化如果同一张转移矩阵要乘很多次不同的指数可以先把底数的 2^k 次幂预计算成一个数组每次查询只需要对二进制位逐位相乘。虽然预计算要花 O(k³ log MAXN) 的时间但多次查询时能省去重复的 base 自乘。8.2 代码风格上的小建议结构体里直接开 MAXN×MAXN 的二维数组会让代码有点“丑”但对竞赛和面试来说这是最可靠的选择。相对地如果在实际工程里使用我会封装成一个 Matrix 类把维度、模数都作为构造参数乘法函数也改成成员函数以降低误用概率。对于刷题场景我反而建议保持精简因为结构体越短越不容易在紧张的比赛环境里写出隐藏 bug。8.3 调试矩阵快速幂时我最常做的事最后分享一个非常个人的习惯每写完一次矩阵快速幂我会先初始化一个 2×2 的矩阵让它是 [[1,2],[3,4]]然后写一个小测试分别计算 M¹、M²、M³和手算的结果对照。这么做不是为了验证数学公式而是为了确认乘法运算符重载和取模逻辑没有错。只要这个小测试通过再换成真正的递推矩阵问题基本就解决了。这个习惯帮我避免了很多次“写完但不知道哪里错了”的尴尬局面。对于矩阵快速幂我实际用下来最大的感受是它不只是一个算法模板更是一种把递推转化为状态空间映射的思维工具。如果你能熟练地写出 2×2 的斐波那契矩阵那你离解决“图上 K 步路径”“马尔可夫链状态分布”“多个递推式耦合”这类问题其实只差一层窗户纸。多拿几道题练手尤其是亲手把递推公式在纸上展开成矩阵再落成 C 代码这个能力就算真正长在你身上了。