很多初学者第一次写求幂函数，十有八九是长这样的：
def naive_pow(a, n): res = 1 for _ in range(n): res *= a return res代码本身没什么毛病，但它经不起大指数的考验：n 一旦上亿，这个循环就跑不动了。我第一次被这个问题卡住是在某次课程设计里，要求计算 2^k 对 10^9+7 取模，k 可以大到 10^18。我当时写了上面这个循环，等了大半天没结果，后来才意识到问题不是程序错了，而是算法复杂度根本没有救。真正该做的，是把指数拆成二进制，用平方代替连乘——这就是快速幂。
这篇文章把整数快速幂和 Strassen 矩阵快速幂一起讲透。前者负责解决"大数 n 次方"的效率问题，后者把同一套分治思想搬到矩阵上，顺手解决斐波那契数列、线性递推、图上路径计数这一大类问题。无论你是准备算法竞赛、应付数据结构考试，还是工作中遇到大指数运算，这篇都值得仔细读一遍。
1. 连乘为什么慢：从二进制视角看指数
1.1 循环连乘的复杂度短板
回到 naive_pow。它的时间复杂度是 O(n)，也就是 n 次乘法。当 n 是 10、100 的时候毫无压力，可一旦 n 进入 10^9 甚至 10^18 的量级，这个复杂度就直接判死刑。做个粗略估算：现代 CPU 一秒钟大概能执行 10^9 次级别的简单运算，n=10^9 的循环至少要跑一秒以上，n=10^18 就意味着几十亿年的计算量，这在工程上等于不可计算。
所有涉及大指数的严肃场景——RSA 类的模幂运算、大数计算、组合数学里的快速求幂——全都依赖快速幂，不可能靠循环硬撑过去。理解这一点之后，再看快速幂的每个设计决策，思路都会清晰很多。
1.2 二进制拆分：13 = 8 + 4 + 1 带来的顿悟
核心洞察其实一句话就能说清：任何一个正整数 n 都可以写成 2 的幂之和。拿 13 举例，13 的二进制是 1101，也就是 13 = 8 + 4 + 1。于是：
a^13 = a^8 · a^4 · a^1
接下来注意到一个递推关系：a^1、a^2、a^4、a^8 每一项都是前一项的平方。也就是说，我不用真的去乘 13 个 a，只需要从 a^1 开始平方三次就能得到所有候选值：
- a^2 = (a^1)²
- a^4 = (a^2)²
- a^8 = (a^4)²
平方三次得到四个候选值，再把二进制位为 1 的那些挑出来相乘：a^8 · a^4 · a^1，乘三次。总共 6 次乘法，比 13 次少了一半多。指数越大，这个差距越夸张。
拿 3^13 演算一遍：13 = 1101₂，我们从 3 开始平方，3² = 9，9² = 81，81² = 6561。二进制位为 1 的是 1、4、8 这三个，所以结果 = 3 × 81 × 6561 = 1594323，和 3^13 完全一致。这个例子很朴素，但已经完整展示了快速幂的全部机制。
1.3 分治公式与复杂度推导
把上面的思路整理成递归公式：
f(a, 0) = 1 f(a, n) = f(a, n/2)² n 为偶数 f(a, n) = f(a, n/2)² · a n 为奇数每一层递归都把 n 砍一半，所以层数是 O(log n)，每层最多做两次乘法，总复杂度 O(log n)。把 n=10^18 代进去，log2(10^18) 约等于 60，也就是说大概 60 层递归、一两百次乘法就能出结果。从 10^18 次乘法降到约 10^2 次乘法，差了整整 16 个数量级。
这个公式看起来简单，但它其实是整篇文章的地基。注意推导过程从头到尾只用了乘法的结合律和"乘以 1 不变"这两条性质，完全没有依赖任何"数"特有的规则。这意味着，只要一个运算系统满足结合律并且存在单位元，这套分治办法就成立——矩阵就是这个系统里最重要的一个。
2. 整数快速幂的两种写法与三个边界坑
2.1 递归版和迭代版，我推荐后者
递归版代码最短，逻辑也最贴合分治公式：
def fast_pow_recursive(a, n): if n == 0: return 1 half = fast_pow_recursive(a, n // 2) if n % 2 == 0: return half * half return half * half * a但我实际写代码，几乎总是用迭代版。原因有两个。第一，递归调用本身有栈帧开销，虽然 log2(10^18) 只有 60 层，远不会爆掉 Python 默认的递归深度，但工程上能省一点是一点。第二，迭代版的循环结构直接对应"扫描指数的二进制位"，可读性反而更好，调试起来也直观：
def fast_pow(a, n): result = 1 while n > 0: if n & 1: result *= a a *= a n >>= 1 return result循环体里就三件事：看最低位是否为 1、决定要不要把当前的 a 乘进 result、然后让 a 自乘、n 右移一位。循环结束时 result 恰好就是所有二进制位为 1 的那些 a^(2^i) 的乘积。这个写法建议大家直接背下来，它几乎是所有快速幂变体的公共骨架。
2.2 取模版本：两个防溢出细节
真实场景里几乎没人直接求完整幂，因为 a^n 指数增长，几轮平方就超过 64 位整数上限。最常见的需求是求 a^n mod m，模运算对乘法是友好的——先乘再取模和边乘边取模结果一样，所以每步都能压模：
long long mod_pow(long long a, long long n, long long mod) { long long result = 1 % mod; a %= mod; while (n > 0) { if (n & 1) { result = (__int128)result * a % mod; } a = (__int128)a * a % mod; n >>= 1; } return result; }这里有两个特别容易翻车的细节。第一，result 初始化写成 1 % mod 而不是 1，专门处理 mod = 1 的极端情况：任何数对 1 取模都是 0，但 1 % 1 才是 0，直接初始化成 1 会返回错误结果。第二，乘法用 __int128 中转，因为两个接近 10^9 的数相乘会到 10^18 量级，long long 上限约 9.2×10^18，虽然单次乘法通常不爆，但多个累加就可能出问题，用 128 位中转最稳妥。
如果你的编译器不支持 __int128，可以用"快速乘"把乘法拆成加法和位运算来避免溢出。原理很简单：乘法本质上就是多次加法，而取模运算对加法同样保序，所以把乘数按二进制拆开累加，每次加完立刻取模，就能让中间值始终保持在安全范围。这是以时间换空间的写法，但关键时候能救命。
2.3 三个边界：n=0、负数指数、底数为 0
我调试快速幂踩过的坑，几乎都集中在边界条件上。
- n=0：任何非零底数的 0 次方是 1，循环体不执行，result 必须初始化为 1。如果当初把 result 初始化成 0，n=0 时直接返回错误答案。
- 负数指数：整数快速幂的 while 循环面对负数是死循环，因为右移永远不为 0。处理方式取决于场景：实数域直接 1.0 / fast_pow(a, -n)；模质数域则借助费马小定理，先算 a^(p-2) 得到 a 的逆元，再对逆元做快速幂，这样就能计算负指数的模意义结果。
- 底数 a 是 mod 的倍数：此时 a % mod = 0，只要 n > 0 结果必然是 0；但 n = 0 时结果是 1。所以函数入口处最好先判断 n == 0，再做 a %= mod，避免底数取模变成 0 后把幂的结果全部带偏。
另外有一条实战经验：底数很大时，一定要在进循环前先 a %= mod。我曾在某个项目里忘了这行，底数是一个几千位的随机数，前几轮 result 直接溢出，结果怎么都不对，排查了很久才发现 bug 出现在函数第一行，而不是算法逻辑。这种低级错误只要形成习惯就不会再犯。
整数快速幂的坑处理完，下一步就是把同一个模板搬到矩阵上。你会发现代码几乎不需要改动，变的只是"数"变成了"矩阵"。
3. Strassen 矩阵快速幂：从数到矩阵的分治推广
3.1 为什么矩阵可以原样套用快速幂
回顾第一节末尾那句话：快速幂的所有推导只依赖乘法结合律和单位元存在。矩阵乘法满足结合律，(AB)C = A(BC)，而且 n 阶单位矩阵 E 满足 EA = AE = A。所以对矩阵 A 求 A^n，公式和整数完全同构：
A^0 = E A^n = (A^(n/2))² n 为偶数 A^n = (A^(n/2))² · A n 为奇数唯一的新增注意点是矩阵乘法不满足交换律，所以"乘 A"的方向要统一。快速幂实现里我们始终用 res = res * A 这个方向，不要写成 A * res。原理在于：展开过程中底数始终是同一个矩阵 A，结合律保证括号怎么加都对；但一旦交换了顺序，结果就可能完全不一样。这个坑对刚接触矩阵快速幂的人非常容易踩，一定要在代码里固化写乘法时的参数顺序。
3.2 "Strassen"这个称呼：先说清两种算法
标题里的"Strassen"值得单独解释一下，因为中文资料对这两个概念经常混着用。
一方面，1969 年提出的 Strassen 矩阵乘法算法，使用分治策略把 2×2 分块矩阵的乘法从 8 次实数乘法降到 7 次，从而把矩阵乘法复杂度从 O(n³) 优化到约 O(n^2.807)。另一方面，很多教程和题解里的"Strassen 法求矩阵幂"，指的是"对指数做二进制分解、对矩阵做平方"的这套快速幂框架，它本身不涉及 7 次乘法那套优化，只是名字挂在这里。
这两种算法一个针对指数分治，一个针对矩阵乘法本身分治，容易混淆。我的建议是先把基础版矩阵快速幂吃透——内层使用标准 O(k³) 矩阵乘法，整体复杂度 O(k³ log n)，这覆盖了 95% 以上的实际需求。等真遇到 k 特别大（比如上百）且矩阵不太稀疏时，再把内层乘法替换成 Strassen 乘法，把复杂度降到 O(k^2.807 log n)。注意 Strassen 乘法常数非常大，k 小于某个阈值（一般几百）时反而比标准乘法更慢，盲目上 Strassen 属于典型的负优化。
提示：算法竞赛和大多数工程场景里矩阵规模 k 基本不超过 100，标准矩阵乘法完全够用。不要为了"听起来高级"而选择 Strassen 乘法，先把基础模板写对。
3.3 完整实现与实现细节
下面是一个可直接复制的 C++ 矩阵快速幂通用模板：
#include <vector> using namespace std; using Matrix = vector<vector<long long>>; const long long MOD = 1000000007LL; Matrix mat_mul(const Matrix& A, const Matrix& B) { int n = A.size(); int m = B.size(); int p = B[0].size(); Matrix C(n, vector<long long>(p, 0)); for (int i = 0; i < n; ++i) { for (int k = 0; k < m; ++k) { if (A[i][k] == 0) continue; long long v = A[i][k]; for (int j = 0; j < p; ++j) { C[i][j] = (C[i][j] + v * B[k][j]) % MOD; } } } return C; } Matrix mat_pow(Matrix A, long long exp) { int n = A.size(); Matrix res(n, vector<long long>(n, 0)); for (int i = 0; i < n; ++i) res[i][i] = 1; while (exp > 0) { if (exp & 1) res = mat_mul(res, A); A = mat_mul(A, A); exp >>= 1; } return res; }细节挨个说清楚。res 初始化为单位矩阵，对角线全 1，对应整数版本里的 result = 1，这是"什么都没乘"时的单位元。res = mat_mul(res, A) 对应"二进制位为 1 时乘入 A"，A = mat_mul(A, A) 对应"平方"，三个语句的顺序和方向都不能乱。
mat_mul 里我加了 if (A[i][k] == 0) continue，对稀疏矩阵是很好的常数优化。斐波那契这种小矩阵里 0 元素占比高，实测能省两成左右的时间。要注意的是 MOD=10^9+7 时，v * B[k][j] 的结果约 10^18，还在 long long 范围内；如果 mod 更大，这里同样要换成 __int128 中转。
4. 矩阵快速幂的三类典型应用
4.1 斐波那契：最经典的入门案例
二阶线性递推 F(n) = F(n-1) + F(n-2) 可以写成矩阵乘法的形式：
[F(n) ] [1 1] [F(n-1)] [F(n-1)] = [1 0] × [F(n-2)]反复代入右侧，可以得到：
[F(n) ] [1 1]^(n-1) [F(1)] [F(n-1)] = [1 0] [F(0)]于是求第 n 项斐波那契数就变成了求矩阵 M = [[1,1],[1,0]] 的 n-1 次幂，复杂度从 O(n) 直接降到 O(log n)。我第一次看到这个转化时觉得特别惊艳：一个看似只能逐个递推的问题，通过"状态转移矩阵"的抽象，变成了一个可以快速幂的问题。
下面给一份完整可跑的 Python 参考：
def mat_mul(A, B, mod): n, m, p = len(A), len(B), len(B[0]) C = [[0] * p for _ in range(n)] for i in range(n): for k in range(m): if A[i][k] == 0: continue for j in range(p): C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % mod return C def mat_pow(A, exp, mod): n = len(A) res = [[0] * n for _ in range(n)] for i in range(n): res[i][i] = 1 while exp > 0: if exp & 1: res = mat_mul(res, A, mod) A = mat_mul(A, A, mod) exp >>= 1 return res def fib(n, mod=10**9 + 7): if n == 0: return 0 M = [[1, 1], [1, 0]] T = mat_pow(M, n - 1, mod) return T[0][0]注意 fib(0) 必须单独判断：如果 n = 0 时还去算 M^(-1)，快速幂遇到负数指数会直接出问题。这个 bug 我在帮同学看代码时见过不止一次，几乎都是 n 的边界没处理好。
4.2 任意 k 阶线性递推：伴侣矩阵的构造
斐波那契只是特例。任何形如 F(n) = c1·F(n-1) + c2·F(n-2) + ... + ck·F(n-k) 的 k 阶线性齐次递推，都能用同一个模板解决。构造方法叫伴侣矩阵，结构长这样：
[F(n) ] [c1 c2 c3 ... ck] [F(n-1) ] [F(n-1) ] [1 0 0 ... 0 ] [F(n-2) ] [F(n-2) ] = [0 1 0 ... 0 ] [F(n-3) ] [... ] [................] [.... ] [F(n-k+1) ] [0 0 ... 1 0 ] [F(n-k) ]第一行是递推公式本身，下面每一行负责把状态向量往下"移动"一位。构造这个矩阵不需要手抄：第一行按顺序填系数，然后从第二行开始把 (i, i-1) 位置填 1，其余位置填 0，循环一遍就能完成。
举个例子，三阶递推 F(n) = F(n-1) + F(n-2) + F(n-3)，矩阵就是：
[F(n) ] [1 1 1] [F(n-1)] [F(n-1)] = [1 0 0] × [F(n-2)] [F(n-2)] [0 1 0] [F(n-3)]有了这个模板，大量递推题都可以直接套 mat_pow。如果递推式还带常数项，比如 F(n) = c1·F(n-1) + c（非齐次），只要把常数也扩进状态向量，多加一个"永远等于 1"的维度，同样能转成齐次形式。这种"加维度"的思路在状态机模型里特别常见。
4.3 邻接矩阵的幂：图上有多少条 k 步路径
矩阵快速幂另一个反直觉的应用是图论路径计数。设图 G 的邻接矩阵为 A，A[i][j] = 1 表示存在从 i 到 j 的边。那么 A^k 中第 i 行第 j 列的元素，恰好等于从 i 出发走恰好 k 步到达 j 的路径总数。
证明用归纳法很简洁：A² 的 [i][j] 等于 Σ_t A[i][t]·A[t][j]，这正好是所有经过中间点 t 的两步路径数之和；A³ 就是再乘一次 A，把两条路径拼接成三条。所以"求 k 步路径条数"本质上就是"求矩阵 k 次幂"，直接套 mat_pow 即可。
这个性质在工程里的应用比想象中广泛：社交网络中的"k 度关系"挖掘、推荐系统里用户到商品的间接关联强度、知识图谱中两个实体之间的路径依赖度，都可以抽象成邻接矩阵的幂。当然现实场景图很大，矩阵往往是稀疏的，需要配合稀疏矩阵存储和幂迭代算法，这里不展开，但核心思想是完全一致的。
5. 实战选型：什么时候用快速幂，什么时候别硬上
5.1 三种做法的时间复杂度对照
| 求 a^n 的方式 | 乘法次数 | n=10^12 时的大致量级 | 适用场景 |
|---|---|---|---|
| 朴素循环连乘 | n | 10^12 | n 较小，几乎用不上 |
| 整数快速幂 | 约 2·log2(n) | 约 80 | 绝大多数大指数求幂 |
| 矩阵快速幂（k 阶矩阵） | 约 2·log2(n)·k³ | 约 80·k³ | 线性递推、路径计数 |
n = 10^12 时，朴素循环要一万亿次乘法，现代 CPU 也得跑到地老天荒；整数快速幂大概 80 次乘法，瞬间出结果。矩阵快速幂多出来的 O(k³) 来自矩阵乘法本身，这就是为什么遇到递推题要先问一句：k 是多大？
如果 k 是 1000，k³ 就是 10^9，叠加 log n 的 40，整体是 4×10^10 次乘法量级。这种规模硬上矩阵快速幂就不太明智，要么利用矩阵稀疏性做优化，要么换多项式取模等其他方法。
5.2 快速幂救不了的场景
快速幂不是银弹。我总结过三类"硬上反而更亏"的情况。
第一类是指数 n 很小。n 小于几十时，直接循环连乘可能比矩阵快速幂更快，因为矩阵构造和函数调用本身有固定开销，省下的几次乘法不足以回本。第二类是递推式带非线性或卷积项，比如 F(n) = F(n-1) + Σ F(i)·F(n-i)，这种没法写成"常数矩阵乘以状态向量"的形式，需要生成函数、CDQ 分治或 FFT 优化这些更专门的工具。第三类是矩阵规模过大，k 上万时 O(k³ log n) 完全不现实，必须利用循环矩阵、稀疏矩阵、带状矩阵的特殊结构做乘法优化，否则连内存都放不下。
另外补充一个细节：有些递推可以用快速幂但没必要。比如 n 只有 10^5 而且要输出所有中间项时，线性递推一遍扫过去是 O(nk)，矩阵快速幂反而因为要反复构造矩阵而显得笨重。选择算法的顺序永远是先看数据范围，再看能否归约到已知模板，最后才是考虑常数优化。
5.3 两个实测有效的常数优化
最后分享两个我实测有效的优化技巧。
第一个，快速幂主循环里不要把平方写成 if/else 两个分支。很多新手会写成"如果 n 是奇数，先乘再平方；否则只平方"，其实平方无论如何都要做，只需要在 n 为奇数时把当前的 A 乘进 res，然后无条件执行 A = mat_mul(A, A)。写成下面这样最干净：
while (exp > 0) { if (exp & 1) res = mat_mul(res, A); A = mat_mul(A, A); exp >>= 1; }第二个，对于固定的小矩阵，比如 2×2 斐波那契矩阵，手写结构体展开乘法比通用三层循环快一个数量级。通用模板胜在可读性和灵活性，但如果程序里要跑成千上万次矩阵快速幂，固定小矩阵值得手写展开：
struct Mat2 { long long a00, a01, a10, a11; }; Mat2 mul2(const Mat2& x, const Mat2& y, long long mod) { return { (x.a00 * y.a00 + x.a01 * y.a10) % mod, (x.a00 * y.a01 + x.a01 * y.a11) % mod, (x.a10 * y.a00 + x.a11 * y.a10) % mod, (x.a10 * y.a01 + x.a11 * y.a11) % mod, }; }这套东西我用了很多年，最大的体会是：快速幂真正值钱的不是那十几行代码，而是"把指数拆成二进制、把重复连乘换成分治平方"这个思维模型。整数用它，矩阵用它，状态机转移、线性递推、路径计数最后都归一成同一个模板。把这个模板吃透，你会发现自己看许多难题的视角都不一样了——很多看似复杂的模型，本质上不过是一次求幂而已。
