矩阵快速幂
矩阵快速幂
1. 算法简介
矩阵快速幂是一种用于高效计算矩阵的 n 次幂的算法，其核心思想是基于 快速幂（前置芝士）。
其时间复杂度为 O(log n)，远快于朴素的 O(n)方法。
矩阵乘法满足结合律：(A×B)×C=A×(B×C) 这使得矩阵乘方可以用相同的“折半指数”的思路快速计算
矩阵快速幂常用于：
- 斐波那契数列求第 n项
- 动态规划转移优化
- 图论（如状态转移）
2. 基础公式
给定一个矩阵 A，我们要求出：
\[A^n = \underbrace{A \times A \times \cdots \times A}_{n\text{ 次}}
\]
矩阵乘法公式：
\[C_{i,j} = \sum_{k=1}^{r} A_{i,k} \cdot B_{k,j}
\]
矩阵乘法演示
我们来看一个矩阵乘法的具体计算过程：
设有两个 矩阵：
\[A =
\begin{bmatrix}
1 & 2 \\
3 & 4
\end{bmatrix}
\quad
B =
\begin{bmatrix}
5 & 6 \\
7 & 8
\end{bmatrix}
\]
它们的乘积：
\[C = A \times B =
\begin{bmatrix}
19 & 22 \\
43 & 50
\end{bmatrix}
\]
计算过程
\[每个元素 C_{i,j} 是 A 的第 i 行 与 B 的第 j 列 做点积得到的。
\]
第 1 行 × 第 1 列：
\[C_{1,1} = 1 \times 5 + 2 \times 7 = 5 + 14 = 19
\]
第 1 行 × 第 2 列：
\[C_{1,2} = 1 \times 6 + 2 \times 8 = 6 + 16 = 22
\]
第 2 行 × 第 1 列：
\[C_{2,1} = 3 \times 5 + 4 \times 7 = 15 + 28 = 43
\]
第 2 行 × 第 2 列：
\[C_{2,2} = 3 \times 6 + 4 \times 8 = 18 + 32 = 50
\]
最终结果矩阵：
\[A \times B =
\begin{bmatrix}
19 & 22 \\
43 & 50
\end{bmatrix}
\]
构造例题（斐波那契数列）
目标
我们要求的是斐波那契数列第 f(n) 项。
斐波那契定义：
\[f(n) = f(n-1) + f(n-2)
\]
我们希望使用 矩阵快速幂 来快速求解 f(n)，将时间复杂度降为 O(log n)。
1. 构造基本转移矩阵
从斐波那契递推式中，我们可以构造如下矩阵关系：
\[\begin{bmatrix}
f(n) \\
f(n-1)
\end{bmatrix}
=
\begin{bmatrix}
1 & 1 \\
1 & 0
\end{bmatrix}
\times
\begin{bmatrix}
f(n-1) \\
f(n-2)
\end{bmatrix}
\]
这就是核心的状态转移矩阵（记作 B）：
\[B =
\begin{bmatrix}
1 & 1 \\
1 & 0
\end{bmatrix}
\]
2. 迭代构造过程（推理）
根据上面的定义，我们令：
\[C(n) =
\begin{bmatrix}
f(n) \\
f(n-1)
\end{bmatrix}
\]
由于：
\[C(n) = B \times C(n-1)
\]
我们可以不断迭代推理：
\[\begin{align*}
C(n) &= B \cdot C(n-1) \\
&= B \cdot (B \cdot C(n-2)) \\
&= B^2 \cdot C(n-2) \\
&= B^3 \cdot C(n-3) \\
&\cdots \\
&= B^{n-1} \cdot C(1)
\end{align*}
\]
而初始向量为：
\[C(1) =
\begin{bmatrix}
f(1) \\
f(0)
\end{bmatrix}
=
\begin{bmatrix}
1 \\
0
\end{bmatrix}
\]
3. 最终公式总结
\[C(n) = B^{n-1} \cdot
\begin{bmatrix}
f(1) \\
f(0)
\end{bmatrix}
\]
因此，我们只需要快速计算矩阵 B 的 n - 1 次幂即可！
Code
void solve() {
i64 n;
cin >> n;
if (n <= 2) {
cout << 1 << '\n';
return;
}
Matrix T(2);
T.mat[0][0] = 1;
T.mat[0][1] = 1;
T.mat[1][0] = 1;
T.mat[1][1] = 0;
Matrix V(2);
V.mat[0][0] = 1;
V.mat[1][0] = 0;
V.mat[0][1] = 0;
V.mat[1][1] = 0;
Matrix ans = T.power(n - 1) * V;
cout << (ans.mat[0][0] + MOD) % MOD << '\n';
}
Template
const i64 MOD = 1e9 + 7;
class Matrix {
public:
vector<vector<i64>> mat;
Matrix(int N) { mat.resize(N, vector<i64>(N, 0)); }
// 输出矩阵
void print() {
for (i64 i = 0; i < mat.size(); i++) {
for (i64 j = 0; j < mat[i].size(); j++) {
cout << mat[i][j] << " ";
}
cout << endl;
}
}
static Matrix identity(int N) {
Matrix res(N);
for (i64 i = 0; i < N; i++) {
res.mat[i][i] = 1;
}
return res;
}
// 矩阵加法
Matrix operator+(const Matrix &other) const {
int N = mat.size();
Matrix res(N);
for (i64 i = 0; i < N; i++) {
for (i64 j = 0; j < N; j++) {
res.mat[i][j] = (mat[i][j] + other.mat[i][j] + MOD) % MOD;
}
}
return res;
}
// 矩阵减法
Matrix operator-(const Matrix &other) const {
int N = mat.size();
Matrix res(N);
for (i64 i = 0; i < N; i++) {
for (i64 j = 0; j < N; j++) {
res.mat[i][j] = (mat[i][j] - other.mat[i][j] + MOD) % MOD;
}
}
return res;
}
// 矩阵乘法
Matrix operator*(const Matrix &other) const {
int N = mat.size();
Matrix res(N);
for (i64 i = 0; i < N; i++) {
for (i64 j = 0; j < N; j++) {
res.mat[i][j] = 0;
for (i64 k = 0; k < N; k++) {
res.mat[i][j] =
(res.mat[i][j] + mat[i][k] * other.mat[k][j] + MOD) % MOD;
}
}
}
return res;
}
// 矩阵与数相乘
Matrix operator*(const i64 &x) const {
int N = mat.size();
Matrix res(N);
for (i64 i = 0; i < N; i++) {
for (i64 j = 0; j < N; j++) {
res.mat[i][j] = mat[i][j] * x % MOD;
}
}
return res;
}
Matrix power(i64 exp) const {
int N = mat.size();
Matrix res = identity(N);
Matrix base = *this;
while (exp) {
if (exp & 1) res = res * base;
base = base * base;
exp >>= 1;
}
return res;
}
// 矩阵转置
Matrix transpose() const {
int N = mat.size();
Matrix res(N);
for (i64 i = 0; i < N; i++) {
for (i64 j = 0; j < N; j++) {
res.mat[j][i] = mat[i][j];
}
}
return res;
}
bool invert(Matrix &out) const {
int N = mat.size();
Matrix res = *this;
out = identity(N);
for (i64 i = 0; i < N; i++) {
i64 j = i;
while (j < N && res.mat[j][i] == 0) j++;
if (j == N) return false; // 如果某列没有非零元素，矩阵不可逆
if (i != j) {
swap(res.mat[i], res.mat[j]);
swap(out.mat[i], out.mat[j]);
}
i64 inv = modInverse(res.mat[i][i], MOD);
for (i64 k = 0; k < N; k++) {
res.mat[i][k] = res.mat[i][k] * inv % MOD;
out.mat[i][k] = out.mat[i][k] * inv % MOD;
}
for (i64 j = 0; j < N; j++) {
if (j == i) continue;
i64 factor = res.mat[j][i];
for (i64 k = 0; k < N; k++) {
res.mat[j][k] =
(res.mat[j][k] - res.mat[i][k] * factor % MOD + MOD) % MOD;
out.mat[j][k] =
(out.mat[j][k] - out.mat[i][k] * factor % MOD + MOD) % MOD;
}
}
}
return true;
}
// 求逆元
i64 modInverse(i64 a, i64 m) const {
i64 x, y;
exgcd(a, m, x, y);
return (x % m + m) % m;
}
static i64 exgcd(i64 a, i64 b, i64 &x, i64 &y) {
if (b == 0) {
x = 1;
y = 0;
return a;
}
i64 g = exgcd(b, a % b, y, x);
y -= a / b * x;
return g;
}
};
浙公网安备 33010602011771号