企业官网建设流程全解析

你有没有遇到过这种情况：一个递推式，比如斐波那契数列，第 (10^{18}) 项你要算出来，取模 (10^9+7)。暴力循环肯定不行，每秒跑一亿次也得等几百年。通项公式倒是好看，但里面有根号和浮点除法，大数精度一塌糊涂。这时候就该请出矩阵快速幂了。

这东西是线性代数和算法竞赛的第一次深度碰撞，也是我当年入门时被反复“虐”后彻底吃透的一个知识点。它的核心思路非常朴素：把递推关系写成矩阵乘法，然后沿用普通快速幂的二进制分治思想，把 O(n) 的时间压到 O(log n)。一旦掌握，斐波那契这种入门题只是开始，线性递推、图论路径计数、马尔可夫链状态演化、部分动态规划优化，全都可以用同一套框架解决。

这篇博文会从线性代数基础开始讲，结合完整代码和可验证的例子，把矩阵快速幂的构造、实现、扩展和排错一次性讲透。适合正在学算法的新手，也适合已经会写模板但总在应用上卡壳的朋友。

1. 这东西到底解决什么问题

先别急着看代码，先把问题场景看清楚。矩阵快速幂不是炫技，它解决的是一类非常具体的问题：已知一个线性递推关系，要求第 n 项，且 n 极大。

1.1 从一道题说起：暴力为什么扛不住

斐波那契数列是最典型的例子。定义很简单：

[ F_0 = 0,\quad F_1 = 1,\quad F_n = F_{n-1} + F_{n-2} ]

现在问 (F_{10^{18}} \bmod (10^9+7)) 是多少。你可能会想，O(n) 的循环不是挺快吗？但 (10^{18}) 次循环，即使每次操作只需要 1 纳秒，也要超过 31 年。所以 O(n) 在 (n) 达到 (10^{18}) 这个量级时，和 O(2^n) 没有本质区别，都是不可接受的。

还有一条路是通项公式。斐波那契通项里有 ((1+\sqrt{5})/2)，你要在模意义下处理这个无理数，就得引入二次剩余，实现起来很麻烦，而且很容易被浮点误差坑。通项公式在数学分析里很漂亮，但在工程实现里并不是首选。

这时就该转换视角：这个递推关系本质上是一个向量被同一个矩阵反复作用。只要这个矩阵的幂能快速算出来，递推项就能直接求。这就是矩阵快速幂的切入点。

1.2 普通快速幂到矩阵快速幂的跨越

普通快速幂大家都很熟。计算 (a^b)，不是连乘 b 次，而是把 b 拆成二进制，通过不断平方来减少乘法次数。比如 (a^{13} = a^{8+4+1} = a^8 \cdot a^4 \cdot a^1)，只需要 4 次乘法。

矩阵快速幂的逻辑一模一样，只是把“数”换成了“矩阵”：

[ M^k = M^{2^0 \cdot bit0} \cdot M^{2^1 \cdot bit1} \cdot M^{2^2 \cdot bit2} \cdots ]

每次迭代把当前矩阵平方，遇到二进制位为 1 就把结果累乘进去。

但这里有一个很关键的前提：为什么可以这样拆？因为矩阵乘法满足结合律，也就是 ((A \times B) \times C = A \times (B \times C))。有了结合律，我们才能任意调整乘法顺序，把一个巨大的幂拆成若干个小幂的乘积。而矩阵乘法不满足交换律，所以拆出来的因子顺序不能乱，好在求幂时所有因子都是同一个矩阵 (M)，天然不会出现顺序问题。这一点在后面构造转移矩阵时还要再强调。

2. 线性代数基础：先把矩阵乘法吃透

矩阵快速幂最底层的零件是矩阵乘法。如果这一步建得不稳，后面所有东西都是空中楼阁。

2.1 矩阵乘法规则回顾

两个矩阵 (A)（(m \times n)）和 (B)（(n \times p)）相乘，得到矩阵 (C)（(m \times p)），其中：

[ C_{i,j} = \sum_{k=1}^{n} A_{i,k} \cdot B_{k,j} ]

也就是说，C 的第 i 行第 j 列，等于 A 的第 i 行和 B 的第 j 列对应位置相乘再求和。内层维度必须一致，否则无法相乘。

