矩阵快速幂加速递推、优化dp

前置知识

快速幂：就是普通快速幂，看完代码其实你会发现矩阵快速幂和普通快速幂基本上没有区别。
线性代数基础：主要就是矩阵乘法的运算，外加一些矩阵乘法的运算律，如结合律，这是满足快速幂的基本条件。
动态规划：没什么可说的。

原理介绍

矩阵快速幂的原理是将矩阵乘法和快速幂相结合。其核心思想是，当需要计算矩阵 $A$ 的 $n$ 次幂时，借鉴快速幂的思路，把指数 $n$ 转化为二进制，然后通过矩阵乘法来计算。
说白了，只要你把矩阵乘法重载一下 $\ast$，就和普通快速幂没有区别，原理也和普通快速幂类似， 具体代码如下。

代码分析

constexpr int M = 2;\\定义矩阵大小
using Matrix = std::array<std::array<int, M>, M>;\\定义矩阵类型

void init(Matrix &a) {\\初始化为单位矩阵
    for (int i = 0; i < M; i++) {
        a[i][i] = 1;
    }
}

Matrix operator * (const Matrix &a, const Matrix &b) {\\重载MAtrix的乘法运算符
    Matrix ans{};
    for (int i = 0; i < M; i++) {
        for (int j = 0; j < M; j++) {
            for (int k = 0; k < M; k++) {
                ans[i][j] = (ans[i][j] + a[i][k] * b[k][j] % mod) % mod;
            }
        }
    }
    return ans;
}

Matrix power(Matrix a, int b) {\\矩阵快速幂，和普通快速幂基本没有区别
    Matrix ans{};
    init(ans);
    while (b) {
        if (b & 1) {
            ans = ans * a;
        }
        a = a * a;
        b >>= 1;
    }
    return ans;
}

例题分析

1.斐波那契数列

题面很简单，输出斐波那契数列的第 $n(n<2^{63})$ 项对 $10^9+7$取模。
这里先给出递推矩阵
$$\left[\begin{array}{c}f(n)\\ f(n-1)\end{array}\right]={\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]}^{n-1}\times \left[\begin{array}{c}f(1)\\ f(0)\end{array}\right]$$
根据公式可以直接写出代码。

代码

void solve() {
    int n;
    std::cin >> n;
    Matrix a;
    a[0][0] = 1, a[0][1] = 1, a[1][0] = 1, a[1][1] = 0;
    Matrix ans = power(a, n - 1);
    std::cout << ans[0][0] << '\n';
}

下面，讲解一下这个递推矩阵是怎么来的。
首先先写出第 $i$ 项与第 $i-1$ 项的关系。
$$\left[\begin{array}{c}f(n)\\ ?\end{array}\right]=\left[\begin{array}{cc}?& ?\\ ?& ?\end{array}\right]\times \left[\begin{array}{c}f(n-1)\\ ?\end{array}\right]$$
根据斐波那契数列递推公式可以知道，$f(n)=f(n-1)+f(n-2)$。
所以可以推出右侧矩阵的 $?$ 应为 $f(n-2)$。
$$\left[\begin{array}{c}f(n)\\ ?\end{array}\right]=\left[\begin{array}{cc}?& ?\\ ?& ?\end{array}\right]\times \left[\begin{array}{c}f(n-1)\\ f(n-2)\end{array}\right]$$
所以左侧矩阵的 $?$ 为 $f(n-1)$ 。
$$\left[\begin{array}{c}f(n)\\ f(n-1)\end{array}\right]=\left[\begin{array}{cc}?& ?\\ ?& ?\end{array}\right]\times \left[\begin{array}{c}f(n-1)\\ f(n-2)\end{array}\right]$$
我们设递推矩阵为
$$\left[\begin{array}{cc}{x}_1& x_2\\ x_3& x_4\end{array}\right]$$
可以得到
$$\left[\begin{array}{c}f(n)\\ f(n-1)\end{array}\right]=\left[\begin{array}{cc}{x}_1& x_2\\ x_3& x_4\end{array}\right]\times \left[\begin{array}{c}f(n-1)\\ f(n-2)\end{array}\right]$$
根据矩阵乘法计算可以得到
$$\left[\begin{array}{c}f(n)\\ f(n-1)\end{array}\right]=\left[\begin{array}{c}{x}_1\times f(n-1)+x_2\times f(n-2)\\ x_3\times f(n-1)+x_4\times (f(n-2)\end{array}\right]$$
很容易得到$x_1=1,x_2=1,x_3=1,x_4=0$。
因而得到如下
$$\begin{gathered} \left[\begin{array}{c}f(n)\\ f(n-1)\end{array}\right]=\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]\times \left[\begin{array}{c}f(n-1)\\ f(n-2)\end{array}\right] \\ \left[\begin{array}{c}f(n-1)\\ f(n-2)\end{array}\right]=\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]\times \left[\begin{array}{c}f(n-2)\\ f(n-3)\end{array}\right] \\ ⋮ \\ \left[\begin{array}{c}f(2)\\ f(1)\end{array}\right]=\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]\times \left[\begin{array}{c}f(1)\\ f(0)\end{array}\right] \end{gathered}$$
所以递推矩阵为
$$\left[\begin{array}{c}f(n)\\ f(n-1)\end{array}\right]=\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]\times \left[\begin{array}{c}f(n-1)\\ f(n-2)\end{array}\right]=\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]\times \left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]\times \left[\begin{array}{c}f(n-2)\\ f(n-3)\end{array}\right]=\cdots ={\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]}^{n-1}\times \left[\begin{array}{c}f(1)\\ f(0)\end{array}\right]$$

