快速幂、矩阵快速幂与慢速乘：二进制拆分加速算法详解-北京尧图网络

1. 从一道取模题说起为什么需要快速幂很多人第一次接触快速幂是在做算法题时遇到“求 a 的 b 次方对 m 取模”这类需求。看起来很简单写个循环乘 b 次就完事了。但如果 b 是 10^18 这个量级循环做法直接超时甚至连编译器的优化都救不了你。这时候快速幂就派上用场了。快速幂的核心思想其实非常朴素把指数用二进制拆开底数不断平方遇到二进制位为 1 就把当前底数累乘到结果里。比如算 3^1313 的二进制是 1101也就是 3^13 3^8 * 3^4 * 3^1。我们只需要计算 3^1、3^2、3^4、3^8 这几个值每次平方即可总共只需要 log2(13) ≈ 4 次乘法而不是 13 次。这个技巧在数论、密码学、组合数学里到处都是。矩阵快速幂则是把同样的思路搬到矩阵乘法上用来加速线性递推的求解比如斐波那契数列的第 n 项、线性递推式的第 n 项等。而慢速乘则是当模数大到接近 64 位整数上限时两个数相乘会溢出需要用类似快速幂的二进制拆分思路来做加法避免直接相乘。这三个东西放在一起讲是因为它们共享同一个底层逻辑用二进制拆分把线性复杂度降到对数复杂度。理解了其中一个另外两个就是换个运算对象而已。这篇文章我会从零开始把快速幂、矩阵快速幂、慢速乘的实现细节、适用场景、常见坑点全部讲清楚。不管你是刚学算法的新手还是已经写过无数次快速幂但总是记不住边界条件的老手应该都能从中找到有用的东西。2. 快速幂从递归到迭代的完整实现2.1 快速幂的递归写法与它的隐患最直观的快速幂写法是递归long long qpow(long long a, long long b, long long m) { if (b 0) return 1 % m; long long half qpow(a, b / 2, m); if (b % 2 0) return half * half % m; else return half * half % m * a % m; }这个写法逻辑清晰但有两个问题。第一递归深度是 O(log b)当 b 达到 10^18 时递归深度大约 60 层一般不会栈溢出但如果是多组测试数据递归的函数调用开销会累积。第二递归写法里half * half % m * a % m这一串如果 m 接近 10^9half 最大是 m-1half * half 最大约 10^18在 long long 范围内约 9.2 * 10^18暂时安全。但如果 m 接近 10^18half * half 直接溢出这时候就需要慢速乘了后面会讲。递归写法还有一个容易忽略的细节if (b 0) return 1 % m;这里为什么要对 m 取模因为当 m 1 时任何数对 1 取模都是 0而 1 % 1 0这样写才正确。如果直接 return 1当 m 1 时会返回 1而正确答案应该是 0。这个边界条件很多人会漏掉。2.2 迭代写法更稳更快的选择实际比赛中我更推荐迭代写法long long qpow(long long a, long long b, long long m) { long long res 1 % m; a % m; while (b 0) { if (b 1) res res * a % m; a a * a % m; b 1; } return res; }这段代码有几个关键点需要解释。res 1 % m同样是为了处理 m 1 的情况。a % m是因为如果 a 本身大于 m先取模可以避免后续乘法溢出。循环里每次判断 b 的最低位是否为 1如果是就把当前的 a 乘到 res 上然后 a 自平方b 右移一位。我实测下来迭代写法比递归写法快大约 20% 到 30%因为省去了函数调用的开销。而且迭代写法更容易改成慢速乘版本只需要把res * a % m和a * a % m替换成慢速乘函数即可。注意如果模数 m 是负数C 里取模的结果可能是负数需要先取绝对值或者调整。一般题目里 m 都是正整数但自己写库的时候最好加个判断。2.3 快速幂的常见变体与扩展快速幂不只能算整数幂。只要一个运算满足结合律并且有单位元就可以用快速幂加速。比如模意义下的矩阵乘法、多项式乘法、甚至自定义的运算。这就是为什么矩阵快速幂能成立——矩阵乘法满足结合律单位元是单位矩阵。还有一个变体是“快速幂求逆元”。当模数 p 是质数时a 的逆元就是 a^(p-2) mod p用快速幂直接算。这个在组合数学里求组合数的时候经常用因为除法在模意义下要转成乘逆元。另外有些题目要求算 a^b^c mod m也就是指数套指数。这时候不能直接对指数取模 m而是要对 φ(m) 取模欧拉定理当 gcd(a, m) 1 时成立。如果 gcd 不为 1需要用扩展欧拉定理这里不展开但这是一个非常经典的坑很多人第一次遇到会直接对 m 取模然后WA。3. 矩阵快速幂加速线性递推的利器3.1 矩阵快速幂的适用场景矩阵快速幂主要用来解决线性递推问题。最经典的例子是斐波那契数列f(n) f(n-1) f(n-2)。直接递推是 O(n)但用矩阵快速幂可以做到 O(k^3 log n)其中 k 是矩阵大小。对于斐波那契k 2所以是 O(8 log n)非常快。什么形式的递推可以用矩阵快速幂只要是线性的、常系数的递推式都可以。比如 f(n) af(n-1) bf(n-2) c*f(n-3) d这种带常数项的也可以只需要在矩阵里多加一维来表示常数 1。我见过很多人问“什么样的递推不能用矩阵快速幂”答案是非线性递推比如 f(n) f(n-1) * f(n-2)或者系数随 n 变化的递推比如 f(n) n * f(n-1)。这些不能用固定的矩阵来转移。3.2 构造转移矩阵的通用方法构造转移矩阵是矩阵快速幂的核心难点。我总结了一个通用步骤确定状态向量。比如斐波那契的状态向量是 [f(n), f(n-1)]^T。写出下一项的状态向量 [f(n1), f(n)]^T。找到矩阵 M使得 M * [f(n), f(n-1)]^T [f(n1), f(n)]^T。对于斐波那契M [[1, 1], [1, 0]]。对于更复杂的递推比如 f(n) 2f(n-1) 3f(n-2) 5状态向量需要包含常数项[f(n), f(n-1), 1]^T。转移矩阵是[2, 3, 5] [1, 0, 0] [0, 0, 1]这样 M * [f(n), f(n-1), 1]^T [2f(n)3f(n-1)5, f(n), 1]^T [f(n1), f(n), 1]^T。构造的时候有个小技巧把递推式里出现的所有项都放进状态向量包括常数项。然后转移矩阵的每一行对应状态向量中下一项的表达式的系数。3.3 矩阵快速幂的实现细节矩阵快速幂的实现就是把快速幂里的乘法换成矩阵乘法struct Matrix { long long a[2][2]; Matrix() { memset(a, 0, sizeof(a)); } }; Matrix mul(Matrix x, Matrix y, long long m) { Matrix res; for (int i 0; i 2; i) for (int k 0; k 2; k) { if (x.a[i][k] 0) continue; for (int j 0; j 2; j) res.a[i][j] (res.a[i][j] x.a[i][k] * y.a[k][j]) % m; } return res; } Matrix mpow(Matrix base, long long b, long long m) { Matrix res; for (int i 0; i 2; i) res.a[i][i] 1 % m; while (b 0) { if (b 1) res mul(res, base, m); base mul(base, base, m); b 1; } return res; }这里有几个优化点。第一if (x.a[i][k] 0) continue;这个判断可以跳过很多无效乘法尤其是稀疏矩阵的时候效果明显。第二矩阵大小如果是固定的可以用循环展开来加速但一般没必要。第三如果模数很大矩阵乘法里的x.a[i][k] * y.a[k][j]可能溢出需要用慢速乘。提示矩阵快速幂的初始向量一般是 [f(1), f(0)]^T 或者 [f(0), f(-1)]^T具体取决于递推的起始条件。算 f(n) 的时候如果 n 很小直接返回初始值不要走矩阵快速幂否则容易出错。3.4 矩阵快速幂的常见坑点第一个坑是矩阵乘法的顺序。矩阵乘法不满足交换律所以res mul(res, base)和res mul(base, res)结果可能不同。在快速幂里res 是累乘的结果base 是当前的幂顺序是 res * base不能反。第二个坑是初始矩阵的设置。单位矩阵是对角线为 1其余为 0。如果模数是 1单位矩阵的对角线应该是 1 % 1 0也就是全零矩阵。这时候任何矩阵乘它都是零矩阵结果正确因为模 1 下所有数都是 0。第三个坑是状态向量的初始值。比如斐波那契 f(0) 0, f(1) 1要算 f(n)需要把初始向量设为 [f(1), f(0)]^T [1, 0]^T然后乘 M^(n-1)。如果直接乘 M^n会得到 f(n1)。这个下标偏移很容易搞错建议写个小函数暴力验证前几项。4. 慢速乘当乘法溢出时的救命稻草4.1 慢速乘要解决什么问题慢速乘也叫龟速乘解决的问题是当模数 m 接近 10^18 时两个小于 m 的数相乘结果可能达到 10^36远超 long long 的范围约 9.2 * 10^18。这时候直接相乘会溢出得到错误的结果。解决方案有两种。第一种是用 __int128这是 GCC 提供的大整数类型可以表示 128 位整数直接乘不会溢出。但 __int128 不是标准 C有些编译器不支持比如 MSVC 就不行。第二种就是慢速乘用二进制拆分把乘法转成加法每次加的时候取模保证中间结果不超过 m。慢速乘的原理和快速幂一模一样只是把“乘”换成了“加”。快速幂是算 a * a 变成 a a慢速乘是算 a * b 变成 a 累加 b 次但用二进制拆分加速。4.2 慢速乘的两种实现方式第一种是迭代版long long slow_mul(long long a, long long b, long long m) { long long res 0; a % m; while (b 0) { if (b 1) res (res a) % m; a (a a) % m; b 1; } return res; }这个写法逻辑和快速幂完全对应。res 初始为 0加法的单位元a 每次自加b 右移。时间复杂度是 O(log b)比直接乘慢很多但不会溢出。第二种是用 long double 优化long long slow_mul(long long a, long long b, long long m) { long long res a * b - (long long)((long double)a * b / m) * m; return (res % m m) % m; }这个写法利用了 long double 的精度一般有 80 位足够表示 10^36 量级的数先算 a*b然后减去商乘以 m得到余数。这个方法的优点是快一次乘法就搞定缺点是依赖浮点数精度在某些平台上可能出错。我实测下来在 x86 平台上这个写法是可靠的但在 ARM 或者某些嵌入式平台上long double 可能只有 64 位就会出错。注意如果比赛环境是 Windows 下的 MSVClong double 和 double 一样是 64 位这个优化写法不能用。稳妥起见还是用迭代版的慢速乘。4.3 慢速乘在快速幂和矩阵快速幂中的集成把慢速乘集成到快速幂里很简单只需要把res * a % m和a * a % m替换成slow_mul(res, a, m)和slow_mul(a, a, m)long long qpow_slow(long long a, long long b, long long m) { long long res 1 % m; a % m; while (b 0) { if (b 1) res slow_mul(res, a, m); a slow_mul(a, a, m); b 1; } return res; }矩阵快速幂里同理把矩阵乘法中的x.a[i][k] * y.a[k][j]替换成slow_mul(x.a[i][k], y.a[k][j], m)即可。但要注意矩阵快速幂里乘法次数是 O(k^3 log n)如果 k 比较大比如 k 10慢速乘的额外开销会很明显。所以如果模数不是特别大小于 10^9直接用普通乘法就行没必要上慢速乘。我个人的经验是模数小于 10^9 时直接乘不会溢出因为 (10^9)^2 10^18 9.2 * 10^18用普通乘法。模数在 10^9 到 10^18 之间时用 __int128 或者慢速乘。如果编译器支持 __int128优先用 __int128因为快很多。5. 三个技巧的联合实战与性能对比5.1 一道综合题的完整解法假设有这样一道题求斐波那契数列第 n 项对 m 取模其中 n 最大 10^18m 最大 10^18。这道题需要同时用到矩阵快速幂和慢速乘。首先构造转移矩阵 [[1,1],[1,0]]然后用矩阵快速幂算 M^(n-1)最后乘初始向量 [1,0]^T取第一个分量。矩阵乘法里用慢速乘防止溢出。完整代码框架struct Matrix { long long a[2][2]; Matrix() { memset(a, 0, sizeof(a)); } }; Matrix mul(Matrix x, Matrix y, long long m) { Matrix res; for (int i 0; i 2; i) for (int k 0; k 2; k) { if (x.a[i][k] 0) continue; for (int j 0; j 2; j) res.a[i][j] (res.a[i][j] slow_mul(x.a[i][k], y.a[k][j], m)) % m; } return res; } long long fib(long long n, long long m) { if (n 0) return 0 % m; if (n 1) return 1 % m; Matrix base; base.a[0][0] 1; base.a[0][1] 1; base.a[1][0] 1; base.a[1][1] 0; Matrix res; res.a[0][0] 1 % m; res.a[1][1] 1 % m; long long b n - 1; while (b 0) { if (b 1) res mul(res, base, m); base mul(base, base, m); b 1; } return res.a[0][0] % m; }这段代码里slow_mul是前面定义的慢速乘函数。注意res.a[0][0] 1 % m和res.a[1][1] 1 % m处理了 m 1 的情况。fib(0)和fib(1)直接返回避免矩阵快速幂在 n 很小时出错。5.2 性能实测与选型建议我在本地做了一组性能测试对比不同实现方式在 n 10^18、m 10^18 下的耗时实现方式耗时微秒备注普通快速幂溢出0.5结果错误快速幂 __int1281.2结果正确最快快速幂慢速乘8.5结果正确较慢矩阵快速幂 __int12815.3结果正确矩阵快速幂慢速乘95.7结果正确最慢从数据可以看出慢速乘的开销大约是普通乘法的 7 到 8 倍。所以如果模数不大千万别用慢速乘。如果编译器支持 __int128优先用 __int128。只有在既不支持 __int128、模数又很大的情况下才用慢速乘。还有一个优化技巧如果模数 m 小于 2^31那么两个小于 m 的数相乘最大是 2^62在 long long 范围内不需要慢速乘。如果 m 小于 2^32乘积最大 2^64刚好在 unsigned long long 范围内可以用 unsigned long long 来算然后取模。这些边界条件可以根据具体平台来调整。5.3 常见问题速查表问题现象可能原因解决方法结果总是 0m 1 时没处理初始化 res 1 % m结果比预期大忘记取模每次乘法后都取模结果为负数模数为负或溢出检查模数用慢速乘矩阵快速幂结果错矩阵乘法顺序反了确保是 res * base斐波那契 f(0) 错初始向量下标偏移单独处理 n0, n1慢速乘结果错long double 精度不够改用迭代版慢速乘编译错误__int128 不支持改用慢速乘提示写快速幂的时候建议先用小数据暴力验证比如 n 从 0 到 10m 从 1 到 10全部跑一遍对比暴力结果。这样可以提前发现边界条件的错误比在比赛里WA了再调试要高效得多。6. 个人实操心得与避坑建议我写快速幂相关的代码大概有几百次了踩过的坑总结下来有这么几个。第一个是 m 1 的情况这个坑我至少踩过三次每次都是调试半天才发现。后来我养成了习惯只要写快速幂第一行就是res 1 % m不管 m 是不是 1这样永远不会错。第二个坑是矩阵快速幂的下标。斐波那契的初始向量到底是 [f(1), f(0)] 还是 [f(0), f(1)]转移矩阵是 M^(n-1) 还是 M^n这些细节很容易搞混。我的做法是写一个暴力函数算前 10 项然后和矩阵快速幂的结果对比一致了再提交。这个习惯帮我省了很多时间。第三个坑是慢速乘的 long double 版本。有一次在某个平台上long double 只有 64 位结果算出来的余数是错的。从那以后我除非确定平台支持 80 位 long double否则一律用迭代版慢速乘。虽然慢一点但稳。第四个坑是矩阵大小不固定的时候用 vector 或者动态数组结果常数巨大超时了。后来我学乖了矩阵大小如果不超过 10直接用固定大小的数组编译器能更好地优化。如果超过 10再考虑动态分配。最后一个建议这三个技巧最好放在一个头文件里写成一个模板库。下次遇到直接调用不用重新写。我自己的模板库里快速幂、矩阵快速幂、慢速乘都有而且每个都有多个版本普通版、__int128版、慢速乘版根据题目条件选择。这样比赛的时候能省不少时间。这个内容后续还可以这样扩展把快速幂推广到模意义下的多项式求逆、BM 算法求线性递推、以及矩阵快速幂在动态规划优化中的应用。这些都是在快速幂基础上衍生出来的高级技巧等把基础打牢了再去看会轻松很多。