用生活化的类比：矩阵乘法像是一个工厂里两条流水线的衔接。A 的每一行代表一种原料配比要求，B 的每一列代表一种机器的产出参数，两者相乘得到最终产出。这也是为什么矩阵乘法天然适合描述“经过多个步骤后的组合变化”。

这个规则是整个应用的地基。比如图论里，邻接矩阵相乘时，每一项 (C_{i,j} = \sum A_{i,k} \cdot B_{k,j}) 其实是在枚举所有中间节点 k，把所有“从 i 到 k，再从 k 到 j”的方案数累加起来。矩阵乘法的意义在不同场景下会“长”出不同的解释，这也是它好用得不得了的原因。

2.2 单位矩阵：快速幂里最容易忽略的细节

整数乘法里有 1，任何数乘 1 等于它本身。矩阵里也有对应的东西，就是单位矩阵 (I)，它是一个方阵，主对角线全为 1，其余全是 0。

[ I = \begin{bmatrix} 1 & 0 & 0 \ 0 & 1 & 0 \ 0 & 0 & 1 \end{bmatrix} ]

单位矩阵在快速幂里的角色是乘法初始结果。普通快速幂里，结果 res 初始化为 1；矩阵快速幂里，res 初始化为单位矩阵 I。我第一次写代码时就在这里栽过跟头，把 res 初始化成全零矩阵，结果算出来全是 0，那真是怎么查都查不出毛病。

单位矩阵的存在也很有意思：它是矩阵乘法的“中性元素”，就像加法里的 0、乘法里的 1。构造单位矩阵时不要手写一堆 if，直接用循环判断 i == j 赋最稳妥，高阶维数时也不容易错。

2.3 结合律是快速幂的地基

为什么快速幂能把 (M^k) 拆开算？核心就是结合律。它意味着当你有一长串矩阵相乘时，先算哪一部分都不影响最终结果。快速幂正是利用这一点，把 (M^{2^0}, M^{2^1}, M^{2^2}, \ldots) 这些预先算好的结果，按二进制位重新组合成最终答案。

这里有个细节值得单独拎出来：矩阵乘法没有交换律。也就是说 (A \times B \neq B \times A)。所以写快速幂时，结果矩阵 res 和当前矩阵 base 的乘法顺序不能随便调换。我习惯统一写res = mul(res, base)，保持“res 在左、base 在右”。如果你中途换了风格写成res = mul(base, res)，虽然有时候碰巧结果一样，但完全没有数学保障，尤其是状态向量和转移矩阵配合时，方向错了会让结果整个崩掉。

所以在实现之前，先把乘法方向约定清楚：要么全程行向量左乘，要么全程列向量右乘。混用是调试时最折磨人的错误来源。

3. 核心实操：用矩阵快速幂实现斐波那契

理论说了半天，现在开始动手。这一节我会从递推式到矩阵构造，再到完整代码，一步步把过程拆开。

3.1 怎么把递推公式变成矩阵形式

斐波那契递推是 (F_n = F_{n-1} + F_{n-2})。这个式子的问题在于，它需要前两项来计算当前项，是个“二阶”依赖。为了把它变成矩阵乘法，我们可以把相邻两项放进一个向量：

[ v_n = \begin{bmatrix} F_n \ F_{n-1} \end{bmatrix} ]

然后观察从 (v_{n-1}) 到 (v_n) 发生了什么：

[ \begin{bmatrix} F_{n} \ F_{n-1} \end{bmatrix}

\begin{bmatrix} F_{n-1} + F_{n-2} \ F_{n-1} \end{bmatrix}

\begin{bmatrix} 1 & 1 \ 1 & 0 \end{bmatrix} \begin{bmatrix} F_{n-1} \ F_{n-2} \end{bmatrix} ]

所以转移矩阵就是：

[ M = \begin{bmatrix} 1 & 1 \ 1 & 0 \end{bmatrix} ]

于是：

[ v_n = M^{n-1} \cdot v_1 ]

其中 (v_1 = [F_1, F_0]^T = [1, 0]^T)。如果初始条件从 (F_0, F_1) 开始，也可以写成 (v_n = M^n \cdot v_0)，其中 (v_0 = [F_0, F_{-1}]^T)。我建议每个问题都自己推一遍初始向量，不要死记公式，因为初始下标不同，指数就会差 1，结果就会错位。

