1. 从“矩阵快速幂不会”到“力扣模板直接秒”先说一下我为什么要写这题。去年某次集训队友随机抽了一道信奥题让我讲抽到的就是 P3390 模板。当时我第一反应是“这题不是背模板就行了吗”结果真让我手写的时候居然卡在了单位矩阵初始化和乘法顺序上。那种感觉挺尴尬的——平时用现成模板用惯了真要自己从头写反倒不踏实。所以这篇博文不是那种“欢迎来到矩阵快速幂教程”的套路文而是一个踩过坑、翻过车、最后把模板吃透的人把整个过程重新梳理一遍。题目本身只有一个核心任务给定一个 n×n 矩阵 A 和一个很大的指数 kk 最大到 1e9 级别求 A 的 k 次方并对 1e97 取模。看起来简单但它背后牵扯出两个东西——矩阵乘法怎么实现才不溢出、快速幂怎么写才能处理“指数巨大”的情况。如果你正在学信奥或者刷题刷到递推和优化问题这篇文章应该能帮你省下至少一个晚上的懵逼时间。我也会把“为什么这样写”讲清楚而不是丢一份代码让你背了事。2. 整体设计与思路拆解2.1 为什么“暴力连乘”一定不行先说说最简单的思路直接循环 k 次每次做一次 n×n 的矩阵乘法。如果 n2、k4那很快但要是 n100、k1e9那总计算量就是 1e9 × 100^3数字大得离谱。光是乘法次数就爆炸更别说还要取模。所以核心思路只有一个——快速幂。快速幂本身不是一个新东西几乎所有学过算法的人都见过整数的快速幂比如算 a^13不是乘 13 次而是把 13 拆成二进制1101对应 a^8 × a^4 × a^1乘三次加上一次幂运算就够了。矩阵快速幂的思想一模一样只不过底数从标量 a 换成了矩阵 A“乘法”换成了矩阵乘法。这里我画一个思维链条方便你理解整数快速幂能高效计算 a^k是因为满足结合律 (a^m)^n a^(mn)。矩阵乘法也满足结合律所以同样可以把指数拆成二进制用“倍增”的方式来算。2.2 指数二进制拆解的实际含义举个例子算 A^13。13 的二进制是 1101也就是 13 8 4 1。那 A^13 A^8 × A^4 × A^1。怎么得到 A^8、A^4、A^1很简单从 A^1 开始不断自乘A^2 A×AA^4 A^2 × A^2A^8 A^4 × A^4。每轮只需要一次矩阵乘法。整个过程的复杂度是 O(log k) 次矩阵乘法每次矩阵乘法是 O(n^3)。总复杂度 O(n^3 log k)在 n100、k1e9 的约束下稳稳通过。这还没完。如果只是“能算”还不够还要考虑“怎么把代码写得不容易出错”。矩阵快速幂的模板代码看起来短真正写对的人不多原因就在于它对初始化和乘法细节非常敏感。2.3 为什么我决定用“结构体 全局 n __int128”这个组合我翻过很多题解发现大家写法五花八门有人用类封装有人用动态二维数组有人用 vector有人直接在 main 里写三重循环。这些都能过但作为模板我希望代码具备下面几个特质可读性强变量名清楚逻辑分层初学者对着看能懂。不容易出错尽量少依赖“容易写错”的写法。性能足够在 n100 时不会因为常数太大被卡。所以我基本采用“结构体 固定数组”的方案理由也很直白结构体可以直接返回函数签名干净固定数组不需要动态分配避免内存泄漏和越界问题全局 n 配合结构体里的二维数组写起来流畅、跑得快。至于__int128这个属于“想一步到位就别偷懒”的选择。后面我会专门讲这里先记住一点在取模类题目中用long long直接做乘法累加是存在溢出风险的用__int128可以一次性累加完全部乘积再取模既干净又不容易踩坑。3. 核心细节解析与实操要点3.1 单位矩阵整个模板最容易出错的第一步很多人第一次写矩阵快速幂都是照着整数快速幂的代码克隆的。整数快速幂里res 初始化为 1矩阵快速幂里res 应该初始化为“矩阵意义上的 1”也就是单位矩阵——对角线为 1、其余为 0 的矩阵。为什么必须是单位矩阵因为单位矩阵 E 满足 E×A A×E A相当于整数里的 1。如果你把 res 初始化为零矩阵那不管怎么乘结果都还是零矩阵。这个错误非常隐蔽因为 k1 时正确答案应该是原矩阵而零矩阵直接输出全 0没跑几个样例根本发现不了。下面是我的init_identity写法void init_identity() { memset(a, 0, sizeof(a)); for (int i 0; i n; i) a[i][i] 1; }注意两点第一先 memset 清零避免上次使用残留。第二对角线赋值 1 时用的是全局 n而不是数组固定大小。这个细节看起来很无聊但如果是结构体内部的构造函数里写for (int i 0; i n; i)而 n 是全局变量在构造时可能还没读入就会产生未定义行为。这类问题非常难排查。3.2 矩阵乘法循环顺序与累加时机矩阵乘法的定义是C[i][j] A[i][0]*B[0][j] A[i][1]*B[1][j] … A[i][n-1]*B[n-1][j]。代码上最常见的写法是三重循环for (int i 0; i n; i) for (int j 0; j n; j) for (int t 0; t n; t) C.a[i][j] A.a[i][t] * B.a[t][j];这个写法在 n 小、数值不大时没问题。但注意两点数值累加可能溢出A[i][t] 和 B[t][j] 都是 1e97 级别的数两者相乘是 1e18 级别long long 上限大约是 9.2e18。看起来不超但如果 n100累加 100 次就是 1e20绝对爆掉。所以要么边加边取模要么用 __int128 累加。循环顺序影响性能i、j、t 的顺序不同访存局部性完全不同。对于固定数组推荐 i、j、t 外层循环因为内层 t 访问时 A 的行连续、B 的列连续但 B 是按行存储的所以 B[t][j] 实际上是非连续访问。不过 n100 时差距不会天翻地覆习惯上我依旧推荐最直观的 i、j、t。我用 __int128 累加完整行乘积后再取模代码长这样Matrix mul(const Matrix A, const Matrix B) { Matrix C; for (int i 0; i n; i) for (int j 0; j n; j) { __int128 sum 0; for (int t 0; t n; t) sum (__int128)A.a[i][t] * B.a[t][j]; C.a[i][j] (long long)(sum % MOD); } return C; }这里一次乘法运算就把所有加数累加完毕避免了多次取模的开销也彻底避免了 long long 溢出。__int128 在绝大多数现代评测机上都能用只有极老的编译器可能不支持遇到那种情况再用边乘边取模方案。3.3 快速幂主体和整数快速幂一模一样矩阵快速幂的主体基本就是整数快速幂的翻版只是乘法换成矩阵乘法初始值换成单位矩阵。Matrix qpow(Matrix base, long long exp) { Matrix res; res.init_identity(); while (exp) { if (exp 1) res mul(res, base); base mul(base, base); exp 1; } return res; }如果 exp 二进制当前位是 1就把当前的 base 乘到结果里无论当前位是不是 1base 都要自乘因为下一位需要的是 base 的平方。这个逻辑背后对应的是“拆二进制指数”。当我们说 13841实际上是“遇到二进制位为 1 的那一轮把对应幂次乘进去”。值得强调一个细节res mul(res, base)不能写成res mul(base, res)。虽然矩阵乘法不满足交换律写反会导致答案错误。我第一次写就是因为图省事把顺序写成 base 在前结果样例一个都没过。这个坑在里面很常见因为它不像整数乘法那样“怎么乘都行”。4. 实操过程与核心环节实现4.1 实验环境的准备这个题目在绝大多数 C 环境下都能跑不需要额外依赖。我习惯用 Linux 环境下 g 11 编译命令是g -O2 -stdc17 main.cpp -o main-O2 优化对于矩阵运算的循环展开很有帮助信奥和力扣的评测环境通常也会开 O2。如果你在 Windows 上用某些老 IDE优先确认编译器版本支持 C11 以上否则 __int128 可能不兼容。测试直接用题目自带的样例以及自己额外构造的边界数据。比如 n1 的情况、k0 的情况、k1 的情况。n1 时矩阵就一个数快速幂退化为普通幂运算特别适合验证结果。4.2 从头推导一个 n2 的手算样例假设A [[2, 3], [0, 1]]A^2 [[2, 3], [0, 1]] × [[2, 3], [0, 1]] [[4, 9], [0, 1]]A^4 [[4, 9], [0, 1]] × [[4, 9], [0, 1]] [[16, 45], [0, 1]]那 A^5 A^4 × A^1 [[16, 45], [0, 1]] × [[2, 3], [0, 1]] [[32, 93], [0, 1]]手算完之后再跑代码如果输出一致基本可以确定矩阵乘法和幂运算这两个核心模块没问题。这种“先用小规模手算再对拍”的方法比直接丢大样例省心太多。4.3 完整代码下面就是我在实战中一直沿用的模板。注释写得比较细方便你对照上面的讲解看。#include bits/stdc.h using namespace std; typedef long long ll; const ll MOD 1000000007LL; int n; ll k; struct Matrix { ll a[105][105]; Matrix() { memset(a, 0, sizeof(a)); } void init_identity() { memset(a, 0, sizeof(a)); for (int i 0; i n; i) a[i][i] 1; } }; Matrix mul(const Matrix A, const Matrix B) { Matrix C; for (int i 0; i n; i) for (int j 0; j n; j) { __int128 sum 0; for (int t 0; t n; t) sum (__int128)A.a[i][t] * B.a[t][j]; C.a[i][j] (ll)(sum % MOD); } return C; } Matrix qpow(Matrix base, ll exp) { Matrix res; res.init_identity(); while (exp) { if (exp 1) res mul(res, base); base mul(base, base); exp 1; } return res; } int main() { ios::sync_with_stdio(false); cin.tie(nullptr); cin n k; Matrix A; for (int i 0; i n; i) for (int j 0; j n; j) cin A.a[i][j]; Matrix ans qpow(A, k); for (int i 0; i n; i) { for (int j 0; j n; j) cout ans.a[i][j] % MOD (j n - 1 ? \n : ); } return 0; }说两个容易忽略的地方。一个是Matrix()构造函数会清零整个数组所以mul里直接定义Matrix C;再累加即可。另一个是输出格式用三目运算符控制行末换行避免多打空格被评测系统判 Presentation Error。虽然在信奥里通常 PE 也算错但能少一个算一个。5. 常见问题与排查技巧实录5.1 我踩过的三个真实大坑坑一res 初始化成零矩阵这个我在 3.1 里提过。第一次提交时样例过了但 k2 的自测没做交上去全 WA。后来发现是 res 初始化不对。从那以后我形成肌肉记忆写快速幂先写单位矩阵初始化再写循环。坑二乘法顺序写反矩阵乘法不满足交换律AB ≠ BA大多数情况下。我有一段时间为了省事把mul(res, base)写成mul(base, res)结果小样例就错。后来想明白一个道理乘法顺序不是拍脑袋定的它取决于你的递推方向。如果你定义 res 是“已经累积的结果”那 res 在左、base 在右因为 res res * base 才是把当前幂次追加到结果后面。这个顺序一旦定下来整个代码都要保持一致。坑三用 long long 直接累加导致溢出我最初写的乘法是C.a[i][j] (C.a[i][j] A.a[i][t] * B.a[t][j]) % MOD;小样例全对一上大样例就错。后来加打印才发现中间某个 sum 已经超过 long long 上限。原因很简单A[i][t] * B[t][j] 是 1e18 级别再加一个也是 1e18 级别两个加起来就爆了。用 __int128 是治本的办法。5.2 常见问题速查表问题现象可能原因解决方式k0 时输出全 0res 初始化为零矩阵而非单位矩阵把 res 初始化为单位矩阵k1 时输出不对矩阵乘法函数写错或输入顺序不对先单独测乘法手算一个 2×2 的例子大样例答案错误long long 累加溢出使用 __int128 或边乘边取模输出多空格行末多输出空格用 if/三目判断是否为行末编译不通过编译器不支持 __int128退回 long long在循环内及时取模数组越界n 在结构体构造时还未读入用全局 n 或把 n 作为参数传递5.3 如何快速定位是“乘法写错”还是“快速幂写错”如果你发现答案不对不要一上来就通读代码。先设计三个小测试k1输出应该等于原矩阵。如果不对问题不在快速幂而在乘法或者读入。k2用 2×2 矩阵手算核对平方结果。k0输出应该是单位矩阵。如果不对初始化写错了。这三个测试通过后再用题目样例基本能直接定位到出错模块。这个方法我用了很多年几乎每次都能在五分钟内找到问题。5.4 关于取模时机的进一步说明很多人会问sum % MOD应该在哪一步做方案一每次累加都取模C.a[i][j] (C.a[i][j] A.a[i][t] * B.a[t][j]) % MOD;方案二累加完再统一取模__int128 sum 0; for (...) sum (__int128)A.a[i][t] * B.a[t][j]; C.a[i][j] sum % MOD;我推荐方案二理由是代码更干净、性能更好。前提是你用了 __int128 防止累加溢出。如果编译器不支持 __int128那只能用方案一但也需要先把乘法结果转成 long long 后再加且在累加前对结果取模。6. 从这道题延伸出去的三个学习方向6.1 线性递推加速斐波那契数列矩阵快速幂最经典的应用之一是斐波那契数列。斐波那契递推式 f(n) f(n-1) f(n-2) 可以被写成矩阵形式[f(n), f(n-1)] [f(n-1), f(n-2)] × [[1, 1], [1, 0]]进一步可以写成 [f(n), f(n-1)] [1, 0] × [[1, 1], [1, 0]]^(n-1)。于是求第 n 项就变成了矩阵快速幂复杂度从 O(n) 降到 O(log n)。很多题目看着和矩阵无关但一旦发现状态转移是“线性”的就能套这个模板。6.2 图论中的应用邻接矩阵的幂另一类常见套壳是图论。邻接矩阵 G 表示两点之间是否有边相连G^k[i][j] 表示从 i 到 j 是否存在长度为 k 的路径。因为矩阵乘法本质是在做“路径拼接”经过一个中间节点把路径数累积起来。这种题在数数类问题中出现频率很高尤其和取模结合时正好用上矩阵快速幂模板。6.3 整数快速幂和矩阵快速幂的统一认识学到这里你会发现整数快速幂和矩阵快速幂只是“底数乘法”不同。如果把“支持乘法和单位元”的东西抽象成一个概念那广义快速幂还能应用到多项式、模意义下的数、甚至一些特殊定义的操作上。这也是为什么很多高手说“模板不重要理解快速幂的本质才重要”。不过对刚入门的人来说先把矩阵版本写对、写熟比抽象理解更有价值。等你写熟练了再回头看整数和矩阵的共性感会有一种“原来如此”的通透。7. 一些个人经验和最后的小建议这题作为模板题最大的价值不是“背下来”而是让你在一次又一次的“为什么”里把底层逻辑问透。我见过不少人把矩阵快速幂模板抄了十遍但遇到变体题依然不会做原因就是没搞清楚单位矩阵、乘法顺序、结合律这些细节。如果让我给一个学习路径我会推荐先彻底理解整数快速幂再手写矩阵乘法并测试然后自己实现矩阵快速幂最后做几道递推加速题和路径计数题。整个过程大概两三天但收获比背十遍模板都大。调试方面我最后分享一个很土但非常管用的技巧在乘法函数里临时加一个assert或者打印中间矩阵观察 k1、k2 时的输出。一旦发现结果和你手算的不一致立刻能看出是第几行矩阵错了。这个方法比盯着代码空想高效得多。关于取模常量我还建议养成用constexpr long long MOD 1000000007LL;的习惯因为这种常量在编译期就确定且不容易被误改。至于 __int128 的兼容性现代 Linux 下的 g 和主流 OJ 基本都支持真遇到不支持的平台也可以换成“每加一次就取模”的写法只是代码会稍微繁琐一点。希望这篇梳理能帮你少踩几个坑。如果你在写矩阵快速幂时还有什么奇怪的 bug欢迎按上面速查表自己排查一轮大概率能解决大部分问题。