学习笔记：矩阵快速幂

· · 算法·理论

矩阵快速幂

引入

斐波那契数列想必大家都知道，它的递推公式为：

f_n=f_{n-1}+f_{n-2}

要求 f_n 的话，首先想到的是线性递推，这种方法的时间复杂度是 O(n)，在大部分的情况下都是够用的。然而对于某些简单（毒瘤）题，它们的数据规模较大，采用线性递推的方法显然还是不够，这个时候就有必要学习一下矩阵快速幂了。

前置知识

快速幂

定义

快速幂，二进制取幂（Binary Exponentiation，也称平方法），是一个在 O(\log n) 的时间内计算 a^n 的小技巧，而暴力的计算需要 O(n) 的时间。

这个技巧也常常用在非计算的场景，因为它可以应用在任何具有结合律的运算中。其中显然的是它可以应用于模意义下取幂、矩阵幂等运算，我们接下来会讨论。

解释

计算 a 的 n 次方表示将 n 个 a 乘在一起：a^{n} = \underbrace{a \times a \cdots \times a}_{n\text{ 个 a}}。然而当 a,n 太大的时侯，这种方法就不太适用了。不过我们知道：a^{b+c} = a^b \cdot a^c,\,\,a^{2b} = a^b \cdot a^b = (a^b)^2。二进制取幂的想法是，我们将取幂的任务按照指数的 二进制表示 来分割成更小的任务。

过程

首先我们将 n 表示为 2 进制，举一个例子：

3^{13} = 3^{(1101)_2} = 3^8 \cdot 3^4 \cdot 3^1

因为 n 有 \lfloor \log_2 n \rfloor + 1 个二进制位，因此当我们知道了 a^1, a^2, a^4, a^8, \dots, a^{2^{\lfloor \log_2 n \rfloor}} 后，我们只用计算 \Theta(\log n) 次乘法就可以计算出 a^n。

于是我们只需要知道一个快速的方法来计算上述 3 的 2^k 次幂的序列。这个问题很简单，因为序列中（除第一个）任意一个元素就是其前一个元素的平方。

因此为了计算 3^{13}，我们只需要将对应二进制位为 1 的整系数幂乘起来就行了：

3^{13} = 6561 \cdot 81 \cdot 3 = 1594323

将上述过程说得形式化一些，如果把 n 写作二进制为 (n_tn_{t-1}\cdots n_1n_0)_2，那么有：

n = n_t2^t + n_{t-1}2^{t-1} + n_{t-2}2^{t-2} + \cdots + n_12^1 + n_02^0

其中 n_i\in\{0,1\}。那么就有

\begin{aligned} a^n & = (a^{n_t 2^t + \cdots + n_0 2^0})\\\\ & = a^{n_0 2^0} \times a^{n_1 2^1}\times \cdots \times a^{n_t2^t} \end{aligned}

根据上式我们发现，原问题被我们转化成了形式相同的子问题的乘积，并且我们可以在常数时间内从 2^i 项推出 2^{i+1} 项。

这个算法的复杂度是 O(\log n) 的，我们计算了 O(\log n) 个 2^k 次幂的数，然后花费 O(\log n) 的时间选择二进制为 1 对应的幂来相乘。

矩阵乘法

矩阵的乘法是向量内积的推广，矩阵相乘只有在第一个矩阵的列数和第二个矩阵的行数相同时才有意义，设 A 为 P \times M 的矩阵，B 为 M \times Q 的矩阵，设矩阵 C 为矩阵 A 与 B 的乘积，其中矩阵 C 中的第 i 行第 j 列元素可以表示为：

C_{i,j} = \sum_{k=1}^MA_{i,k}B_{k,j}

在矩阵乘法中，结果 C 矩阵的第 i 行第 j 列的数，就是由矩阵 A 第 i 行 M 个数与矩阵 B 第 j 列 M 个数分别 相乘再相加 得到的。这里的 相乘再相加，就是向量的内积。乘积矩阵中第 i 行第 j 列的数恰好是乘数矩阵 A 第 i 个行向量与乘数矩阵 B 第 j 个列向量的内积，口诀为 左行右列。

线性代数研究的向量多为列向量，根据这样的对矩阵乘法的定义方法，经常研究对列向量左乘一个矩阵的左乘运算，同时也可以在这里看出「打包处理」的思想，同时处理很多个向量内积。矩阵乘法满足结合律，不满足一般的交换律。利用结合律，矩阵乘法可以利用快速幂的思想来优化。

具体做法

首先来看一道题。

RABBIT1 - Counting Rabbits

题目大意