这里为什么要用一个 2 维向量而不是直接用单个数字？因为单个数字无法同时记录“上两项的信息”。把状态空间扩展成向量，矩阵就可以在这个状态空间里做线性变换，这就是矩阵能介入递推问题的数学根源。

3.2 代码实现：迭代版矩阵快速幂

下面是完整的 Python 实现。为了性能演示，我把矩阵乘法写成二维列表操作。实际比赛或工程中，可以用 numpy 或固定结构体，但算法骨架是不变的。

MOD = 10**9 + 7 def mat_mul(A, B): n = len(A) m = len(B[0]) p = len(B) C = [[0] * m for _ in range(n)] for i in range(n): for k in range(p): if A[i][k] == 0: continue aik = A[i][k] for j in range(m): C[i][j] = (C[i][j] + aik * B[k][j]) % MOD return C def mat_pow(A, exp): n = len(A) # 单位矩阵初始化 res = [[1 if i == j else 0 for j in range(n)] for i in range(n)] base = A while exp > 0: if exp & 1: res = mat_mul(res, base) base = mat_mul(base, base) exp >>= 1 return res def fib(n): if n == 0: return 0 M = [[1, 1], [1, 0]] Mn = mat_pow(M, n - 1) # Mn * [F1, F0]^T = [Fn, F_{n-1}]^T return Mn[0][0] # 因为 v1 = [1, 0]^T，结果第一行第一列就是 Fn

这段代码里有一个关键点：mat_mul里的循环顺序写成了i -> k -> j，并且先判断A[i][k]是否为零，这是为了利用稀疏性提速，同时对缓存更友好。写矩阵乘法时，我们习惯性地按“行乘列”写三重循环，但实际执行时，内存访问模式会影响常数。i -> k -> j的写法让 B 的访问在内存中更连续，这在矩阵比较大的时候体验非常明显。

3.3 复杂度分析：O(n^3 log k) 是怎么来的

矩阵快速幂的总复杂度由两部分组成：每次矩阵乘法是 O(n^3)（n 是矩阵边长），快速幂需要进行 O(log k) 次乘法。所以总体是 O(n^3 log k)。

以斐波那契为例，n = 2，于是复杂度 O(2^3 log k) = O(8 log k)，也就是 O(log k)。对比原来的 O(k)，这可以说是天壤之别。当 k=10^18 时，log2(k) 大约是 60，也就是说哪怕矩阵乘法常数大一点，也只是几百次 2x2 乘法的事，毫秒级出结果。

空间复杂度是 O(n^2)，用来存矩阵。这里有一个容易误判的地方：快速幂需要保留中间结果吗？不需要。我们只需要二进制分解过程中不断平方的 base 矩阵，以及不断累乘的 res 矩阵。所以每一轮迭代都只需要两个矩阵的内存，和 log k 无关。

3.4 正确性验证：怎么知道自己没写错

写完代码第一步别拿大 n 直接跑，先拿小 n 验证。斐波那契前几项：0, 1, 1, 2, 3, 5, 8, 13, 21, 34。你跑fib(10)应该等于 55。如果对不上，优先检查三个地方：

1. 单位矩阵是否初始化正确，有没有把i == j写成i != j。
2. 转移矩阵 M 的顺序是否和状态向量定义一致。我们是[F_n, F_{n-1}]^T = M * [F_{n-1}, F_{n-2}]^T，如果定义的是行向量[F_n, F_{n-1}]左乘矩阵，结果会差一个转置，数值自然不对。
3. 指数是否偏移。n=1 时应该直接返回 1，如果统一套mat_pow(M, n)，就会变成求 M^1 而不是 M^0，结果错位到第 2 项。

我自己的习惯是：先用暴力递推写一个fib_slow(n)，然后对比fib(0)到fib(30)的全部结果，做成一个自动化小测试。三十项以内暴力计算毫无压力，却能覆盖几乎所有构造错误，比直接跑大样例靠谱得多。

4. 从斐波那契扩展到更多场景

矩阵快速幂的威力不在于斐波那契本身，而在于它背后的抽象：一切线性递推，都可以塞进这个框架。下面这几个场景非常典型，值得一个个看过去。

