本篇将介绍矩阵快速幂

矩阵介绍

在线性代数中有一个很重要的内容——矩阵（Matrix），对于矩阵 AA ，有 nn 行 mm 列，记作 [abcd]\begin{bmatrix} a&b \\ c&d \end{bmatrix}

一般使用 II 表示单位矩阵，乘以任何矩阵答案都为对方，其主对角线上为 11 其余为 00，主对角线是指 Ai,iA_{i,i} 的元素，例如 [1001]\begin{bmatrix} 1&0 \\ 0&1 \end{bmatrix}

矩阵乘法

两个矩阵相乘时，只有第一个矩阵的列数和第二个矩阵的行数相同时才有意义。
所以设矩阵 AA 有 PP 行 MM 列，矩阵 BB 有 MM 行 QQ 列，答案矩阵为 CC

矩阵的乘法是向量内积的推广，向量内积为将向量的每一位相乘后相加
同样，矩阵 Ci,jC_{i,j} 为 AA 的第 ii 行的每个数 与 BB 的第 jj 列的每个数相乘后相加

即为：Ci,j=∑k=1MAi,kBk,jC_{i,j} = \sum_{k=1}^{M}A_{i,k}B{k,j}

注意！乘法具有结合律，但不具有交换律，AA 和 BB 交换后结果不一致

利用结合律，矩阵乘法可以使用 快速幂 的思想来优化，这里我们直接使用例子来说明

矩阵乘法加速递推

线性递推式可以通过倒推矩阵 MM，计算并表示成矩阵乘法的形式
从复杂的递推式变为连乘之后可以使用矩阵快速幂的方式，来加速求解线性递推数列的某一项

我们以 斐波那契数列 题目为例：

斐波那契数列

斐波那契数列是满足如下性质的一个数列：

Fn={1 (n≤2)Fn−1+Fn−2 (n≥3)F_n = \left\{\begin{aligned} 1 \space (n \le 2) \\ F_{n-1}+F_{n-2} \space (n\ge 3) \end{aligned}\right.

请你求出 Fn mod 109+7F_n \bmod 10^9 + 7 的值。
【数据范围】1≤n<2631\le n < 2^{63}。

来源：洛谷P1962

显然，如果使用正常的递推/递归求法，时间条件无法满足，这个时候就要使用矩阵来加速递推运算

我们假设有一个矩阵 M×[FnFn+1]=[Fn+1Fn+2]M \times \begin{bmatrix}F_n\\F_{n+1}\end{bmatrix} = \begin{bmatrix}F_{n+1}\\F_{n+2}\end{bmatrix} ，根据斐波那契递推式 Fn+2=Fn+1+FnF_{n+2} = F_{n+1} + F_n 可以解出矩阵 M=[0111]M = \begin{bmatrix}0&1\\1&1\end{bmatrix}

我们定义初始矩阵为 ans=[11]ans = \begin{bmatrix}1\\1\end{bmatrix} base=[0111]base = \begin{bmatrix}0&1\\1&1\end{bmatrix}，运用 快速幂 的思想，每次运算时对 nn 右移一位，然后如果 nn&1=11 = 1 就将 ans=base×ansans = base \times ans ，结束后将 base=base2base = base^2 。使用 log⁡n\log n 的时间解决运算

参考答案 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
#include <bits/stdc++.h>
#define int long long
using namespace std;
const int M = 1e9+7;

int n;
int nans[3][3], nres[3][3];
int ans[3][3] = {{},{0,1,0},{0,1,0}};

void mul(int a[3][3], int b[3][3], int (&res)[3][3]) {
    res[1][1] = (a[1][1] * b[1][1] + a[1][2] * b[2][1]) % M;
    res[1][2] = (a[1][1] * b[1][2] + a[1][2] * b[2][2]) % M;
    res[2][1] = (a[2][1] * b[1][1] + a[2][2] * b[2][1]) % M;
    res[2][2] = (a[2][1] * b[1][2] + a[2][2] * b[2][2]) % M;
}

void quick(int n, int (&res)[3][3]) {
    while(n) {
        if(n&1) mul(res, ans, nans), memcpy(ans, nans, sizeof(ans));;
        mul(res, res, nres);
        memcpy(res, nres, sizeof(res));
        n >>= 1;
    }
}

signed main() {
    cin >> n;
    int res[3][3] = {{},{0,0,1},{0,1,1}};
    if (n >= 2) quick(n-2, res);
    cout << ans[2][1] % M;
    return 0;
}

如果我们将其打包为结构体，可以制作矩阵乘法的模板

矩阵快速幂算法模板 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
#include <bits/stdc++.h>
using namespace std;
const int N = 101, M = 1e9+7;

int n;
long long k;

struct mat {
    int c[N][N];
    mat operator* (const mat a) const {
        mat res;
        for(int i = 0;i < N;i++) {
            for(int j = 0;j < N;j++) {
                for(int l = 0;l < N;l++) {
                    res.c[i][j] = (res.c[i][j] + 1LL * c[i][l] * a.c[l][j]) % M;
                }
            }
        }
        return res;
    }
    mat() {memset(c, 0, sizeof(c));}
} ans, a;

void quick(long long k) {
    while(k) {
        if(k&1) ans = ans * a;
        a = a * a;
        k >>= 1;
    }
}

signed main(){
    cin >> n >> k;
    for(int i = 0;i < n;i++) ans.c[i][i] = 1;
    for(int i = 0;i < n;i++) for(int j = 0;j < n;j++) cin >> a.c[i][j];
    quick(k);
    for(int i = 0;i < n;i++) {
        for(int j = 0;j < n;j++) cout << ans.c[i][j] << " ";
        cout << endl;
    }
    return 0;
}