$$F_n = \left\{\begin{aligned} 1 \space (n \le 1) \\ F_{n-1}+F_{n-2} \space (n\ge 2) \end{aligned}\right.$$ $1\le T\le 100$，$1\le n\le 2147483647$，$1\le m\le 20$。 #### 思路 显然是斐波那契数列，在斐波那契数列（Fibonacci Sequence）当中，$F_1 = F_2 = 1$，$F_i = F_{i - 1} + F_{i - 2}(i \geq 3)$。 注意到 $1\le n\le 2147483647$，如果用线性递推的话，时间复杂度 $O(n)$ 稳稳 TLE，可以考虑矩阵快速幂优化递推。 #### 推导过程 设 $F(n)$ 表示一个 $1 \times 2$ 的矩阵 $\left[ \begin{array}{ccc}F_n & F_{n-1} \end{array}\right]$。我们希望根据 $F(n-1)=\left[ \begin{array}{ccc}F_{n-1} & F_{n-2} \end{array}\right]$ 推出 $F(n)$。 试推导一个转移矩阵 $\text{base}$，使 $F(n-1) \times \text{base} = F(n)$，即 $\left[\begin{array}{ccc}F_{n-1} & F_{n-2}\end{array}\right] \times \text{base} = \left[ \begin{array}{ccc}F_n & F_{n-1} \end{array}\right]$。 怎么推导转移矩阵呢？因为 $F_n=F_{n-1}+F_{n-2}$，所以 $\text{base}$ 矩阵第一列应该是 $\left[\begin{array}{ccc} 1 \\ 1 \end{array}\right]$，这样在进行矩阵乘法运算的时候才能令 $F_{n-1}$ 与 $F_{n-2}$ 相加，从而得出 $F_n$。同理，为了得出 $F_{n-1}$，矩阵 $\text{base}$ 的第二列应该为 $\left[\begin{array}{ccc} 1 \\ 0 \end{array}\right]$。 综上所述：$\text{base} = \left[\begin{array}{ccc} 1 & 1 \\ 1 & 0 \end{array}\right]$ 原式化为 $\left[\begin{array}{ccc}F_{n-1} & F_{n-2}\end{array}\right] \times \left[\begin{array}{ccc} 1 & 1 \\ 1 & 0 \end{array}\right] = \left[ \begin{array}{ccc}F_n & F_{n-1} \end{array}\right]$。 转化为代码，应该怎么求呢？ 定义初始矩阵 $\text{ans} = \left[\begin{array}{ccc}F_2 & F_1\end{array}\right] = \left[\begin{array}{ccc}1 & 1\end{array}\right], \text{base} = \left[\begin{array}{ccc} 1 & 1 \\ 1 & 0 \end{array}\right]$。那么，$F_n$ 就等于 $\text{ans} \times \text{base}^{n-1}$ 这个矩阵的第一行第一列元素，也就是 $\left[\begin{array}{ccc}1 & 1\end{array}\right] \times \left[\begin{array}{ccc} 1 & 1 \\ 1 & 0 \end{array}\right]^{n-1}$ 的第一行第一列元素。 注意到矩阵乘法不满足交换律，所以一定不能写成 $\left[\begin{array}{ccc} 1 & 1 \\ 1 & 0 \end{array}\right]^{n-1} \times \left[\begin{array}{ccc}1 & 1\end{array}\right]$ 的第一行第一列元素。另外，对于 $n \leq 1$ 的情况，直接输出 $1$ 即可，不需要执行矩阵快速幂。 为什么要乘上 $\text{base}$ 矩阵的 $n-1$ 次方而不是 $n$ 次方呢？因为 $F_1$ 是不需要进行矩阵乘法就能求的。也就是说，如果只进行一次乘法，就已经求出 $F_2$ 了。如果还不是很理解为什么幂是 $n-1$，建议手算一下。 利用矩阵快速幂，我们可以将时间复杂度由 $O(tn)$ 优化到 $O(t\log n)$，稳过。 **注意每一组数据处理完之后都要初始化矩阵。** ## 代码 ```c++ #include <iostream> #include <cstring> #define int long long // 不开 long long 见祖宗 using namespace std; int T, n, m; struct Matrix{ int a[3][3]; Matrix(){memset(a, 0, sizeof(a));} // 初始化矩阵 Matrix operator*(const Matrix &b)const{ // 利用重载运算符定义矩阵乘法 Matrix res; for(int i = 1 ; i <= 2 ; i ++) for(int j = 1 ; j <= 2 ; j ++) for(int k = 1 ; k <= 2 ; k ++) res.a[i][j] = (res.a[i][j] + a[i][k] * b.a[k][j]) % m; return res; } }ans, base; // ans 为答案矩阵，base 为转移矩阵 signed main(){ ios::sync_with_stdio(false);cin >> T; while(T--){ // 循环处理每一组数据 cin >> n >> m;m = 1 << m;n -= 1; if(n <= 1){puts("1");continue;} // 对于 n <= 1 的情况可以直接输出，无需调用矩阵快速幂 base.a[1][1] = 1;base.a[1][2] = 1; base.a[2][1] = 1;base.a[2][2] = 0; ans.a[1][1] = 1;ans.a[1][2] = 1; // 初始化答案矩阵以及转移矩阵 while(n > 0){ // 快速幂优化矩阵乘法 if(n & 1)ans = ans * base; base = base * base;n >>= 1; } cout << ans.a[1][1] << endl; // 输出答案 } return 0; } // 完结撒花～～ ``` ## 参考题目 这里给出一些相关题目。 [P1226 【模板】快速幂 | 取余运算](https://www.luogu.com.cn/problem/P1226) [P3390 【模板】矩阵快速幂](https://www.luogu.com.cn/problem/P3390) [P1939 【模板】矩阵加速（数列）](https://www.luogu.com.cn/problem/P1939) [P1349 广义斐波那契数列](https://www.luogu.com.cn/problem/P1349) [P1306 斐波那契公约数](https://www.luogu.com.cn/problem/P1306) [P1397 [NOI2013] 矩阵游戏](https://www.luogu.com.cn/problem/P1397)