4.1 线性递推通吃：从二阶到高常系数

假设递推式是：

[ a_n = 2a_{n-1} + 3a_{n-2} + 4a_{n-3} ]

目标是求 (a_n)。做法是扩大状态向量：

[ v_n = \begin{bmatrix} a_n \ a_{n-1} \ a_{n-2} \end{bmatrix} ]

然后观察递推关系：

[ v_n = \begin{bmatrix} 2 & 3 & 4 \ 1 & 0 & 0 \ 0 & 1 & 0 \end{bmatrix} \cdot \begin{bmatrix} a_{n-1} \ a_{n-2} \ a_{n-3} \end{bmatrix} ]

看到套路没有？转移矩阵的第一行就是递推系数，下面是一个单位矩阵的“移位结构”。每一轮，原来的 (a_{n-1}) 下沉成为新的 (a_{n-2})，原 (a_{n-2}) 下沉成为新的 (a_{n-3})。这是经典的高阶线性递推状态转移矩阵构造法。

一般地，k 阶递推：

[ a_n = c_1 a_{n-1} + c_2 a_{n-2} + \cdots + c_k a_{n-k} ]

状态向量为 ([a_n, a_{n-1}, \ldots, a_{n-k+1}]^T)，转移矩阵第一行为 ([c_1, c_2, \ldots, c_k])，左上三角第二行开始是单位阵的移位形式，其余都是 0。

注意：系数顺序千万别反。常见错误是把 ([c_1, c_2, \ldots, c_k]) 写成 ([c_k, \ldots, c_2, c_1])，然后结果全错。建议先在草稿纸上展开 (v_{n+1}) 和 (M v_n)，一行一行对清楚再写代码。

4.2 图论：邻接矩阵的 k 次幂到底意味着什么

另一个非常惊艳的应用是图论里的路径计数。给定一个有向图的邻接矩阵 A，A[i][j] = 1 表示 i 到 j 有一条边。那么 (A^2[i][j]) 表示什么？

按矩阵乘法展开：

[ A^2[i][j] = \sum_{k} A[i][k] \cdot A[k][j] ]

这相当于枚举中间节点 k，把所有“i → k → j”的长度为 2 的路径数加起来。推广一下，(A^m[i][j]) 就是从 i 到 j 长度为 m 的路径总数。如果 A 里存的是边权，那乘法规约要相应调整，比如加法换成取 min，就变成 Floyd-Warshall 式的语义；如果存的是方案数，那就直接乘法和累加。

这里快速幂的意义就大了：如果要求长度为 k 的所有点对路径数，直接做 k 次循环是 O(k n^3)，而矩阵快速幂是 O(n^3 log k)。当 k 很大时，这个差距是决定性的。

不过要提醒一句：路径数量会呈指数增长，如果不取模很快就会爆。实际竞赛中通常会给一个模数，或者要求只数可行性的布尔值（此时用逻辑运算替代乘法和加法，效果等价于传递闭包的多次复合）。

4.3 概率与状态转移：马尔可夫链

再看概率场景。假设一个系统有 n 个状态，(P) 是状态转移概率矩阵，(P[i][j]) 表示从状态 i 到状态 j 的概率，每一行之和为 1。如果初始状态分布是向量 (v_0)，那么经过 1 步后的分布是 (v_1 = v_0 \times P)，经过 k 步后是：

[ v_k = v_0 \times P^k ]

当你需要算第 (10^9) 步的分布时，模拟 (10^9) 次显然不现实，但矩阵快速幂可以直接把 (P^k) 算出来。

这里用行向量还是列向量？取决于你约定的状态分布是行向量还是列向量。如果 (v) 是行向量，转移矩阵写在右边：(v' = v P)。如果 (v) 是列向量，转移矩阵写在左边：(v' = P v)。两种风格都行，但必须统一，然后整个项目一直用这个约定。我见过太多人在这一步混了，导致加一个状态之后整个结果像是随机数。

还有个实用小技巧：概率转移矩阵很多是稀疏的，也就是大多数元素是 0。这时候不要傻乎乎地用三层循环满矩阵乘法，可以先判断< 0.0000001的项跳过。虽然复杂度理论上没变，但实际常数能小好几倍。

