矩阵快速幂详解：从递推式到转移矩阵的构造与模板

这周的专题训练轮到我讲矩阵快速幂。题单放出去的时候几个刚入队的新人一脸茫然——快速幂不是早就学过了吗整数取模的快速幂他们确实会写但一听到矩阵两个字立马觉得是什么高深玄学。其实矩阵快速幂一点都不玄它就是快速幂的底数从整数换成了矩阵仅此而已。这篇文章就按我给集训队新人拆解的顺序来写先搞清楚它解决什么问题再看标准模板长什么样最后把转移矩阵的构造练到条件反射。无论你是刚接触算法竞赛的初学者还是已经会背板子但一遇到变式就卡壳的老手顺着这条思路推一遍收获会比单纯刷题大得多。1. 从一个递推式说起暴力循环为什么会被卡死1.1 斐波那契一次循环的重复劳动先看最经典的递推式F(n) F(n-1) F(n-2)边界 F(0)0, F(1)1。如果让你求第 n 项大多数人第一反应是写个循环long long fib(int n) { if (n 1) return n; long long a 0, b 1, c; for (int i 2; i n; i) { c a b; a b; b c; } return b; }这个循环的复杂度是 O(n)。n 如果只有 10^6随便跑但题目经常给 n 10^18这时候循环次数是十亿亿次量级。我给学生算过一笔账假设一秒钟跑 3×10^9 次简单运算10^18 次也要跑大约十年。显然不可行。所以问题不是能不能写出来而是怎么把复杂度降下去。仔细看这个循环你会发现每一轮干的事完全相同把旧的 {F(n-1), F(n-2)} 变成新的 {F(n), F(n-1)}。这是一个固定的线性变换。既然操作固定就有加速空间——方法就是矩阵快速幂。1.2 矩阵表示把状态和转移拆开把递推的中间状态看成一个二维向量v(n) [F(n), F(n-1)]^T那么 v(n) 可以由 v(n-1) 通过一个矩阵乘出来[ F(n) ] [ 1 1 ] [ F(n-1) ] [ F(n-1) ] [ 1 0 ] × [ F(n-2) ]验证一下第一行 F(n-1) F(n-2) F(n)第二行 F(n-1)。完全正确。令 M [[1, 1], [1, 0]]于是v(2) M × v(1) v(3) M × v(2) M^2 × v(1) ... v(n) M^(n-1) × v(1)到这里问题从循环 n 次变成了计算矩阵的 n-1 次幂。矩阵乘法本身还是 O(维数^3)但幂次可以用快速幂做到 O(log n)。二维矩阵的维数是 2一共只要算大约 60 次乘法这个加速效果是降维打击式的。关键是理解这个思路任何每一步都在做同一个线性变换的递推都可以写成矩阵形式然后用矩阵快速幂加速。这是整个专题的核心思想后面所有例子都围绕这句话展开。2. 快速幂的底层逻辑二进制拆解与幂次合并2.1 从整数版快速幂说起要理解矩阵快速幂先把整数版的快速幂彻底吃透。整数的快速幂代码很短很多人背下来了但说不上来为什么对。目标是求 a^b mod p。朴素写法是乘 b 次但 b 用二进制表示后可以把乘法次数压到 O(log b)。原理是拆指数。假设 b 13二进制是 1101也就是 13 8 4 1。那么a^13 a^(841) a^8 × a^4 × a^1关键点在于a^8、a^4、a^1 这些项可以递推得到每做一次平方指数就翻倍。a^1 平方变成 a^2再平方变成 a^4再平方变成 a^8。所以我们只需要从低位到高位扫描 b 的二进制位把需要的那几项乘起来。标准代码typedef long long ll; ll qpow(ll a, ll b, ll mod) { ll res 1; a % mod; while (b 0) { if (b 1) res res * a % mod; a a * a % mod; b 1; } return res; }2.2 从整数到矩阵底数换掉骨架不变现在把目光从整数移到矩阵。整数的乘法换成矩阵乘法整数的单位元 1 换成单位矩阵 I循环逻辑完全一样。a^b 中 a 是整数矩阵快速幂中 base 是矩阵。为什么这个替换成立因为快速幂只依赖乘法的结合律。矩阵乘法满足结合律(A×B)×C A×(B×C)所以二进制的指数拆分在矩阵里同样成立。我在教学中喜欢用一个类比整数快速幂是用翻倍代替逐个乘矩阵快速幂就是用矩阵的翻倍代替逐次递推。底数从数字变成了矩阵但不断平方底数、按二进制位累乘结果的骨架没有任何变化。这样一来只要你会写整数快速幂矩阵快速幂的代码对你来说就只是把 res 从1换成I把 a 从整数换成矩阵把乘法函数替换成矩阵乘法。剩下的全是细节。3. 矩阵快速幂的标准模板结构体、乘法与幂运算的坑3.1 结构体与矩阵乘法的实现竞赛里写矩阵快速幂一般用结构体封装。这里给出一个我常用的、比较稳的模板Ctypedef long long ll; const ll MOD 1000000007LL; const int SZ 105; // 矩阵最大边长按题目调整 struct Mat { int n, m; // 行数, 列数 ll a[SZ][SZ]; Mat(int n 0, int m 0) : n(n), m(m) { memset(a, 0, sizeof(a)); } }; Mat mul(const Mat A, const Mat B) { Mat C(A.n, B.m); for (int i 0; i A.n; i) for (int k 0; k A.m; k) if (A.a[i][k]) for (int j 0; j B.m; j) C.a[i][j] (C.a[i][j] A.a[i][k] * B.a[k][j]) % MOD; return C; }这里有两个细节要解释。第一个是乘法循环的先后顺序。教科书里通常写 i, j, k 三层循环但实际跑起来 i, k, j 的顺序更快。原因是内存访问的局部性更好A 的第 i 行拿到后k 变化时 B 的访问是逐行扫的不会反复跳跃。虽然矩阵乘法复杂度还是 O(n^3)但常数能小不少。对于 100×100 的矩阵这个优化立竿见影。第二个是取模的位置。C.a[i][j] (C.a[i][j] A.a[i][k] * B.a[k][j]) % MOD;这行每次累加后立刻取模。两个小于 MOD 的数相乘乘积最大约 10^18还在 long long 范围内所以安全。如果累计好几项才取一次模long long 可能会爆这是后文要强调的重点。3.2 幂运算主函数单位矩阵是关键接下来是幂运算Mat qpow(Mat A, ll b) { Mat res(A.n, A.n); for (int i 0; i A.n; i) res.a[i][i] 1; // 单位矩阵 I while (b 0) { if (b 1) res mul(res, A); A mul(A, A); b 1; } return res; }初学者最容易在这里翻车res 初始化成零矩阵结果算出来全是 0。单位矩阵在矩阵乘法里的地位就是1在普通乘法里的地位任何矩阵乘单位矩阵结果还是它自己。所以 res 必须是对角线全 1 的 I。调用方式也很直白。还是斐波那契Mat M(2, 2); M.a[0][0] 1; M.a[0][1] 1; M.a[1][0] 1; M.a[1][1] 0; Mat res qpow(M, n - 1); // v(n) M^(n-1) * v(1) ll ans res.a[0][0]; // F(n) 结果矩阵第一行第一列前提是初始向量取 v(1) [F(1), F(0)]^T [1, 0]^T。这样 res 的第一行第一列就是 F(n)。如果初始向量定义成 v(0) [F(0), F(1)]^T幂次和对答案的取值都会变需要重新算。这种下标与幂次错位的问题第 6 节还会专门讲。4. 转移矩阵构造从递推式到矩阵的通用套路4.1 状态向量选哪几项一个经验法则很多选手能默写模板但一到给个递推式自己构造矩阵就懵了。构造矩阵的本质是回答两个问题状态向量里放哪些数转移矩阵每一行怎么写我的经验法则是状态向量里放下一轮要用的所有量。递推式右边出现什么左边下一轮还需要什么就统统放进去。以二阶线性递推为例F(n) a×F(n-1) b×F(n-2)为了从 v(n-1) 推出 v(n)v(n-1) 至少要包含 F(n-1) 和 F(n-2)。v(n) 应该包含 F(n) 和 F(n-1)。其中 F(n-1) 只是把旧向量里的 F(n-1) 平移过来靠转移矩阵第二行的 [1, 0] 实现F(n) 是两项的线性组合靠第一行的 [a, b] 实现。所以v(n) [F(n) ] M [ a b ] [F(n-1)] [ 1 0 ]每一行对应 v(n) 的一个分量每一列对应 v(n-1) 的一个分量。写矩阵时先确定状态向量有哪些分量然后一行一列地对照填系数。4.2 带常数项的递推式往向量里塞一个常数竞赛题不会只给这么简单的递推。常见的升级是带常数项比如F(n) a×F(n-1) b×F(n-2) c问题的关键常数项不是 F 的分量但矩阵乘法只能做线性组合没法直接加一个 c。解决办法是给状态向量增加一个恒等于 1 的维度。取 v(n-1) [F(n-1), F(n-2), 1]^T那么F(n) a*F(n-1) b*F(n-2) c*1 v(n) [F(n) ] M [ a b c ] [F(n-1)] [ 1 0 0 ] [1 ] [ 0 0 1 ]最后一行的 [0, 0, 1] 保证了常数维永远保持为 1。这样额外加一个常数的操作就被吸收进了矩阵乘法。这个技巧非常重要很多看起来无法矩阵化的递推其实只是缺一个常数为 1 的维度。我见过不少人卡在这里其实只要想明白乘法只能做带系数的求和不能凭空加数思路就通了。4.3 带多项式的递推式把 n 的多项式也拆成状态再难一点的题目递推式会出现 n 的多项式比如F(n) F(n-1) F(n-2) n^2这时候不能只放 F(n-1)、F(n-2)、1因为 n^2 会随 n 变化。但从 n 到 n1n^2、n、1 的变化是确定的(n1)^2 n^2 2n 1 (n1) n 1 1 1所以把 n^2、n、1 也放进状态向量。设v(n) [F(n), F(n-1), (n1)^2, n1, 1]^T注意我这里为了推导方便让 v(n) 里包含 (n1)^2这样 v(n-1) 里就包含 n^2。转移矩阵是F(n) F(n-1) F(n-2) n^2 - 系数 [1, 1, 1, 0, 0] F(n-1) F(n-1) - 系数 [1, 0, 0, 0, 0] (n1)^2 n^2 2n 1 - 系数 [0, 0, 1, 2, 1] n1 n 1 - 系数 [0, 0, 0, 1, 1] 1 1 - 系数 [0, 0, 0, 0, 1]写成矩阵就是M [ 1 1 1 0 0 ] [ 1 0 0 0 0 ] [ 0 0 1 2 1 ] [ 0 0 0 1 1 ] [ 0 0 0 0 1 ]到这里你应该能感觉到构造矩阵的核心不是背某个固定形式而是数清楚下一轮需要哪些量并把这些量的递推关系都写出来。只要能写成线性组合就能进矩阵。4.4 构造完怎么验证先手算再对拍构造完矩阵强烈建议先手算验证一次不要直接交题。方法很简单取 n 很小的值比如 n2把初始向量 v(1) 乘上 M看得到的 v(2) 是否符合递推式。例如斐波那契v(1) [F(1), F(0)]^T [1, 0]^T。M×v(1) [1, 1]^T对应 [F(2), F(1)]^TF(2)1正确。对于复杂递推手算一遍能发现八成错误。剩下的交给对拍——写一个暴力递推函数和矩阵快速幂结果对比跑几千组随机小数据。对拍的具体写法放在第 6 节。5. 邻接矩阵与路径计数矩阵快速幂的另一类应用5.1 恰好走 k 步的路径方案数矩阵快速幂不止用于递推数列图论里也有一个非常经典的应用路径计数。给定一个有向无权图用邻接矩阵 M 表示M[i][j] 1 表示存在 i 到 j 的边。那么 M^k[i][j] 恰好等于从 i 到 j 走恰好 k 步的路径条数。为什么以 k2 为例M^2[i][j] sum over t of M[i][t] × M[t][j]每一项 M[i][t]×M[t][j] 表示i 先到 tt 再到 j的一条路径t 就是中间节点。对所有中间节点求和就是 i 到 j 两步路径的总数。数学归纳法可以继续推到任意 k。这个结论非常漂亮因为它把枚举路径这种组合爆炸的问题变成了矩阵幂这种 O(log k) 的运算。唯一要注意的是路径数随 k 增长极快通常题目会要求在模意义下输出。5.2 扩展自环、无向图与不超过 k 步使用的时候有几个细节要小心。第一自环。如果图允许自环邻接矩阵 M[i][i] 要设为 1这样停在原地的路径也会被正确计算。第二无向图。无向图的邻接矩阵是对称的但路径计数不再区分方向。比如一个只有一条边的两点图从 1 到 2 走 1 步的方案数是 1走 2 步回到 1 的方案数是 11→2→1这些结果直接用矩阵算都对。但要注意无向图中来回走会导致路径数增长更快取模前数值会很大。第三不超过 k 步的问题。如果题目问的是从 i 到 j 走不超过 k 步的方案数直接求 M^1 M^2 ... M^k 是可行的但直接算需要 O(k log k) 次矩阵乘法k 一大又卡住。标准技巧是构建一个增广矩阵[ M I ] [ 0 I ]这个分块矩阵的 k 次幂的右上角恰好是 M^1 M^2 ... M^k。这个技巧在我带训练的时候考得不多但值得知道因为很多计数题本质是它的变体。5.3 大维度稀疏图的优化方向路径计数的矩阵维度是图的点数 n。如果 n 在 20 以下O(n^3 log k) 完全没问题n 到 100 还行n 到 1000 甚至更高普通矩阵乘法就扛不住了。这时候要利用稀疏性。图如果可以简化到每个点出边很少可以用遍历边集的方式代替矩阵乘法枚举中间节点 t对 t 的每条出边做累加。复杂度会降到 O(m log k)m 是边数。不过这种优化一般比赛里不常用我更愿意把它当作一种知道有这条路的知识储备。6. 专题训练中的高频错误与调试技巧6.1 单位矩阵初始化遗忘症前面提过res 要初始化为单位矩阵 I不是零矩阵。这是新人在这个专题里出现频率最高的错误没有之一。为什么不对因为零矩阵在乘法里扮演的是0的角色任何矩阵乘零矩阵结果都是零矩阵。如果你把 res 写成全 0那么无论 b 多大结果始终是 0然后你的答案就永远输出 0。这个 bug 的关键特征是跑小样例有时恰好没错一放大 n 就全错。我建议每个选手在模板里把单位矩阵初始化的那几行做成肌肉记忆。for (int i 0; i A.n; i) res.a[i][i] 1;6.2 取模与溢出long long 能装下的极限第二个高频问题是溢出。很多人写着写着就忘了取模或者取模太晚。long long 的最大值大约是 9.22×10^18。MOD 取 10^97 时两个小于 MOD 的数相乘乘积约 10^18刚好在 long long 范围内不会爆。但如果代码写成C.a[i][j] A.a[i][k] * B.a[k][j]; // 累加多次后才取模这里 C.a[i][j] 会连续累加多个乘积比如矩阵维度是 100一层循环下来就加了 100 个约 10^18 的数必然溢出。正确写法是每加一项就取模一次或者用 unsigned long long 配合周期取模。记住矩阵乘法里没有攒着再取模这种操作除非你确定维数很小且每一项都远小于极限。6.3 幂次与下标的错位第三个问题是幂次写错。求斐波那契第 n 项时v(n) M^(n-1) × v(1)不是 M^n × v(1)。很多人把 n 直接丢进快速幂导致结果从 F(n1) 开始或者干脆算错。我推荐一个自查方法先规定清楚初始向量对应哪一代。如果你的状态向量是 v(0) [F(0), F(1)]^T那 F(n) (M^n × v(0))[0]如果你的状态向量是 v(1) [F(1), F(0)]^T那 F(n) (M^(n-1) × v(1))[0]。写代码前先把式子写在注释里再动手写能避免很多无意义的 WA。还有 n 0 的边界情况。有的题目 n 从 0 开始如果直接做 n-1 次幂就会得到负数指数。这时要么特判要么把初始向量重新定义成 v(0) 形式。边界测试一定要覆盖 n0, n1, n2 三档。6.4 对拍最笨也最有效的调试法我在专题训练里反复强调对拍。矩阵快速幂的题目一旦 WA肉眼很难看出是构造矩阵错还是模板错最好的办法就是写一个暴力递推函数然后随机生成小 n 去对比。暴力函数可以是这样的ll fib_bruteforce(int n) { ll f[100] {0}; f[1] 1; for (int i 2; i n; i) f[i] f[i-1] f[i-2]; return f[n]; }主函数里随机测 n 0..100比较矩阵快速幂结果和暴力结果。如果暴力都对、矩阵结果错说明矩阵构造或调用逻辑有问题如果连暴力都不对说明你递推式本身就写错了。对拍脚本不需要多复杂一个随机数生成器 两个函数 一个循环就够。对小样例多测几组之后再提交到题库里。这一步能帮你省下大量反复提交试错的时间非常值得养成习惯。7. 带完这期专题我再啰嗦三件事第一模板一定要自己默写过一遍不要直接复制。默写过程中你会把单位矩阵初始化先取模再相加行与列的对应关系这些细节刻进脑子。等到比赛的时候这些细节不会背叛你。第二构造转移矩阵的能力靠手推。我建议每个专题至少找五道不同风格的题一阶含常数、二阶含多项式、带前缀和的递推、邻接矩阵路径计数、带分块矩阵的求和。每道题都在纸上把状态向量写出来把矩阵填出来再用小数据验证。练完这五道你对线性变换的理解会上升一个台阶。第三一个小技巧写矩阵前先写目标公式。比如你想让 F(n) 某行 × 某向量先把这一行该填的系数竖着列出来再去填代码。这样行和列不容易弄反。我在训练时见过太多人把转移矩阵写成转置如果先在纸上写清楚这个问题能减少一半以上。矩阵快速幂本质上不是玄学它只是把重复的线性变换压缩成了快速幂。理解到这一层你遇到递推式时就不会再害怕反而会想这个递推能不能也写成一个矩阵能想到这里专题训练的目的就达到了。