矩阵乘法应用
矩阵快速幂¶
题意：给定一个 $n\times n$ 的矩阵 $A$，求 $A^k$。矩阵每个元素对 $10^9+7$ 取模。
- $n\leq 10^2,k\leq 10^{12}$。
由于矩阵乘法满足结合律，因此可以使用快速幂优化。
例如：$A^{13}=A^{8}\cdot A^{4}\cdot A^{1}$。也就是利用倍增维护出矩阵 $A$ 的 $2$ 的次幂的结果即可加速。
我们类比着之前整数快速幂的写法去写即可。这是之前的写法：
ll ksm(ll a, ll b)
{
ll res = 1;
while (b)
{
if (b & 1)
res = res * a % mod;
a = a * a % mod;
b >>= 1;
}
return res;
}
那么在做矩阵快速幂时，首先初始化矩阵 res 为单位矩阵 $I$。
单位矩阵保证矩阵乘法结果不变。确保第一次乘法结果与原矩阵一致。
接下来利用重载运算符技巧实现矩阵快速幂即可。
时间复杂度：$O(n^3\log{k})$。
struct matrix
{
int mat[N][N];
matrix() // 构造函数
{
memset(mat, 0, sizeof(mat));
}
};
matrix operator * (const matrix &a, const matrix &b) // 矩阵乘法
{
matrix c;
for (int k = 1; k <= n; k++)
{
for (int i = 1; i <= n; i++)
{
for (int j = 1; j <= n; j++)
{
c.mat[i][j] = (c.mat[i][j] + 1ll * a.mat[i][k] * b.mat[k][j]) % mod;
}
}
}
return c;
}
matrix operator ^ (matrix &a, ll b)
{
matrix res;
for (int i = 1; i <= n; i++) res.mat[i][i] = 1; // 单位矩阵
while (b)
{
if (b & 1) res = res * a; // 注意不要写成 a * res
a = a * a;
b >>= 1;
}
return res;
}
矩阵乘法加速递推¶
斐波那契数列定义：
\[
\begin{cases}
F_1 = 1, \; F_2 = 1 \\
F_n = F_{n-1} + F_{n-2}, \quad n \ge 3
\end{cases}
\]
目标：快速求第 $n$ 项，当 $n$ 很大（比如 $10^{18}$）时，直接递推太慢，需要矩阵加速。
方法一：列向量法¶
1. 构造列向量¶
把两项打包成列向量：
\[
\mathbf{v}_i =
\begin{bmatrix} F_i \\ F_{i-1} \end{bmatrix}
\]
2. 构造递推矩阵¶
我们希望：
\[
A\cdot \mathbf{v}_i=\mathbf{v}_{i+1}
\]
因此可知 $A$ 是一个 $2\times 2$ 的矩阵：即就是
\[
\begin{bmatrix} a & b \\ c & d \end{bmatrix}\cdot \begin{bmatrix} F_i \\ F_{i-1} \end{bmatrix}=\begin{bmatrix} F_{i+1} \\ F_{i} \end{bmatrix}
\]
按照 列向量线性组合视角 ，可以得到：
\[
F_i \begin{bmatrix} a \\ c \end{bmatrix}+F_{i-1}\begin{bmatrix} c \\ d \end{bmatrix}=\begin{bmatrix} F_{i+1} \\ F_{i} \end{bmatrix}
\]
即得到方程组：
\[\begin{cases}
a\cdot F_i+c\cdot F_{i-1}=F_{i+1}\\
c\cdot F_i+d\cdot F_{i-1}=F_{i}
\end{cases}\]
易得矩阵 $A$ 是：
\[
A = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}
\]
3. 快速求第 $n$ 项¶
\[
\mathbf{v}_n = A^{n-2} \mathbf{v}_2, \quad
\mathbf{v}_2 = \begin{bmatrix} F_2 \\ F_1 \end{bmatrix} = \begin{bmatrix} 1 \\ 1 \end{bmatrix}.
\]
- 利用 矩阵快速幂（$O(\log n)$）即可得到 $F_n$。
- 注意当 $n>2$ 在使用矩阵快速幂。
参考代码
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
constexpr int mod = 1e9 + 7;
ll n;
struct matrix
{
int mat[3][3];
matrix() // 构造函数
{
memset(mat, 0, sizeof(mat));
}
};
matrix operator * (const matrix &a, const matrix &b) // 矩阵乘法
{
matrix c;
for (int k = 1; k <= 2; k++)
{
for (int i = 1; i <= 2; i++)
{
for (int j = 1; j <= 2; j++)
{
c.mat[i][j] = (c.mat[i][j] + 1ll * a.mat[i][k] * b.mat[k][j]) % mod;
}
}
}
return c;
}
matrix operator ^ (matrix &a, ll b)
{
matrix res;
for (int i = 1; i <= 2; i++) res.mat[i][i] = 1; // 单位矩阵
while (b)
{
if (b & 1) res = res * a; // 注意不要写成 a * res
a = a * a;
b >>= 1;
}
return res;
}
int main()
{
ios::sync_with_stdio(false), cin.tie(0);
cin >> n;
if (n <= 2)
{
cout << 1;
return 0;
}
matrix F;
F.mat[1][1] = F.mat[2][1] = 1; // 列向量
matrix a;
a.mat[1][1] = a.mat[1][2] = a.mat[2][1] = 1; // 转移矩阵
a = a ^ (n - 2);
a = a * F;
cout << a.mat[1][1];
return 0;
}
方法二：行向量法¶
1. 构造行向量¶
把两项打包成行向量：
\[
\mathbf{w}_i = \begin{bmatrix} F_{i} & F_{i-1} \end{bmatrix}
\]
2. 构造递推矩阵¶
我们希望：
\[
\mathbf{w}_i B=\mathbf{w}_{i+1}
\]
首先依然可以确定 $B$ 是一个 $2\times 2$ 的矩阵。即就是
\[
\begin{bmatrix} F_i & F_{i-1} \end{bmatrix} \cdot \begin{bmatrix} a & b \\ c & d \end{bmatrix}=\begin{bmatrix} F_{i+1} & F_{i} \end{bmatrix}
\]
按照 行向量线性组合视角 ，可以得到：
\[\begin{bmatrix}
a\cdot F_i+c\cdot F_{i-1} & b\cdot F_i+d\cdot F_{i-1}
\end{bmatrix}=\begin{bmatrix} F_{i+1} & F_{i} \end{bmatrix}
\]
因此矩阵 $B$ 是：
\[
B = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}.
\]
3. 快速求第 $n$ 项¶
\[
\mathbf{w}_n = \mathbf{w}_2 B^{n-2}, \quad
\mathbf{w}_2 = \begin{bmatrix} F_2 & F_1 \end{bmatrix} = \begin{bmatrix} 1 & 1 \end{bmatrix}.
\]
- 同样使用 矩阵快速幂，$O(\log n)$ 得到 $F_n$。
参考代码
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
constexpr int mod = 1e9 + 7;
ll n;
struct matrix
{
int mat[3][3];
matrix() // 构造函数
{
memset(mat, 0, sizeof(mat));
}
};
matrix operator * (const matrix &a, const matrix &b) // 矩阵乘法
{
matrix c;
for (int k = 1; k <= 2; k++)
{
for (int i = 1; i <= 2; i++)
{
for (int j = 1; j <= 2; j++)
{
c.mat[i][j] = (c.mat[i][j] + 1ll * a.mat[i][k] * b.mat[k][j]) % mod;
}
}
}
return c;
}
matrix operator ^ (matrix &a, ll b)
{
matrix res;
for (int i = 1; i <= 2; i++) res.mat[i][i] = 1; // 单位矩阵
while (b)
{
if (b & 1) res = res * a; // 注意不要写成 a * res
a = a * a;
b >>= 1;
}
return res;
}
int main()
{
ios::sync_with_stdio(false), cin.tie(0);
cin >> n;
if (n <= 2)
{
cout << 1;
return 0;
}
matrix F;
F.mat[1][1] = F.mat[1][2] = 1; // 行向量
matrix a;
a.mat[1][1] = a.mat[1][2] = a.mat[2][1] = 1; // 转移矩阵
a = a ^ (n - 2);
a = F * a;
cout << a.mat[1][1];
return 0;
}
总结对比¶
- 核心思想：把递推关系转化为矩阵乘法，把连续项打包，重复矩阵乘法即可快速得到任意项。
- 列向量和行向量方法本质相同，只是乘法方向和矩阵形式不同。
线段树维护矩阵乘法¶
例题一：[THUSC 2017] 大魔法师¶
题目中包含六种修改操作。核心观察是：
这些操作都可以转化为对列向量
\[
\begin{bmatrix} A \\ B \\ C \\ 1 \end{bmatrix}
\]
的线性变换，即等价于左乘某个 $4\times 4$ 矩阵。
各操作的矩阵表示
操作一：$A_i \gets A_i + B_i$
\[
\begin{bmatrix}
1 & 1 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1 \\
\end{bmatrix}
\begin{bmatrix} A \\ B \\ C \\ 1 \end{bmatrix}=
\begin{bmatrix} A+B \\ B \\ C \\ 1 \end{bmatrix}
\]
操作二、三：与操作一类似，略。
操作四：$A_i \gets A_i + v$
\[
\begin{bmatrix}
1 & 0 & 0 & v \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1 \\
\end{bmatrix}
\begin{bmatrix} A \\ B \\ C \\ 1 \end{bmatrix}=
\begin{bmatrix} A+v \\ B \\ C \\ 1 \end{bmatrix}
\]
操作五：$B_i \gets B_i \cdot v$
\[
\begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & v & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1 \\
\end{bmatrix}
\begin{bmatrix} A \\ B \\ C \\ 1 \end{bmatrix}=
\begin{bmatrix} A \\ B\cdot v \\ C \\ 1 \end{bmatrix}
\]
操作六：$C_i \gets v$
\[
\begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 0 & v \\
0 & 0 & 0 & 1
\end{bmatrix}
\begin{bmatrix} A \\ B \\ C \\ 1 \end{bmatrix}=
\begin{bmatrix} A \\ B \\ v \\ 1 \end{bmatrix}
\]
解题方法
由此可知，所有操作都可以抽象为矩阵乘法。因此可以在线段树中维护每个区间的矩阵标记，通过区间矩阵乘法实现批量修改。
注意要点：
- 矩阵乘法常数较大：$4^3=64$ 次运算，但由于矩阵中有大量 $0$ 元素，可以通过 特判/循环展开 降低常数。
- 矩阵不满足交换律：在线段树标记下传时，必须严格区分“谁在左、谁在右”，即保证操作顺序正确。
👉 这样，问题就转化为一个 线段树 + 矩阵乘法 的模板题。
例题二：[NOIP2022] 比赛¶
题目要求：给定 $q$ 次询问，计算：
\[
\sum_{i=l}^r \sum_{j=i}^r \max(a_i,\cdots,a_j)\cdot \max(b_i,\cdots,b_j)
\]
思路分析
首先，将所有询问离线，并按右端点 $r$ 排序。
考虑固定右端点 $r$ 时的情况：
定义
\[
x_i = \max(a_i,\cdots,a_r), \quad y_i = \max(b_i,\cdots,b_r)
\]
那么查询结果可转化为：
\[
\sum_{i=l}^r x_i \cdot y_i
\]
暴力解法（20分）
我们可以直接维护 $x,y,f$ 三个数组：
- $x[i],y[i]$：表示区间后缀最大值
- $f[i]$：维护 $\sum x_i \cdot y_i$ 的历史和
复杂度为 $O(n^2+nq)$，如下代码：
参考代码
#include <bits/stdc++.h>
#define int unsigned long long
using namespace std;
constexpr int N = 2.5e5 + 5;
int T, n, a[N], b[N], q;
signed main()
{
ios::sync_with_stdio(false), cin.tie(0);
cin >> T >> n;
for (int i = 1; i <= n; i++)
cin >> a[i];
for (int i = 1; i <= n; i++)
cin >> b[i];
cin >> q;
vector<vector<pair<int, int>>> que(n + 1);
for (int i = 1; i <= q; i++)
{
int l, r;
cin >> l >> r;
que[r].push_back({l, i});
}
vector<int> x(n + 1), y(n + 1), f(n + 1), ans(q + 1);
for (int r = 1; r <= n; r++)
{
for (int i = 1; i <= r; i++)
{
x[i] = max(x[i], a[r]);
y[i] = max(y[i], b[r]);
f[i] += x[i] * y[i];
}
for (auto [l, i] : que[r])
{
for (int j = l; j <= r; j++)
{
ans[i] += f[j];
}
}
}
for (int i = 1; i <= q; i++) cout << ans[i] << "\n";
return 0;
}
优化思路
观察发现，查询本质上是 区间求和，同时 $f$ 的更新是 在历史和上累加。
如果仅靠线段树维护区间和，难点在于：需要同步维护历史和、$x,y$ 的变化。
启发来自 单调栈：
- $a_r$ 只有在成为某个后缀最大值时才会更新 $x_i$；
- 同理，$b_r$ 也只有在成为某个后缀最大值时才会更新 $y_i$。
因此，更新过程可以抽象为 区间赋值。
考虑在线段树中维护以下五类信息：
ans：历史和- $\sum xy$
- $\sum x$
- $\sum y$
- 区间长度
len
矩阵维护转移
对区间的操作，可以用矩阵左乘统一维护：
- 区间赋值 $x=v$
\[
\begin{bmatrix}
1 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & v & 0\\
0 & 0 & 0 & 0 & v\\
0 & 0 & 0 & 1 & 0\\
0 & 0 & 0 & 0 & 1\\
\end{bmatrix}
\cdot
\begin{bmatrix}
ans \\ \sum xy \\ \sum x \\ \sum y \\ len
\end{bmatrix}
\]
- 区间赋值 $y=v$
\[
\begin{bmatrix}
1 & 0 & 0 & 0 & 0\\
0 & 0 & v & 0 & 0\\
0 & 0 & 1 & 0 & 0\\
0 & 0 & 0 & 0 & v\\
0 & 0 & 0 & 0 & 1\\
\end{bmatrix}
\cdot
\begin{bmatrix}
ans \\ \sum xy \\ \sum x \\ \sum y \\ len
\end{bmatrix}
\]
- 更新历史和
\[
\begin{bmatrix}
1 & 1 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0\\
0 & 0 & 1 & 0 & 0\\
0 & 0 & 0 & 1 & 0\\
0 & 0 & 0 & 0 & 1\\
\end{bmatrix}
\cdot
\begin{bmatrix}
ans \\ \sum xy \\ \sum x \\ \sum y \\ len
\end{bmatrix}
\]
这样，每次操作都可以在线段树中完成，查询时 $O(\log n)$。
虽然常数较大，但可通过 循环展开、内联优化 等方式降低常数。
参考代码
for (int r = 1; r <= n; r++)
{
while (top1 && seq[r].first > seq[stk1[top1]].first) top1--;
while (top2 && seq[r].second > seq[stk2[top2]].second) top2--;
matrix op = get(seq[r].first, 1);
modify(1, 1, n, stk1[top1] + 1, r, op);
op = get(seq[r].second, 2);
modify(1, 1, n, stk2[top2] + 1, r, op);
op = get(1, 3);
modify(1, 1, n, 1, r, op);
for (auto [l, i] : que[r]) {
ans[i] = query(1, 1, n, l, r).mat[0][0];
}
stk1[++top1] = r, stk2[++top2] = r;
}