4.4 动态规划优化视角

最后说一个更难发现的应用：把一部分 DP 优化成矩阵快速幂。

有一类 DP 的转移是“线性且周期固定”的。比如说经典的爬楼梯问题，dp[i] = dp[i-1] + dp[i-2]，和斐波那契没区别。更一般的，如果 dp[i] 只和 dp[i-1], dp[i-2], ..., dp[i-k] 线性相关，那么整个 DP 的递推就可以看成是一个向量状态下的线性变换。

举一个稍微复杂一点的例子：有些题目把 dp 数组当作向量，每一步的转移都一样，那么计算第 n 步 dp 状态就等价于计算 (dp_0 \times T^n)。这和马尔可夫链的形式本质上是一样的。

这时候，判断一道 DP 能不能用矩阵快速幂，就看两点：

1. 状态转移是否线性，即新值是否由旧状态的线性组合得到。
2. 转移次数是否巨大，小范围就直接循环了，不必用矩阵。

很多同学见到 DP 就想着“矩阵优化”，反而把小问题搞复杂了。正确姿势是：先暴力递推验证思路，再考虑优化。如果 n 只有 (10^5)，普通 DP 完全扛得住，没必要碰矩阵；只有当 n 到了 (10^{12}) 以上，矩阵快速幂才是不可替代的选择。

5. 写矩阵快速幂时踩过的坑

这一部分全是我自己实际排错的经验。矩阵快速幂代码量不大，但错误类型非常集中。下面这些问题，可以说每一条都有无数人踩过。

5.1 模运算：取模的时机和位置

在算法题里，结果通常要求对某个模数取模。取模时机看起来是小事，但做错了会拿到一个“看起来没问题”的错误答案。

我的原则是：矩阵乘法内部，每次累加都取模，而不是最后统一取模。原因很简单，中间累加结果很可能是一个巨大的数字。在 Python 里大整数不会有溢出，但速度会显著变慢；在 C++ 里则直接 UB（未定义行为），long long 一溢出就全完蛋。

for (int i = 0; i < n; i++) { for (int k = 0; k < n; k++) { if (!a[i][k]) continue; long long aik = a[i][k]; for (int j = 0; j < n; j++) { c[i][j] = (c[i][j] + aik * b[k][j]) % MOD; } } }

这里有个细节：即使aik * b[k][j]两个数都在 (10^9) 左右，乘积也到 (10^{18}) 量级，long long 勉强装得下，但如果再加到 c[i][j] 上，可能就爆了。所以严格的写法是每乘一次就取模，再把结果加到 c[i][j] 上取模。这样每一步结果都控制在 MOD 以内，不会累积。

5.2 单位矩阵初始化错误

单位矩阵初始化成零矩阵，是新手最常见的问题。零矩阵乘任何矩阵都还是零矩阵，结果就是全 0，而且跑得快，不报错，很容易被忽略。

检查方法不是盯着代码看，而是写一个自测：mat_pow(A, 1)应该等于 A，mat_pow(A, 0)应该等于单位矩阵。这两种情况都不满足的话，说明初始化或者 mul 函数有问题。

另外，不同编程语言里单位矩阵的写法不一样。C++ 里可以用嵌套 for 赋值，也可以用memset清零后再把对角线置 1。Python 里可以用列表推导式，但要注意每个子列表都要独立创建，不能在循环里复制同一个 list，否则修改一行会影响到其它行。这个 Python 坑我用[[0] * n] * n踩过一次，改一行全变了，排查了半小时。

5.3 递归深度与迭代选择

快速幂有两种实现，递归和迭代。递归版本代码简洁：

mat qpow(mat a, long long b) { if (b == 0) return unit(a.n); mat tmp = qpow(a, b >> 1); tmp = tmp * tmp; if (b & 1) return tmp * a; return tmp; }

但问题在于递归深度是 O(log k)，理论上并不深，k 到 (10^{18}) 时也才 60 层左右，好像没什么问题。然而在实际比赛环境中，如果每层都拷贝矩阵，资源开销很大。更重要的是，如果某个平台栈空间有限，60 层递归虽然不深，但矩阵拷贝的临时对象多，可能触发内存问题。

