【矩阵性质】矩阵加速递推
本文介绍线性代数中一个非常重要的内容——矩阵（Matrix）的一个重要性质：矩阵加速递推
同时本文已经更新至：矩阵（Matrix）系统介绍篇
斐波那契数列（Fibonacci Sequence）大家应该都非常的熟悉了。在斐波那契数列当中，\(F_1 = F_2 = 1\)，\(F_i = F_{i - 1} + F_{i - 2}(i \geq 3)\)。
如果有一道题目让你求斐波那契数列第 \(n\) 项的值，最简单的方法莫过于直接递推了。但是如果 \(n\) 的范围达到了 \(10^{18}\) 级别，递推就不行了，稳 TLE。考虑矩阵加速递推。
设 \(Fib(n)\) 表示一个 \(1 \times 2\) 的矩阵 \(\left[ \begin{array}{ccc}F_n & F_{n-1} \end{array}\right]\)。我们希望根据 \(Fib(n-1)=\left[ \begin{array}{ccc}F_{n-1} & F_{n-2} \end{array}\right]\) 推出 \(Fib(n)\)。
试推导一个矩阵 \(\text{base}\)，使 \(Fib(n-1) \times \text{base} = Fib(n)\)，即 \(\left[\begin{array}{ccc}F_{n-1} & F_{n-2}\end{array}\right] \times \text{base} = \left[ \begin{array}{ccc}F_n & F_{n-1} \end{array}\right]\)。
怎么推呢？因为 \(F_n=F_{n-1}+F_{n-2}\)，所以 \(\text{base}\) 矩阵第一列应该是 \(\left[\begin{array}{ccc} 1 \\ 1 \end{array}\right]\)，这样在进行矩阵乘法运算的时候才能令 \(F_{n-1}\) 与 \(F_{n-2}\) 相加，从而得出 \(F_n\)。同理，为了得出 \(F_{n-1}\)，矩阵 \(\text{base}\) 的第二列应该为 \(\left[\begin{array}{ccc} 1 \\ 0 \end{array}\right]\)。
综上所述：\(\text{base} = \left[\begin{array}{ccc} 1 & 1 \\ 1 & 0 \end{array}\right]\) 原式化为 \(\left[\begin{array}{ccc}F_{n-1} & F_{n-2}\end{array}\right] \times \left[\begin{array}{ccc} 1 & 1 \\ 1 & 0 \end{array}\right] = \left[ \begin{array}{ccc}F_n & F_{n-1} \end{array}\right]\)
转化为代码，应该怎么求呢？
定义初始矩阵 \(\text{ans} = \left[\begin{array}{ccc}F_2 & F_1\end{array}\right] = \left[\begin{array}{ccc}1 & 1\end{array}\right], \text{base} = \left[\begin{array}{ccc} 1 & 1 \\ 1 & 0 \end{array}\right]\)。那么，\(F_n\) 就等于 \(\text{ans} \times \text{base}^{n-2}\) 这个矩阵的第一行第一列元素，也就是 \(\left[\begin{array}{ccc}1 & 1\end{array}\right] \times \left[\begin{array}{ccc} 1 & 1 \\ 1 & 0 \end{array}\right]^{n-2}\) 的第一行第一列元素。
注意，矩阵乘法不满足交换律，所以一定不能写成 \(\left[\begin{array}{ccc} 1 & 1 \\ 1 & 0 \end{array}\right]^{n-2} \times \left[\begin{array}{ccc}1 & 1\end{array}\right]\) 的第一行第一列元素。另外，对于 \(n \leq 2\) 的情况，直接输出 \(1\) 即可，不需要执行矩阵快速幂。
为什么要乘上 \(\text{base}\) 矩阵的 \(n-2\) 次方而不是 \(n\) 次方呢？因为 \(F_1, F_2\) 是不需要进行矩阵乘法就能求的。也就是说，如果只进行一次乘法，就已经求出 \(F_3\) 了。如果还不是很理解为什么幂是 \(n-2\)，建议手算一下。
下面是求斐波那契数列第 \(n\) 项对 \(10^9+7\) 取模的示例代码（核心部分）。
const int mod = 1000000007;
struct Matrix {
int a[3][3];
Matrix() { memset(a, 0, sizeof a); }
Matrix operator*(const Matrix &b) const {
Matrix res;
for (int i = 1; i <= 2; ++i)
for (int j = 1; j <= 2; ++j)
for (int k = 1; k <= 2; ++k)
res.a[i][j] = (res.a[i][j] + a[i][k] * b.a[k][j]) % mod;
return res;
}
} ans, base;
void init() {
base.a[1][1] = base.a[1][2] = base.a[2][1] = 1;
ans.a[1][1] = ans.a[1][2] = 1;
}
void qpow(int b) {
while (b) {
if (b & 1) ans = ans * base;
base = base * base;
b >>= 1;
}
}
int main() {
int n = read();
if (n <= 2) return puts("1"), 0;
init();
qpow(n - 2);
println(ans.a[1][1] % mod);
}
这是一个稍微复杂一些的例子。
我们发现，\(f_n\) 和 \(f_{n-1}, f_{n-2}, n\) 有关，于是考虑构造一个矩阵描述状态。
但是发现如果矩阵仅有这三个元素 \(\begin{bmatrix}f_n& f_{n-1}& n\end{bmatrix}\) 是难以构造出转移方程的，因为乘方运算和 \(+1\) 无法用矩阵描述。
于是考虑构造一个更大的矩阵。
我们希望构造一个递推矩阵可以转移到
转移矩阵即为
有了上面的基础能很快做出 又见斐波那契
图源来自网络上一位博主
这个相比普通的斐波那契数列多了后面四项，看一眼数据范围，使用普通的 \(O(n)\) 的算法肯定会超时
因此这里需要使用矩阵快速幂（斐波那契数列的项数n一旦过大，就要考虑快速幂，普通算法时间空间都开销太大）。
使用矩阵快速幂的一个关键问题就是矩阵递推式。
参考普通快速幂那一片博客最后面的那个扩展式，就可以得到下面这个递推式了：
然后通过计算等价替换可得出该矩阵A：
下面只需要把普通斐波那契数列的构造由2*2的矩阵换为6*6的即可。
using ll = long long;
const int mod = 1e9 + 7;
const int N = 6;
ll tmp[N][N], res[N][N];
ll n;
void mul(ll a[][N], ll b[][N]) {
memset(tmp, 0, sizeof(tmp));
for (int i = 0; i < N; ++i)
for (int j = 0; j < N; ++j)
for (int k = 0; k < N; ++k) {
tmp[i][j] += a[i][k] * b[k][j] % mod;
if (tmp[i][j] >= mod) tmp[i][j] -= mod;
}
for (int i = 0; i < N; ++i)
for (int j = 0; j < N; ++j)a[i][j] = tmp[i][j];
}
void pow(ll a[][N]) {
memset(res, 0, sizeof(res));
for (int i = 0; i < N; ++i)res[i][i] = 1;
for (; n; n >>= 1) {
if (n & 1)mul(res, a);
mul(a, a);
}
}
void solve() {
cin >> n, n--;
ll ans[N][N] = { 1, 1, 1, 1, 1, 1,
1, 0, 0, 0, 0, 0,
0, 0, 1, 3, 3, 1,
0, 0, 0, 1, 2, 1,
0, 0, 0, 0, 1, 1,
0, 0, 0, 0, 0, 1
};
pow(ans);
ll sum = 0;
ll x[N] = {1, 0, 8, 4, 2, 1};
for (int i = 0; i < N; ++i) {
sum += (res[0][i] * x[i]) % mod;
if (sum >= mod) sum -= mod;
}
cout << sum << "\n";
}
