OI算法——矩阵加速递推

在开始之前，确保你至少已经学会了矩阵乘法的计算方法。

矩阵加速递推（后边简称矩阵递推）充分利用了初等矩阵的性质，可以将原本耗时间费空间的函数递归、优化但复杂度较高的记忆化搜索进一步加速为复杂度仅 $\mathcal O(m^3\log n)$ 的对数级别算法（ $m$ 一般为 ${2}$ 或 ${3}$ ）。在数据范围很大时可以考虑使用。

矩阵乘法是如何起作用的#

初等矩阵——概念及用法，花一点篇幅来复习一下：

初等行变换：

交换矩阵某两行
将某一行的元素全部乘以一个非零数
将某一行的非零倍加到另一行上

非常简单，甚至我们化简多元方程式都会用到上边的变换。事实上，如果把一次方程组的系数和常数项按一定顺序排列起来，将会得到一个系数矩阵，系数矩阵经过一系列化简和反代也可以解出原方程组的解。

考虑一个单位矩阵 $I=\begin{bmatrix}1&0&0\\0&1&0\\0&0&1\end{bmatrix}$ ，它很特殊，和实数运算的 $1$ 类似：任何矩阵左乘和右乘单位矩阵所得的乘积矩阵都是这个矩阵本身。顺便说一句，矩阵乘法一般不符合乘法交换律，但单位矩阵乘法除外，也就是说当矩阵长和宽符合要求时，就有 $AI=IA=A$ 成立。

那么一个初等矩阵就是单位矩阵只进行一次初等行变换得到的矩阵。进行变换法则的第几条就是第几类初等矩阵。一般来说，矩阵递推里的转移矩阵不属于初等矩阵，因为它通常会经过不止一次的初等行变换。递推时，如果初始矩阵是一列数，那么一般选择左乘转移矩阵；如果是一行数字，就选择右乘。（具体情况具体分析）

接下来进入正题：

假如有下边这个递推关系

f(x)=\begin{cases} 1&x=1 \\2&x=2 \\f(x-1)+f(x-2)&x>2 \end{cases}

初始情况下，我们向矩阵中放入两个元素： $A=\begin{bmatrix}1&2\end{bmatrix}$ 。我们希望构造出一个转移矩阵使得转移后的结果就是前两个元素相加的和（根据递推方程），做一次乘法，矩阵变为 $\begin{bmatrix}2&3\end{bmatrix}$ ；再来一次，变成 $\begin{bmatrix}3&5\end{bmatrix}$ 。不难发现，当做了 $n$ 次乘法后，矩阵变为 $\begin{bmatrix}f(n+1)&f(n+2)\end{bmatrix}$ 。所以假设我们的转移矩阵是 $M$ ， $f(n)=s_{11}$ ，其中矩阵 $S=AM^{n-1}$ 。就基本搞清楚了矩阵递推的原理了。

矩阵结构体#

函数I()的功能是构造单位矩阵，后面会涉及到定义它的原因。整体思路就是定义二维数组存放矩阵元素，声明时对内部元素自动置零，以及构造单位矩阵

1
struct Matrix {
2
  ll mat[N + 1][N + 1];
3

4
  Matrix() {
5
      memset(mat, 0, sizeof mat);
6
  }
7

8
  void I() {
9
            memset(mat, 0, sizeof mat);
10
      for (int i = 1; i <= N; i++) mat[i][i] = 1;
11
  }
12
};

矩阵快速幂#

和实数运算一样，矩阵自乘也可以用二进制快速幂的方式快速求解，复杂度是 $\mathcal O(\log n)$ 的。

下边是实数快速幂的代码（带取模）：

1
int qpow(int a, int b) {
2
  int res = 1;
3
  while (b) {
4
    if (b & 1) res = (ll) res * a % MOD;
5
    a = a * a % MOD;
6
    b >>= 1;
7
  }
8
  return res;
9
}

原理不再赘述。如果想要把它改造成适用于矩阵的快速幂算法，我们就需要实现这两个运算：

矩阵乘法
矩阵置一

对于第一点，我们根据矩阵乘法的定义可以很轻松写出代码。一般来说，重载运算符是一个很方便的办法（注意如果重载在结构体内就需要在函数定义时加上friend友元访问权）：