我个人的习惯是一律写迭代版。迭代版代码同样短，不需要考虑递归，也不会有栈溢出风险，还能在同一轮循环里顺手做二进制位判断。模板一旦固定下来，以后写任何快速幂都能直接套。

5.4 矩阵维度与下标混乱

矩阵维度搞错是一个很无语但很高发的问题。有些时候，转移矩阵是 (n \times n)，状态向量是 (n \times 1) 或 (1 \times n)，两者维度必须匹配。结果矩阵 res 的维度应等于转移矩阵的维度，不是状态向量的维度。写mat_mul时如果把 res 和 base 的维度混在一起，极端情况下会出现运行时错误，在 Python 里可能直接 IndexError，在 C++ 里就是越界访问。

下标从 0 开始还是从 1 开始，也要统一。如果你在纸上用 1-index 推导，写代码时却忘了减 1，所有下标错位，结果会让人完全摸不着头脑。建议先确定状态向量的定义，把它写在注释里，然后所有推导都围绕这个注释进行。

5.5 性能与内存优化小技巧

矩阵比较大时，常数优化就不是可选项而是必须项。

第一个优化是跳过 zero 乘法。转移矩阵和 base 在迭代过程中会很快变“密”，但相乘时如果 A[i][k] 为 0，可以跳过内层 j 循环，省掉一大批不必要操作。这个优化的收益取决于矩阵的稀疏程度，稀疏矩阵下收益极其可观。

第二个优化是固定维度模板。在 C++ 里，与其用vector<vector<long long>>，不如用static long long a[N][N]，或者直接写一个固定 N 的结构体。好处是内存连续，访问速度大幅提升，还能避免 vector 嵌套带来的多次指针跳转。模板里 N 是编译期常量，可以这样定义：

template<int N> struct Matrix { long long a[N][N]; Matrix() { memset(a, 0, sizeof(a)); } void unit() { memset(a, 0, sizeof(a)); for (int i = 0; i < N; i++) a[i][i] = 1; } Matrix operator*(const Matrix& other) const { Matrix c; for (int i = 0; i < N; i++) { for (int k = 0; k < N; k++) { if (a[i][k] == 0) continue; for (int j = 0; j < N; j++) { c.a[i][j] = (c.a[i][j] + a[i][k] * other.a[k][j]) % MOD; } } } return c; } };

第三个优化是乘法子块或循环展开。不过这通常是在 N 较大时才值得做。2x2 或 3x3 的小矩阵，直接裸乘法都行，优化反而增加代码复杂度。写算法题时，先确认 N 的范围再决定优化等级，不要一开始就上重型优化。

6. 关于这个技术的个人心得

讲了这么多，最后聊几句我自己的实践经验，希望能帮你少走弯路。

矩阵快速幂的价值有很大一部分来自“构造转移矩阵”的过程。不要只看模板，一定要亲手把斐波那契、三阶递推、邻接矩阵路径计数这三类题各推一遍。推完之后你会发现，大多数能用矩阵快速幂的题，本质上都是这三类的变体：要么状态扩展成向量，要么矩阵带有特殊结构。

在实际项目中，我很少遇到需要手写矩阵快速幂的场景，因为很多框架里已经有现成的线性代数库。但理解这个算法有一个额外的好处：你会更清醒地意识到，很多看似复杂的问题，本质上就是“同一个变换反复执行”。这种抽象思维比单纯记住一个模板更有价值。

如果你在学的时候反复卡壳，我建议给每个构造步骤配一个“最小测试用例”。不要只想不做，拿 n=4 或 n=5 的实际小例子，用笔手算一遍矩阵乘法，再和代码输出对比。这类问题最容易出错的不是算法逻辑，而是构造方案的某个小细节，比如向量顺序、指数偏移、二维 vs 一维。小例子能把这些错误全部揪出来。

最后再分享一个小技巧：写完模板后，把它封装成一个独立函数，然后长期保存。以后遇到新题，如果识别出递推/路径计数/状态转移的特征，直接往里套，省下大量调试时间。封装时把取模、维度、乘法顺序都固定好，顶部注释里写清楚状态向量是行还是列、指数从几开始。这样哪怕半年后再看，也能一眼回忆起来。这个模板我用了很多年，每次用它解决问题，都还会感叹一次：线性代数在代码里的力量，真的比想象中大得多。