矩阵快速幂#

矩阵快速幂是算法竞赛中一个非常实用的工具，它的核心价值在于：把 $O(n)$ 的线性递推优化到 $O(\log n)$ ，从而解决 n 极大（比如 $10^{18}$ 量级）的问题。

常见用途有以下几种：

快速求解线性递推数列

图论中的路径计数

状态压缩DP的加速

前置线性代数知识#

矩阵乘法#

NOTE
注意矩阵快速幂中的矩阵一定是一个方阵，这样才有“幂”的概念。
下面的代码是通用的矩阵乘法形式，在矩阵快速幂中有 m == n。

1
 for (int i = 1; i <= n; i++)
2
        for (int j = 1; j <= m; j++)
3
            for (int k = 1; k <= n; k++)
4
                c[i][j] += a[i][k] * b[k][j];

单位矩阵#

主对角线全是 1，其他地方全是 0 的矩阵是单位矩阵，它乘任意矩阵的结果都等于任意矩阵本身。

矩阵快速幂#

由于矩阵乘法满足结合律，一个单位矩阵用快速幂乘上 n 次就是矩阵快速幂了。

特殊地，定义 $A^0$ 为单位矩阵 $I = \begin{bmatrix} 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{bmatrix}$ 。

Code#

1
// Problem: Luogu P3390
2
// Contest: Luogu
3
// URL: https://www.luogu.com.cn/problem/P3390
4
// Time: 2026-05-08 12:56:14
5
#include <bits/stdc++.h>
6
using namespace std;
7

8
// clang-format off
9
#define endl '\n'
10
#define all(x) (x).begin(), (x).end()
11
#define fastio() ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
12
// clang-format on
13

14
using ll = long long;
15
using ull = unsigned long long;
16
using pii = pair<int, int>;
17
using pdd = pair<double, double>;
18
using pll = pair<long long, long long>;
19

20
const int dx[] = {-1, 0, 1, 0, -1, 1, 1, -1};
21
const int dy[] = {0, 1, 0, -1, 1, 1, -1, -1};
22
const int inf = 0x3f3f3f3f;
23
const int N = 101;
24
const ll mod = 1e9 + 7;
25

26
auto matrix_mul(const vector<vector<ll>> &a, const vector<vector<ll>> &b)
27
{
28
    int n = (int) a.size() - 1;
29
    vector<vector<ll>> res(n + 1, vector<ll>(n + 1, 0));
30

31
    for (int i = 1; i <= n; i++)
32
        for (int j = 1; j <= n; j++)
33
            for (int k = 1; k <= n; k++)
34
                res[i][j] = (res[i][j] + a[i][k] * b[k][j]) % mod;
35

36
    return res;
37
}
38

39
auto matrix_power(vector<vector<ll>> m, ll k)
40
{
41
    int n = m.size() - 1;
42

43
    vector<vector<ll>> res(n + 1, vector<ll>(n + 1, 0)); // 注意一开始res是单位矩阵
44

45
    for (int i = 1; i <= n; i++)
46
        res[i][i] = 1;
47

48
    while (k) // 和快速幂相同
49
    {
50
        if (k & 1)
51
            res = matrix_mul(res, m);
52

53
        m = matrix_mul(m, m);
54

55
        k >>= 1;
56
    }
57

58
    return res;
59
}
60

61
void solve()
62
{
63
    int n;
64
    ll k;
65
    cin >> n >> k;
66

67
    vector<vector<ll>> a(n + 1, vector<ll>(n + 1));
68

69
    for (int i = 1; i <= n; i++)
70
        for (int j = 1; j <= n; j++)
71
            cin >> a[i][j];
72

73
    auto res = matrix_power(a, k);
74

75
    for (int i = 1; i <= n; i++)
76
        for (int j = 1; j <= n; j++)
77
            cout << res[i][j] << " \n"[j == n];
78
}
79

80
int main()
81
{
82
    fastio();
83

84
    int T = 1;
85
    // cin >> T;
86

87
    while (T--)
88
        solve();
89

90
    return 0;
91
}

矩阵加速#

了解了矩阵快速幂的基本算法以后，我们可以在很多情况下利用矩阵快速幂实现矩阵加速。

矩阵加速只适用于固定关系的递推式，也就是说不适用于存在条件转移的情况。

固定关系的一维一阶递推表达式，可以用基本的乘法快速幂解决，时间复杂度 $O(\log n)$ ，可以推广到一维 $k$ 阶。

例如有一个递推表达式为 $f(n)=3 \times f(n-1)$ ，我们在求每一项时只要知道前一项的值就可以了，这样的表达式称作一维一阶的递推式。

固定关系的一维 $k$ 阶递推表达式，可以用矩阵快速幂解决，时间复杂度 $O(\log n \times k^3)$ 。