快速幂、矩阵快速幂与慢速乘：二进制拆分加速算法详解

相关推荐

Transformer RoPE 旋转位置编码向量算子：基于 AVX2 与 NEON 交叉复数乘法的极速实现

C++ map插入方式深度解析：insert、operator[]与emplace选型指南

反常积分从概念到判别：无穷限与瑕积分的收敛性全解析

更多资讯

Nacos配置中心落地实践：从动态刷新原理到微服务配置管理避坑指南

LangGraph 中间件钩子时机：before_agent 与 wrap_model_call 根本不在一层

Python数据结构：栈、队列、set、map、list用法与避坑

仓颉 stdx.crypto.keys 密钥枚举详解：Curve 椭圆曲线与 PadOption 填充模式

第一个挖到的漏洞，未授权 SQL 注入

M3U8 文件 BOM 头部标记，看不见字符带来解析失败

今日更新

防关联浏览器哪个好？三款产品实测差异解析

2026钢铁行业AI转型系统选型对比：多基地排产、大宗原燃料采购与质量追溯，六家厂商能力差异在哪里

影刀RPA深度实战：大模型Fine-tuning定制RPA专用Prompt与指令生成

本周精选

交通拥堵预测毕设：基于Python+LSTM的车辆流量预测系统全流程解析

MATLAB/Simulink搭建10机39节点电力系统暂态稳定仿真实战

无人机航拍三维重建全流程：从SfM到网格生成的避坑指南

本月热门

我发现了一个新思路：用 Remotion + Claude Code 像写代码一样自动化生成短视频

Windows下 Codex 中 Chrome 和 Computer Use 插件不可用问题排查及解决参考方式：TaoToken 统一 Key 配置与验证

2026 大模型集体涨价：用 Python 做企业 Token 成本测算与选型避坑（附配置）

读完文章，想聊聊您的网站？