【矩阵性质】矩阵加速递推

本文介绍线性代数中一个非常重要的内容——矩阵（Matrix）的一个重要性质：矩阵加速递推

同时本文已经更新至：矩阵（Matrix）系统介绍篇

---

斐波那契数列（Fibonacci Sequence）大家应该都非常的熟悉了。在斐波那契数列当中，$F_1=F_2=1$，$F_i=F_{i-1}+F_{i-2}(i\ge 3)$。

如果有一道题目让你求斐波那契数列第 $n$ 项的值，最简单的方法莫过于直接递推了。但是如果 $n$ 的范围达到了 ${10}^{18}$ 级别，递推就不行了，稳 TLE。考虑矩阵加速递推。

设 $Fib(n)$ 表示一个 $1\times 2$ 的矩阵 $\left[\begin{array}{cc}{F}_n& F_{n-1}\end{array}\right]$。我们希望根据 $Fib(n-1)=\left[\begin{array}{cc}{F}_{n-1}& F_{n-2}\end{array}\right]$ 推出 $Fib(n)$。

试推导一个矩阵 $\text{base}$，使 $Fib(n-1)\times \text{base}=Fib(n)$，即 $\left[\begin{array}{cc}{F}_{n-1}& F_{n-2}\end{array}\right]\times \text{base}=\left[\begin{array}{cc}{F}_n& F_{n-1}\end{array}\right]$。

怎么推呢？因为 $F_n=F_{n-1}+F_{n-2}$，所以 $\text{base}$ 矩阵第一列应该是 $\left[\begin{array}{c}1\\ 1\end{array}\right]$，这样在进行矩阵乘法运算的时候才能令 $F_{n-1}$ 与 $F_{n-2}$ 相加，从而得出 $F_n$。同理，为了得出 $F_{n-1}$，矩阵 $\text{base}$ 的第二列应该为 $\left[\begin{array}{c}1\\ 0\end{array}\right]$。

综上所述：$\text{base}=\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]$ 原式化为 $\left[\begin{array}{cc}{F}_{n-1}& F_{n-2}\end{array}\right]\times \left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]=\left[\begin{array}{cc}{F}_n& F_{n-1}\end{array}\right]$

转化为代码，应该怎么求呢？

定义初始矩阵 $\text{ans}=\left[\begin{array}{cc}{F}_2& F_1\end{array}\right]=\left[\begin{array}{cc}1& 1\end{array}\right],\text{base}=\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]$。那么，$F_n$ 就等于 $\text{ans}\times {\text{base}}^{n-2}$ 这个矩阵的第一行第一列元素，也就是 $\left[\begin{array}{cc}1& 1\end{array}\right]\times {\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]}^{n-2}$ 的第一行第一列元素。

注意，矩阵乘法不满足交换律，所以一定不能写成 ${\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]}^{n-2}\times \left[\begin{array}{cc}1& 1\end{array}\right]$ 的第一行第一列元素。另外，对于 $n\le 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);
}

---

这是一个稍微复杂一些的例子。

$$\begin{gathered} f_1=f_2=0 \\ {f}_n=7f_{n-1}+6f_{n-2}+5n+4\times 3^n \end{gathered}$$

我们发现，$f_n$ 和 $f_{n-1},f_{n-2},n$ 有关，于是考虑构造一个矩阵描述状态。

但是发现如果矩阵仅有这三个元素 $\left[\begin{array}{ccc}{f}_n& f_{n-1}& n\end{array}\right]$ 是难以构造出转移方程的，因为乘方运算和 $+1$ 无法用矩阵描述。

于是考虑构造一个更大的矩阵。

$$\left[\begin{array}{ccccc}{f}_n& f_{n-1}& n& 3^n& 1\end{array}\right]$$

我们希望构造一个递推矩阵可以转移到

$$\left[\begin{array}{ccccc}{f}_{n+1}& f_n& n+1& 3^{n+1}& 1\end{array}\right]$$

转移矩阵即为

$$\left[\begin{array}{ccccc}7& 1& 0& 0& 0\\ 6& 0& 0& 0& 0\\ 5& 0& 1& 0& 0\\ 12& 0& 0& 3& 0\\ 5& 0& 1& 0& 1\end{array}\right]$$

---

有了上面的基础能很快做出 又见斐波那契

图源来自网络上一位博主

这个相比普通的斐波那契数列多了后面四项，看一眼数据范围，使用普通的 $O(n)$ 的算法肯定会超时

因此这里需要使用矩阵快速幂（斐波那契数列的项数n一旦过大，就要考虑快速幂，普通算法时间空间都开销太大）。

使用矩阵快速幂的一个关键问题就是矩阵递推式。

参考普通快速幂那一片博客最后面的那个扩展式，就可以得到下面这个递推式了：

img

然后通过计算等价替换可得出该矩阵A：

img

下面只需要把普通斐波那契数列的构造由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";
}