1
Matrix operator *(const Matrix &l, const Matrix &r) {
2
  Matrix res;
3
  for (int i = 1; i <= N; i++) {
4
    for (int j = 1; j <= N; j++) {
5
      for (int k = 1; k <= N; k++) {
6
        res.mat[i][j] += (l.mat[i][k] * r.mat[k][j]);
7
        res.mat[i][j] %= MOD;
8
      }
9
    }
10
  }
11
}

那么“置一”是什么呢？

相当于实数快速幂里的int res = 1;，试想一个全新定义的零矩阵（元素全为零）拿去计算乘法，最终的结果总会是零。因此我们就需要找到和 $1$ 作用相同的矩阵，赋上初值，计算出来的结果才是正确的。于是很自然想到了单位矩阵，结构体中的I()函数用于将矩阵变为一个主对角线元素全为 $1$ 、其他元素均为 $0$ 的矩阵。

类比实数快速幂，矩阵快速幂是这样的：

1
Matrix qpow(Matrix a, int b) {
2
  Matrix res;
3
  res.I();
4
  while (b) {
5
      if (b & 1) res = res * a;
6
      a = a * a;
7
      b >>= 1;
8
  }
9
  return res;
10
}

此时就不能图省事用*=运算符了（除非你另外重载，但这样会更麻烦）

转移矩阵的构造#

转移矩阵可谓是矩阵递推题目的灵魂所在，合理地构造转移矩阵可以达到事半功倍的效果。接下来通过几个例子来深入探究转移矩阵的构造方法：

洛谷 P1962 斐波那契数列#

题目传送门：这里

题目难度：普及+/提高