例如，Fibonacci 数列的递推表达式是 $f(n)=f(n-1)+f(n-2)$ ，这就是一个一维二阶的递推式。一维二阶的递推式一定存在一个 2x2 的关系矩阵。

使用矩阵加速的难点在于构造初始的关系矩阵（也就是“单位矩阵”），关系矩阵的第1列直接由递推表达式的系数决定，剩下的项可以通过前面的初始项代入求解。

矩阵加速模板#

将矩阵乘法和矩阵快速幂封装在一个结构体里，方便调用模板。

使用重载运算符的方式进行矩阵乘法，下面进行一些代码实现的说明：

当写下 A * B 时，实际上是 A.operator*(B)，左参数会通过一个隐藏的指针 this 进行隐式调用，可以使用 this -> mat 来访问。右参数是显式传入的。
矩阵乘法中循环顺序为 i, k, j 可以加快计算速度。
assert()用于断言表达式，如果表达式错误，程序会输出错误信息并终止，否则正常运行。详见：assert - cppreference.cn - C++参考手册

1
struct Matrix
2
{
3
    int r, c;
4
    vector<vector<ll>> mat;
5

6
    // 构造函数：初始化为 r 行 c 列的零矩阵
7
    Matrix(int r, int c) : r(r), c(c), mat(r, vector<ll>(c, 0)){}
8

9
    // 重载乘法运算符
10
    Matrix operator*(const Matrix &other) const
11
    {
12
        // 结果矩阵的行数等于左矩阵的行数，列数等于右矩阵的列数
13
        Matrix res(r, other.c);
14

15
        for (int i = 0; i < r; i++)
16
        {
17
            for (int k = 0; k < c; k++)
18
            {
19
                // 如果乘数为 0，直接跳过内层循环
20
                if (mat[i][k] == 0)
21
                    continue;
22
                for (int j = 0; j < other.c; j++)
23
                    res.mat[i][j] = (res.mat[i][j] + mat[i][k] * other.mat[k][j]) % mod;
24
            }
25
        }
26
        return res;
27
    }
28

29
    Matrix power(ll k) const // 矩阵快速幂
30
    {
31
        assert(r == c);
32
        Matrix res(r, c);
33

34
        for (int i = 0; i < r; i++) // 初始化单位矩阵
35
            res.mat[i][i] = 1;
36

37
        Matrix base = *this;
38
        while (k > 0)
39
        {
40
            if (k & 1)
41
                res = res * base;
42
            base = base * base;
43
            k >>= 1;
44
        }
45
        return res;
46
    }
47
};

Fibonacci#

例题链接：P1962 斐波那契数列 - 洛谷

我们首先假设关系矩阵是 $\begin{bmatrix} a & c \\ b & d \end{bmatrix}$ ，由于关系式 $f(n)=f(n-1)+f(n-2)$ ，等号右侧的两个数系数都是 1，直接确定 a = 1, b = 1。由于初始矩阵是 $\begin{bmatrix} 1 & 0 \end{bmatrix}$ ，也就是 Fibonacci 数列的第 1 项和第 0 项，我们已知它乘上关系矩阵以后，一定可以得到 Fibonacci 数列的第 2 项和第 1 项，因此直接模拟一次矩阵乘法求出 c 和 d。

最终的矩阵递推式： $\begin{bmatrix} f_i \\ f_{i-1} \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \cdot \begin{bmatrix} f_{i-1} \\ f_{i-2} \end{bmatrix}$ ，注意特判边界情况：求前两项直接输出 1 即可。

Code#

1
// Problem: Luogu P1962
2
// Contest: Luogu
3
// URL: https://www.luogu.com.cn/problem/P1962
4
// Time: 2026-05-08 13:36:52
5
#include <bits/stdc++.h>
6
using namespace std;
7

8
// clang-format off
9
#define endl '\n'
10
#define all(x) (x).begin(), (x).end()
11
#define fastio() ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
12
// clang-format on
13

14
using ll = long long;
15
using ull = unsigned long long;
16
using pii = pair<int, int>;
17
using pdd = pair<double, double>;
18
using pll = pair<long long, long long>;
19

20
const int dx[] = {-1, 0, 1, 0, -1, 1, 1, -1};
21
const int dy[] = {0, 1, 0, -1, 1, 1, -1, -1};
22
const int inf = 0x3f3f3f3f;
23
const int N = 0;
24
const ll mod = 1e9 + 7;
25