2. 矩阵加速（数列）

已知一个数列 $a$，它满足：
a x = { 1 x ∈ { 1 , 2 , 3 } a x − 1 + a x − 3 x ≥ 4 a_x= \begin{cases} 1 & x \in\{1,2,3\}\\ a_{x-1}+a_{x-3} & x \geq 4 \end{cases} ax​={1ax−1​+ax−3​​x∈{1,2,3}x≥4​

求 $a$ 数列的第 $n$ 项对 $10^9+7$ 取余的值。

类似于上面斐波那契数列，依然设递推矩阵
$$\left[\begin{array}{c}f(n)\\ f(n-1)\\ f(n-2)\end{array}\right]=\left[\begin{array}{ccc}{x}_1& x_2& x_3\\ x_4& x_5& x_6\\ x_7& x_8& x_9\end{array}\right]\times \left[\begin{array}{c}f(n-1)\\ f(n-2)\\ f(n-3)\end{array}\right]$$
很容易求出
$$\left[\begin{array}{ccc}{x}_1& x_2& x_3\\ x_4& x_5& x_6\\ x_7& x_8& x_9\end{array}\right]=\left[\begin{array}{ccc}1& 0& 1\\ 1& 0& 0\\ 0& 1& 0\end{array}\right]$$

代码(这里换了一个模板)

void solve() {
    int n;
    cin >> n;
    ModMatrix a(3, 3);
    a[0][0] = 1, a[0][2] = 1, a[1][0] = 1, a[2][1] = 1;
    cout << a.pow(n - 1)[0][0] << '\n';
}

3.E2. Trails (Medium)

$m$ 个小屋，每个小屋 $i$ 有 $s_i$ 条短路、$l_i$ 条长路连湖边。
每日规则：当前小屋 $\to$ 湖 $\to$ 任意小屋（可回原屋），且两条路至少 $1$ 条是短路。从小屋 $1$ 出发，连续走 $n$ 天，求总路线组合数（对 $10^9+7$ 取模）。
首先需要先分析出转移方程。
$$f(i,j)=\sum _{k=1}^{m}f(i-1,k)\times (s_k\cdot s_j+s_k\cdot l_j+l_k\cdot s_j)$$
转移方程有了，就可以设出递推矩阵了，这里假设 $m=3$ 。
$$\left[\begin{array}{c}f(n,1)\\ f(n,2)\\ f(n,3)\end{array}\right]=\left[\begin{array}{ccc}{x}_{1,1}& x_{1,2}& x_{1,3}\\ x_{2,1}& x_{2,2}& x_{2,3}\\ x_{3,1}& x_{3,2}& x_{3,3}\end{array}\right]\times \left[\begin{array}{c}f(n-1,1)\\ f(n-1,2)\\ f(n-1,3)\end{array}\right]$$
剩下的就是解方程了
$$\{\begin{array}{l}f(n,1)=x_{1,1}\times f(n-1,1)+x_{1,2}\times f(n-1,2)+x_{1,3}\times f(n-1,3)\\ f(n,2)=x_{2,1}\times f(n-1,1)+x_{2,2}\times f(n-1,2)+x_{2,3}\times f(n-1,3)\\ f(n,2)=x_{3,1}\times f(n-1,1)+x_{3,2}\times f(n-1,2)+x_{3,3}\times f(n-1,3)\end{array}$$
不难可以看出 $x_{i,j}=s_i\cdot s_j+s_i\cdot l_j+l_i\cdot s_j$

代码

void solve() {
    int m, n;
    cin >> m >> n;
    vector<int> s(m);
    for (int i = 0; i < m; i++) {
        cin >> s[i];
    }
    vector<int> l(m);
    for (int i = 0; i < m; i++) {
        cin >> l[i];
    }
    ModMatrix c(m, m);
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < m; j++) {
            c[i][j] = s[i] * s[j] + s[i] * l[j] + l[i] * s[j];
        }
    }
    auto ans = c.pow(n);
    int res = 0;
    for (int i = 0; i < m; i++) {
        res = (res + ans[0][i]) % mod;
    }
    cout << res << '\n';
}