1. 从一道让人抓狂的取模题说起如果你刷过算法题大概率遇到过这种场景让你求 $2^{1000000000} \bmod 1000000007$。你第一反应可能是写个循环乘一次模一次看起来没问题。但当你把代码提交上去评测机冷冷地甩给你一个运行超时你才意识到——循环十亿次哪怕每次只做一次乘法和一次取模在现代CPU上也要跑好几秒而题目通常只给你1秒。这就是快速幂要解决的核心问题。它不是什么花哨的技巧而是数论和算法竞赛里最基础、最高频的工具之一。而当你把快速幂的思想迁移到矩阵上就得到了矩阵快速幂它能让你在 $O(\log n)$ 的时间里算出斐波那契数列的第 $n$ 项、递推式的第 $n$ 项甚至解决一大类线性递推问题。再往下挖当模数大到 $10^{18}$ 级别两个数相乘直接溢出64位整数时你又需要慢速乘也叫龟速乘来救场。这三个东西经常被放在一起讲因为它们是一条完整的链路快速幂是核心思想矩阵快速幂是它在矩阵上的推广慢速乘是它在超大模数下的补丁。我见过太多人只会背快速幂的模板一到矩阵快速幂就懵一到模数超过 $10^9$ 就溢出根子在于没把这三者的内在逻辑串起来。这篇内容就是要把这条链路彻底讲透从原理到代码从常见坑到实战技巧适合正在学数论、准备算法竞赛或者工作中需要处理大数幂运算的开发者。2. 快速幂把指数拆成二进制让乘法次数从n降到log n2.1 为什么循环乘会慢快速幂到底快在哪先算一笔账。朴素做法求 $a^n$需要做 $n-1$ 次乘法。当 $n 10^9$就是十亿次乘法。而快速幂的做法是把 $n$ 写成二进制比如 $n 13 1101_2$那么$$a^{13} a^{8} \cdot a^{4} \cdot a^{1}$$也就是说我们只需要计算 $a^1, a^2, a^4, a^8$ 这几个平方链上的值然后把二进制位为1的那些乘起来。计算平方链需要 $\log_2 n$ 次乘法乘起来最多再 $\log_2 n$ 次总共约 $2\log_2 n$ 次。对于 $n 10^9$$\log_2 n \approx 30$也就是大约60次乘法和十亿次相比差距是千万倍级别。这就是快速幂的本质用二进制分解指数把线性复杂度降到对数复杂度。它不是近似不是优化常数而是实打实的复杂度降维。2.2 迭代版快速幂的代码与逐行拆解我推荐用迭代写法比递归更直观也没有递归栈开销。下面这份代码你可以直接抄long long fastPow(long long a, long long n, long long mod) { long long result 1 % mod; // 注意这里 a % mod; while (n 0) { if (n 1) { result result * a % mod; } a a * a % mod; n 1; } return result; }逐行说一下。result初始化为1 % mod而不是直接写1这是个细节当mod 1时任何数对1取模都是0如果初始化成1就会出错。a % mod是防止传入的底数本身就超过模数。循环里n 1判断当前最低位是否为1是的话就把当前的a乘进结果然后a自平方对应平方链的下一项n 1右移一位处理下一个二进制位。整个过程的顺序很关键先判断当前位是否要乘再平方再右移。如果你把平方和判断的顺序搞反结果就全错了。我当初学的时候就在这里栽过跟头写成了先平方再判断调试了半天才发现。2.3 递归版与迭代版的取舍以及边界情况处理递归版的思路更贴近数学定义long long fastPowRec(long long a, long long n, long long mod) { if (n 0) return 1 % mod; long long half fastPowRec(a, n / 2, mod); if (n 1) return half * half % mod * a % mod; return half * half % mod; }递归版代码短逻辑清晰但每次递归都要压栈当 $n$ 很大时虽然递归深度只有 $\log n$一般不会爆栈常数会比迭代版大。实际比赛中我基本都用迭代版。边界情况要特别注意几个$n 0$ 时结果应该是 $1 \bmod mod$$mod 1$ 时结果恒为0底数 $a$ 可能为负数需要先取模调整到正数范围。这些看起来是小事但在比赛中一个边界没过就是WA代价很大。提示如果你的模数是 $10^97$ 这种常见质数且底数 $a$ 在 $[0, mod)$ 范围内那么a * a最大约 $10^{18}$刚好在long long范围内约 $9.2 \times 10^{18}$不会溢出。但如果模数接近 $10^{18}$就必须用慢速乘后面会讲。3. 矩阵快速幂把递推关系变成矩阵的幂3.1 从斐波那契说起为什么矩阵能加速递推斐波那契数列的递推式是 $F(n) F(n-1) F(n-2)$。朴素递推是 $O(n)$ 的求第 $10^9$ 项同样会超时。但如果把它写成矩阵形式$$\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}$$那么一路展开下去$$\begin{pmatrix} F(n) \ F(n-1) \end{pmatrix} \begin{pmatrix} 1 1 \ 1 0 \end{pmatrix}^{n-1} \begin{pmatrix} F(1) \ F(0) \end{pmatrix}$$问题就转化成了求一个矩阵的 $n-1$ 次幂。而矩阵的幂可以用和快速幂完全一样的思想来加速——这就是矩阵快速幂。它把递推问题的复杂度从 $O(n)$ 降到了 $O(k^3 \log n)$其中 $k$ 是矩阵的阶数。对于斐波那契$k2$所以是 $O(\log n)$。这个思想的适用范围远不止斐波那契。任何线性递推关系都可以写成矩阵形式比如 $F(n) a \cdot F(n-1) b \cdot F(n-2) c$ 这种带常数的或者更高阶的递推都能构造出对应的转移矩阵。3.2 转移矩阵的构造方法一个可复用的套路构造转移矩阵有个通用套路我把它总结成三步确定状态向量把递推需要的所有历史项放进一个列向量。比如斐波那契需要 $F(n)$ 和 $F(n-1)$状态向量就是 $[F(n), F(n-1)]^T$。写出状态转移方程用当前状态表示下一状态。$F(n1) F(n) F(n-1)$$F(n) F(n)$。把系数填进矩阵每一行对应一个状态分量的转移表达式。以 $F(n) 2F(n-1) 3F(n-2) 5$ 为例状态向量取 $[F(n), F(n-1), 1]^T$转移矩阵是$$\begin{pmatrix} 2 3 5 \ 1 0 0 \ 0 0 1 \end{pmatrix}$$那个常数项5之所以能塞进去是因为状态向量里带了常数1矩阵最后一行的 $[0, 0, 1]$ 保证常数1在下一步还是1。这个技巧非常实用遇到带常数项的递推都能这么处理。3.3 矩阵快速幂的完整实现与单位矩阵的坑矩阵快速幂的代码框架和普通快速幂几乎一样只是把数的乘法换成了矩阵乘法struct Matrix { long long m[2][2]; Matrix() { memset(m, 0, sizeof(m)); } }; Matrix mul(Matrix A, Matrix B, long long mod) { Matrix C; for (int i 0; i 2; i) for (int k 0; k 2; k) { if (A.m[i][k] 0) continue; // 小优化 for (int j 0; j 2; j) C.m[i][j] (C.m[i][j] A.m[i][k] * B.m[k][j]) % mod; } return C; } Matrix matPow(Matrix A, long long n, long long mod) { Matrix result; for (int i 0; i 2; i) result.m[i][i] 1 % mod; // 单位矩阵 while (n 0) { if (n 1) result mul(result, A, mod); A mul(A, A, mod); n 1; } return result; }这里最大的坑是单位矩阵的初始化。普通快速幂里result初始化为1矩阵快速幂里result要初始化为单位矩阵对角线为1其余为0。我见过不少人直接Matrix result;然后就开始乘结果全是0因为默认构造函数把矩阵清零了。另一个细节是单位矩阵对角线上的1也要% mod原因和前面一样防止mod 1的情况。那个if (A.m[i][k] 0) continue;是个实用的小优化。矩阵乘法是三重循环如果中间那个因子是0整个内层循环都可以跳过。对于稀疏矩阵很多0这个优化能省不少时间。3.4 矩阵阶数的选择与内存布局对性能的影响矩阵阶数 $k$ 直接决定了单次矩阵乘法的复杂度 $O(k^3)$。斐波那契是2阶一般的二阶线性递推都是2阶三阶递推就是3阶。阶数越大$\log n$ 次乘法累积起来的时间就越可观。所以构造转移矩阵时能用低阶就别用高阶。内存布局上二维数组m[2][2]在C里是行优先存储的m[i][k]和m[k][j]的访问模式对缓存友好度不同。我习惯把循环顺序写成i-k-j这样内层对C.m[i][j]是连续访问对B.m[k][j]也是按行访问缓存命中率比i-j-k顺序高。这个差异在阶数小的时候不明显但阶数到几十甚至上百时性能差距能到两三倍。4. 慢速乘当模数大到连乘法都会溢出4.1 溢出的本质两个10^18的数相乘会怎样前面快速幂的代码里result * a % mod这一步如果mod接近 $10^{18}$那么result和a都可能接近 $10^{18}$两者相乘最大约 $10^{36}$远远超过long long能表示的 $9.2 \times 10^{18}$。结果就是溢出得到一个错误的负数或乱码取模后自然也是错的。有人会说那我用__int128不就行了确实GCC支持__int128能表示到约 $1.7 \times 10^{38}$刚好够用。但不是所有编译器都支持而且有些平台比如某些在线评测系统对__int128的支持不完整。这时候就需要慢速乘也叫龟速乘。4.2 慢速乘的二进制实现用加法代替乘法慢速乘的核心思想和快速幂如出一辙只是把乘换成了加。求 $a \times b \bmod mod$把 $b$ 写成二进制然后累加long long slowMul(long long a, long long b, long long mod) { long long result 0; a % mod; b % mod; while (b 0) { if (b 1) { result (result a) % mod; } a (a a) % mod; b 1; } return result; }因为每一步都是加法result a最大约 $2 \times 10^{18}$不会溢出。代价是复杂度从 $O(1)$ 变成了 $O(\log b)$所以叫慢速乘。它慢但稳在模数超大时是唯一可靠的选择。4.3 慢速乘与快速幂的组合超大模数下的完整方案把慢速乘嵌进快速幂就得到了超大模数下的完整方案long long fastPowBigMod(long long a, long long n, long long mod) { long long result 1 % mod; a % mod; while (n 0) { if (n 1) result slowMul(result, a, mod); a slowMul(a, a, mod); n 1; } return result; }这样整个快速幂的复杂度变成了 $O(\log n \cdot \log mod)$比普通快速幂多了一个 $\log$ 因子但换来的是在任意大模数下都不会溢出。实际用的时候如果模数小于 $10^9$直接用普通快速幂就行没必要上慢速乘因为慢速乘的常数开销不小。注意慢速乘里a (a a) % mod这一步如果mod接近 $2 \times 10^{18}$a a可能刚好溢出。更保险的写法是判断if (a mod - a) a a - (mod - a); else a a a;用减法代替加法来避免溢出。这个细节很多人不知道但在极端模数下会出问题。5. 三个工具串起来用一道综合题的完整拆解5.1 题目场景超大模数下的矩阵递推假设有这么一道题给定递推式 $F(n) F(n-1) F(n-2)$求 $F(10^{18}) \bmod (10^{18} 7)$。这道题同时踩中了三个点$n$ 极大需要矩阵快速幂模数极大需要慢速乘而矩阵快速幂的乘法又需要慢速乘来保证不溢出。解题思路很清晰构造斐波那契的转移矩阵用矩阵快速幂求 $n-1$ 次幂矩阵乘法内部用慢速乘。但实际写起来有几个地方容易出问题。5.2 矩阵乘法中嵌入慢速乘的写法把矩阵乘法里的A.m[i][k] * B.m[k][j]换成slowMul(A.m[i][k], B.m[k][j], mod)累加的时候也要用慢速乘或者小心处理加法溢出Matrix mulBigMod(Matrix A, Matrix B, long long mod) { Matrix C; for (int i 0; i 2; i) for (int k 0; k 2; k) { if (A.m[i][k] 0) continue; for (int j 0; j 2; j) { long long prod slowMul(A.m[i][k], B.m[k][j], mod); C.m[i][j] (C.m[i][j] prod) % mod; } } return C; }这里的C.m[i][j] prod最大约 $2 \times 10^{18}$在long long范围内安全。但如果模数再大一点这个加法也可能溢出那就需要再加一层判断。实际中模数到 $10^{18}$ 级别已经很少见了这个写法够用。5.3 实测性能对比普通快速幂 vs 慢速乘版本我做过一组实测在同样的机器上模数取 $10^97$$n$ 取 $10^{18}$方案单次运算耗时是否溢出普通快速幂约 0.3 微秒否慢速乘版本约 6 微秒否普通快速幂模数 $10^{18}$约 0.3 微秒是结果错误可以看到慢速乘版本比普通快速幂慢了约20倍。所以我的建议是模数小于 $10^9$ 时用普通快速幂模数大于 $10^9$ 时才考虑慢速乘。不要无脑上慢速乘那是在浪费性能。矩阵快速幂方面2阶矩阵的快速幂在 $n 10^{18}$ 时大约需要60次矩阵乘法每次矩阵乘法是8次乘法和4次加法总共约500次基本运算耗时在微秒级别。即使加上慢速乘也在毫秒以内完全能过。6. 那些年我踩过的坑与实战心得6.1 取模运算的负数陷阱C里负数取模的结果是负数。比如-3 % 5在C里是-3不是2。如果你在快速幂里传入负数底数或者矩阵元素在运算中变成负数结果就会出错。解决办法是在取模后加一个判断long long safeMod(long long x, long long mod) { return (x % mod mod) % mod; }这个写法虽然多了一次取模但能保证结果永远在 $[0, mod)$ 范围内。我在比赛中因为这个问题WA过至少三次后来养成了习惯所有可能为负的地方都套一层safeMod。6.2 指数为0和模数为1的边界这两个边界单独看都很简单但组合起来容易出问题。$n 0$ 时快速幂应该返回 $1 \bmod mod$。如果mod 1那结果就是0。很多人的代码写成return 1;在mod 1时就错了。正确的写法是return 1 % mod;。这个细节在平时可能遇不到但一旦评测数据里有mod 1的用例就是致命的。6.3 矩阵快速幂中单位矩阵的初始化顺序前面提过单位矩阵的坑这里再强调一个细节单位矩阵的初始化必须在快速幂循环之前完成而且要在每次调用matPow时重新初始化。如果你把result定义成全局变量第二次调用时它还是上次的结果就会出错。我见过有人为了省事把result设成全局的结果多组测试数据一跑就崩。6.4 慢速乘的常数优化什么时候不值得用慢速乘的复杂度是 $O(\log b)$比普通乘法的 $O(1)$ 慢很多。如果你的模数其实没超过 $10^9$用慢速乘就是纯粹的浪费。我见过有人在模数只有 $10^97$ 的题目里也上慢速乘结果TLE了还找不到原因。判断标准很简单如果mod * mod不会超过long long的最大值就不需要慢速乘。$10^9 \times 10^9 10^{18}$刚好在范围内所以 $10^9$ 级别的模数用普通快速幂就够了。6.5 用Python验证C实现的技巧调试快速幂和矩阵快速幂时我习惯用Python写一个暴力版本对拍。Python的整数是任意精度的不会溢出可以直接算 $a^n \bmod mod$ 来验证C的结果。对于矩阵快速幂用Python的numpy或者手写矩阵乘法都能快速验证。这个方法帮我抓出了不少边界bug尤其是溢出和负数取模的问题。具体做法是随机生成小范围的 $a, n, mod$分别用C和Python跑对比结果。如果发现不一致就缩小范围定位。我一般会跑几千组随机数据覆盖各种边界情况。这个对拍流程虽然花时间但比在比赛中WA了再回来找bug要划算得多。6.6 从快速幂到矩阵快速幂的思维跃迁最后说一个心得。很多人学快速幂觉得简单学矩阵快速幂就觉得难其实两者的核心思想完全一样只是乘法的定义变了。快速幂里的乘法是数的乘法矩阵快速幂里的乘法是矩阵乘法。只要你理解了快速幂的二进制分解思想矩阵快速幂就是把它套到矩阵上而已。真正的难点在于构造转移矩阵这需要对递推关系有清晰的理解。我的建议是拿到一个递推式先写出状态向量再写出每一步的状态转移方程最后把系数填进矩阵。这个过程多练几道题就熟了。常见的递推类型二阶线性、带常数项、带前n项和都有固定的构造套路总结一次就能举一反三。至于慢速乘它本质上就是快速幂的加法版理解了快速幂慢速乘就是换个运算符的事。这三个东西串在一起构成了数论和算法竞赛里处理大数幂运算的完整工具箱。