26
struct Matrix
27
{
28
    int r, c;
29
    vector<vector<ll>> mat;
30

31
    Matrix(int r, int c) : r(r), c(c), mat(r, vector<ll>(c, 0))
32
    {
33
    }
34

35
    Matrix operator*(const Matrix &other) const
36
    {
37
        Matrix res(r, other.c);
38

39
        for (int i = 0; i < r; i++)
40
        {
41
            for (int k = 0; k < c; k++)
42
            {
43
                if (mat[i][k] == 0)
44
                    continue;
45
                for (int j = 0; j < other.c; j++)
46
                {
47
                    res.mat[i][j] = (res.mat[i][j] + mat[i][k] * other.mat[k][j]) % mod;
48
                }
49
            }
50
        }
51
        return res;
52
    }
53

54
    Matrix power(ll k) const
55
    {
56
        assert(r == c);
57
        Matrix res(r, c);
58

59
        for (int i = 0; i < r; i++)
60
            res.mat[i][i] = 1;
61

62
        Matrix base = *this;
63
        while (k > 0)
64
        {
65
            if (k & 1)
66
                res = res * base;
67
            base = base * base;
68
            k >>= 1;
69
        }
70
        return res;
71
    }
72
};
73

74
void solve()
75
{
76
    ll n;
77
    cin >> n;
78

79
    if (n <= 2)
80
    {
81
        cout << 1 << endl;
82
        return;
83
    }
84

85
    // 构造转移矩阵（关系矩阵）
86
    Matrix base(2, 2);
87
    base.mat[0][0] = 1;
88
    base.mat[0][1] = 1;
89
    base.mat[1][0] = 1;
90
    base.mat[1][1] = 0;
91

92
    Matrix res = base.power(n - 1);
93

94
    cout << res.mat[0][0] << endl;
95
}
96

97
int main()
98
{
99
    fastio();
100

101
    int T = 1;
102
    // cin >> T;
103

104
    while (T--)
105
        solve();
106

107
    return 0;
108
}

矩阵加速数列2#

洛谷 - P1939 矩阵加速（数列）

构造初始矩阵更通用的方法就是，当明确目标矩阵以后，像这样推一下递推式。

Code#

1
// Problem: Luogu P1939
2
// Contest: Luogu
3
// URL: https://www.luogu.com.cn/problem/P1939
4
// Time: 2026-05-12 15:59:21
5
#include <bits/stdc++.h>
6
using namespace std;
7

8
// clang-format off
9
#define endl '\n'
10
#define all(x) (x).begin(), (x).end()
11
#define fastio() ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
12
// clang-format on
13

14
using ll = long long;
15
using ull = unsigned long long;
16
using pii = pair<int, int>;
17
using pdd = pair<double, double>;
18
using pll = pair<long long, long long>;
19
using i128 = __int128;
20

21
const int dx[] = {-1, 0, 1, 0, -1, 1, 1, -1};
22
const int dy[] = {0, 1, 0, -1, 1, 1, -1, -1};
23
const int inf = 0x3f3f3f3f;
24
const int N = 0;
25
const ll mod = 1e9 + 7;
26

27
struct Matrix
28
{
29
    int r, c;
30
    vector<vector<ll>> mat;
31

32
    Matrix(int r, int c) : r(r), c(c), mat(r, vector<ll>(c, 0))
33
    {
34
    }
35

36
    Matrix operator*(const Matrix &other) const
37
    {
38
        Matrix res(r, other.c);
39

40
        for (int i = 0; i < r; i++)
41
        {
42
            for (int k = 0; k < c; k++)
43
            {
44
                if (mat[i][k] == 0)
45
                    continue;
46
                for (int j = 0; j < other.c; j++)
47
                {
48
                    res.mat[i][j] = (res.mat[i][j] + mat[i][k] * other.mat[k][j]) % mod;
49
                }
50
            }
51
        }
52
        return res;
53
    }
54

55
    Matrix power(ll k) const
56
    {
57
        assert(r == c);
58
        Matrix res(r, c);
59

60
        for (int i = 0; i < r; i++)
61
            res.mat[i][i] = 1;
62

63
        Matrix base = *this;
64
        while (k > 0)
65
        {
66
            if (k & 1)
67
                res = res * base;
68
            base = base * base;
69
            k >>= 1;
70
        }
71
        return res;
72
    }
73
};
74

75
int main()
76
{
77
    fastio();
78

79
    int T = 1;
80
    cin >> T;
81

82
    Matrix base(3, 3);
83
    base.mat[0][0] = 1, base.mat[0][1] = 0, base.mat[0][2] = 1;
84
    base.mat[1][0] = 1, base.mat[1][1] = 0, base.mat[1][2] = 0;
85
    base.mat[2][0] = 0, base.mat[2][1] = 1, base.mat[2][2] = 0;
86

87
    while (T--)
88
    {
89
        ll n;
90
        cin >> n;
91

92
        if (n <= 3)
93
        {
94
            cout << 1 << endl;
95
            continue;
96
        }
97

98
        Matrix res = base.power(n);
99
        cout << res.mat[1][0] << endl;
100
    }
101

102
    return 0;
103
}