大家都知道，斐波那契数列是满足如下性质的一个数列：
$F_n = \left\{\begin{aligned} 1 \space (n \le 2) \\ F_{n-1}+F_{n-2} \space (n\ge 3) \end{aligned}\right.$
请你求出 $F_n \bmod 10^9 + 7$ 的值。
输入一行一个正整数 $n$
输出一行一个整数表示答案。
数据范围：
对于 $60\%$ 的数据， $1\le n \le 92$ ；
对于 $100\%$ 的数据， $1\le n < 2^{63}$ 。

正如这头图所示，我们的初始矩阵是 $A=\begin{bmatrix}1&1\end{bmatrix}$ ，转移矩阵是 $M=\begin{bmatrix}0&1\\1&1\end{bmatrix}$ 。 $\operatorname{fib}_n$ 就是 $AM^{n-1}$ 的第一个元素。~~嘎嘎好用是不是~~

首先，要求出某一项，就必须明确它的前两项。因此我们让初始矩阵填上 $\operatorname{fib}_1=1$ 和 $\operatorname{fib}_2=1$ 。对于转移，我们有一个非常好的小技巧：整体左移——这就好比一个滑动窗口（但不是单调队列那个），斐波那契数列按顺序排列在一起： $\operatorname{fib}_1,\operatorname{fib}_2,\operatorname{fib}_3,\operatorname{fib}_4,\dots\operatorname{fib}_{n-1},\operatorname{fib}_n$ 。最开始这个矩阵框住了 $\operatorname{fib}_1$ 和 $\operatorname{fib}_2$ ，操作一次，它框住 $\operatorname{fib}_2$ 和 $\operatorname{fib}_3$ ，以此类推……每次挪一下，因此做了 $n-1$ 次乘法后就挪到了 $\operatorname{fib}_n$ 的位置。那我们怎么构造这种转移矩阵呢？

所谓左移，就是让 $a_{12}$ 换到 $a_{11}$ 的位置来，根据矩阵乘法系数配对的原理。对于转移矩阵第一列，就是下面这样：

M=\begin{bmatrix} 0& \\1& \end{bmatrix}

第二列也很简单，对原先的转移矩阵的两个元素都配上系数 $1$ 即可：

M=\begin{bmatrix}0&1\\1&1\end{bmatrix}

于是代码就出来了（~~不开long long见祖宗~~）：

1
#include <bits/stdc++.h>
2

3
#define N 15
4
using namespace std;
5

6
typedef long long ll;
7

8
const int MOD = 1e9 + 7;
9

10
struct Matrix {
11
    ll mat[N][N];
12

13
    Matrix() {
14
        memset(mat, 0, sizeof mat);
15
    }
16

17
    void I() {
18
        memset(mat, 0, sizeof mat);
19
        for (int i = 1; i <= 2; i++) mat[i][i] = 1;
20
    }
21
};
22

23
Matrix operator *(const Matrix &l, const Matrix &r) {
24
    Matrix res;
25
    for (int i = 1; i <= 2; i++) {
26
        for (int j = 1; j <= 2; j++) {
27
            for (int k = 1; k <= 2; k++) {
28
                res.mat[i][j] += (l.mat[i][k] * r.mat[k][j]);
29
                res.mat[i][j] %= MOD;
30
            }
31
        }
32
    }
33
    return res;
34
}
35

36
Matrix qpow(Matrix a, ll b) {
37
    Matrix res;
38
    res.I();
39
    while (b) {
40
        if (b & 1) res = res * a;
41
        a = a * a;
42
        b >>= 1;
43
    }
44
    return res;
45
}
46

47
int main() {
48
    ios::sync_with_stdio(false);
49
    cin.tie(nullptr);
50
    cout.tie(nullptr);
51

52
    ll n;
53
    cin >> n;
54

55
    Matrix A, M;
56
    A.mat[1][1] = A.mat[1][2] = 1;
57
    M.mat[1][2] = M.mat[2][1] = M.mat[2][2] = 1;
58

59
    A = A * qpow(M, n - 1);
60

61
    cout << A.mat[1][1] % MOD << endl;
62

63
    return 0;
64
}

总用时： $42ms$ 记录

这种类型的题还有洛谷 P1349 广义斐波那契数列，同样是左移技巧，只不过转移矩阵要稍作变动。

洛谷 P1397 矩阵游戏#

题目传送门：这里

题目难度：提高+/省选

题目来源：NOI 2013

~~NOI 2012 和 2013连着两年都考了矩阵递推，真的强！~~

婷婷是个喜欢矩阵的小朋友，有一天她想用电脑生成一个巨大的 $n$ 行 $m$ 列的矩阵（你不用担心她如何存储）。她生成的这个矩阵满足一个神奇的性质：若用 $F[i,j]$ 来表示矩阵中第 $i$ 行第 $j$ 列的元素，则 $F[i,j]$ 满足下面的递推式:
$\begin{aligned} F[1, 1] &= 1 \\ F[i, j] &=a\times F[i, j-1]+b, &j\neq 1 \\ F[i, 1] &=c\times F[i-1, m]+d, &i\neq 1 \\ \end{aligned}$
递推式中 $a,b,c,d$ 都是给定的常数。
现在婷婷想知道 $F[n,m]$ 的值是多少，请你帮助她。由于最终结果可能很大，你只需要输出 $F[n,m]$ 除以 ${10^9+7}$ 的余数。
输入包含一行，有六个整数 $n,m,a,b,c,d$ 。意义如题所述。
输出包含一个整数，表示 $F[n,m]$ 除以 $10^9+7$ 的余数。
数据范围：
${1\leq n,m\leq10^{1~000~000};~1\leq a,b,c,d\leq10^9}$

这道题需要我们推一个式子，因为递推公式出现了两种情况，我们就需要两个不同的转移矩阵。假设一个为 $M_{ab}$ ，一个为 $M_{cd}$ ，分别对应递推式里的系数，以及一个初始矩阵 $A$ 。

这里出现了常数项，通常选择在初始矩阵中放入一个常量 $1$ ，每次递推将它乘以这个常数、并且需要保证它不被转移矩阵改变，这样一来才能保证递推稳定运行。

那么在初始矩阵的第二列放上常量 $1$ ，以 $M_{ab}$ 为例，要满足递推关系，未知数和 $1$ 的系数分别是 $a$ 和 $b$ ，于是转移矩阵第一列就是 $\begin{bmatrix}a\\b\end{bmatrix}$ ，第二列就是 $\begin{bmatrix}0\\1\end{bmatrix}$ 。

同理有： $A=\begin{bmatrix}1&1\end{bmatrix}$ 、 $M_{ab}=\begin{bmatrix}a&0\\b&1\end{bmatrix}$ 、 $M_{cd}=\begin{bmatrix}c&0\\d&1\end{bmatrix}$ 。

接下来根据题目描述，要想一路推到右下角的 $(n,m)$ ，首先就得把 $(n-1,m)$ 弄出来，而这一行又从 $(n-1,1)$ 递推得来，以此类推……那么每一行最右边的元素和该行第一个元素的关系就是 $AM_{ab}^{m-1}$ ，又因为矩阵共 $n$ 行，每一行的开头还得乘上一个 $M_{cd}$ ，因此公式就是： $A(M_{ab}^{m-1}M_{cd})^{n-1}M_{ab}^{m-1}$ 。再看看数据范围，~~十的一百万次方？？？太抽象了~~，对于这么大的幂，普通的位运算快速幂已经满足不了时限了，于是我们引入一种高级方法——十进制快速幂：

快速幂基于数字的拆位，所以我们可以选择在十进制表示下拆位运算。因此就算是十的一百万次方，应用十进制快速幂就会让复杂度降落不少，因此我们试验这个方法：

十进制矩阵快速幂：

1
Matrix dec_qpow(Matrix a, string b) {
2
  Matrix res;
3
  res.I();
4
  int len = b.length();
5
  while (len) {
6
    int p = b[len - 1] - '0';
7
    if (p) {
8
      for (int i = 1; i <= p; i++) {
9
        res = res * a;
10
      }
11
    }
12
    for (int i = 1; i <= 10; i++) a = a * a;
13
    len--;
14
  }
15
  return res;
16
}

当然可以用二进制快速幂取代中间的循环乘幂，代码会简洁一些：

1
Matrix dec_qpow(Matrix a, string b) {
2
  Matrix res;
3
  res.I();
4
  int len = b.length();
5
  while (len) {
6
    int p = b[len - 1] - '0';
7
    res = res * bin_qpow(a, p);
8
    a = bin_qpow(a, 10);
9
    len--;
10
  }
11
  return res;
12
}

然后我们再来算上面推出的式子，考虑到变量 $n$ 和 $m$ 必须要用字符串读入，可以用高精度的思想对二者进行预处理，将它们的值分别减少 ${1}$ ，然后将处理后的字符串作为参数传入快速幂函数中计算即可。下边是代码：

1
#include <bits/stdc++.h>
2
#define N 5
3
using namespace std;
4

5
const int MOD = 1e9 + 7;
6

7
typedef long long ll;
8

9
struct Matrix {
10
  ll a[N][N];
11

12
  Matrix() {
13
    memset(a, 0, sizeof a);
14
  }
15

16
  void I() {
17
    a[1][1] = a[2][2] = 1;
18
  }
19
} A, M, S;
20

21
Matrix operator *(const Matrix &l, const Matrix &r) {
22
  Matrix res;
23
  for (int i = 1; i <= 2; i++) {
24
    for (int j = 1; j <= 2; j++) {
25
      for (int k = 1; k <= 2; k++) {
26
        res.a[i][j] = (res.a[i][j] + l.a[i][k] % MOD * r.a[k][j] % MOD) % MOD;
27
      }
28
    }
29
  }
30
  return res;
31
}
32

33
Matrix bin_qpow(Matrix a, ll b) {
34
  Matrix res;
35
  res.I();
36
  while (b) {
37
    if (b & 1) res = res * a;
38
    a = a * a;
39
    b >>= 1;
40
  }
41
  return res;
42
}
43

44
Matrix dec_qpow(Matrix a, string b) {
45
  Matrix res;
46
  res.I();
47
  int len = b.length();
48
  while (len) {
49
    int p = b[len - 1] - '0';
50
    for (int i = 1; i <= p; i++) {
51
      res = res * a;
52
    }
53
    a = bin_qpow(a, 10);
54
    len--;
55
  }
56
  return res;
57
}
58

59
string init(string s) {
60
  for (int i = s.length() - 1; i >= 0; i--) {
61
    if (s[i] == '0') s[i] = '9';
62
    else {
63
      s[i]--;
64
      break;
65
    }
66
  }
67
  return s;
68
}
69

70
int main() {
71
  ios::sync_with_stdio(false);
72
  cin.tie(nullptr);
73
  cout.tie(nullptr);
74

75
  string n, m;
76
  ll a, b, c, d;
77

78
  cin >> n >> m >> a >> b >> c >> d;
79

80
  A.a[1][1] = 1, A.a[1][2] = 1;
81
  M.a[2][2] = 1, M.a[1][1] = a, M.a[2][1] = b;
82
  S.a[2][2] = 1, S.a[1][1] = c, S.a[2][1] = d;
83

84
  n = init(n), m = init(m);
85

86
  Matrix T = A * dec_qpow((dec_qpow(M, m) * S), n) * dec_qpow(M, m);
87

88
  cout << T.a[1][1] % MOD << endl;
89
  return 0;
90
}

这段代码只能得80分，最后四个点TLE了。此时有两种解决办法，第一是卡常——将矩阵乘法的三层循环完全展开、内联函数、加法取模，可以将代码运行时间压到 ${1.37}s$ （最慢的点 ${341}ms$ ）记录。

矩阵递推的变式#

广义矩阵乘法#

相信你已经背下了矩阵乘法的模板了，为了避免遗忘，再在这里给出矩阵乘法的一般定义：假设一个 $n\times p$ 的矩阵 $A$ 和一个 $p\times m$ 矩阵 $B$ 相乘，就有：

(AB)_{ij}=\sum\limits_{k=1}^{p}a_{ik}b_{kj}

最终结果 $C=AB$ ，满足 $C$ 为一个 $n\times m$ 矩阵。

那么是不是非得要求和和相加呢？当然不是，对于一般的矩阵乘法，才会像上边一样对乘积求和。广义矩阵乘法不仅限于上边的“对积求和”的规则，它还可以做到“对位与求位或和”。接下来探讨广义矩阵乘法所要满足的条件：

定义运算符 $\bigoplus$ 为“异加”（变种加法）、 $\bigotimes$ 为“异乘”（变种乘法）。假如还是上面的 $A,B$ 矩阵，那么广义矩阵乘法将是这种形式：

(AB)_{ij}=\bigoplus\limits_{k=1}^{p}a_{ik}\otimes b_{kj}

要想让这种运算关系能够支持递推算法，那就必须使算式 $A\times M\times M\times\dots\times M=AM^n$ 成立才行。而这个转化很明显是将连乘的 $M$ 结合为矩阵的幂，因此只要上边的新定义算式满足乘法结合律，就可以被归入广义矩阵乘法的范畴内。假设此时有 $A,B,C$ 三个能够相乘的矩阵，那么：

\begin{aligned} [(AB)C]_{ij}&=\bigoplus\limits_{k_1=1}^{n}(AB)_{ik_1}\otimes C_{k_1j} \\&=\bigoplus\limits_{k_1=1}^{n}\left(\bigoplus\limits_{k_2=1}^{n}a_{ik_2}\otimes b_{k_2k_1}\right)\otimes c_{k_1j} \\ [A(BC)]_{ij}&=\bigoplus_{k_1=1}^{n}A_{ik_1}\otimes(BC)_{k_1j} \\&=\bigoplus\limits_{k_1=1}^{n}a_{ik_1}\otimes\left(\bigoplus\limits_{k_2=1}^{n}b_{k_1k_2}\otimes c_{k_2j}\right) \end{aligned}

换句话说：

\bigoplus\limits_{k_1=1}^{n}\left(\bigoplus\limits_{k_2=1}^{n}a_{ik_2}\otimes b_{k_2k_1}\right)\otimes c_{k_1j}=\bigoplus\limits_{k_1=1}^{n}a_{ik_1}\otimes\left(\bigoplus\limits_{k_2=1}^{n}b_{k_1k_2}\otimes c_{k_2j}\right)

类比一般加法和乘法。当 $\bigotimes$ 运算符合交换律、结合律；并且 $\left(\bigoplus a\right)\otimes b=\bigoplus\left(a\otimes b\right)$ 成立，则原式化简为：

\left[A(BC)\right]_{ij}=\left[(AB)C\right]_{ij}=\bigoplus\limits_{k_1=1}^{n}\bigoplus\limits_{k_2=1}^{n}\left(a_{ik_2}\otimes b_{k_2k_1}\otimes c_{k_1j}\right)

音乐

JustPureH₂O

穷方圆平直之情，尽规矩准绳之用

音乐