动态规划-矩阵加速

¶ 矩阵乘法

$$\left\{ \begin{aligned} 2x + y = 3 \\ 3x + 2y = 9 \end{aligned} \right.$$

可表示成

$$\begin{bmatrix} 2 & 1 \\ 3 & 2 \end{bmatrix} \times \begin{bmatrix} x \\ y \end{bmatrix}= \begin{bmatrix} 2x+y \\ 3x+2y \end{bmatrix}$$

$A$ 是 $m\times n$ 矩阵，$B$ 是 $n\times p$ 矩阵。

$A\times B = C$

$C$ 是 $m\times p$ 矩阵

$$C = (c_{i,j})\\ c_{i,j} =\sum_{k=1}^n a_{i,k}\times b_{k,j}$$

¶ 形象理解

$C =c_{i,j} = A$ 的第 $i$ 行逐个乘 $B$ 的第 $j$ 列

非常形象的解释

¶ 结合律证明

$A$ 是 $m\times n$ 矩阵，$B$ 是 $n\times p$ 矩阵，$C$ 是 $p \times q$ 的矩阵。

$$A\times B \times C = A \times (B \times C)\\ A \times B = D \\ D = (d_{i,j}) \\ d_{i,j} = \sum_{k=1}^{n} a_{i,k}\cdot b_{k,j} \\ D \times C = E \\ \begin{aligned} e_{i,j} &= \sum_{l=1}^p (d_{i,l}\cdot c_{l,j}) \\ &= \sum_{l=1}^p ((\sum_{k=1}^{n} a_{i,k}\cdot b_{k,l}) \times c_{l,j}) \\ &= \sum_{l=1}^p (\sum_{k=1}^{n} a_{i,k}\cdot b_{k,l}\cdot c_{l,j}) \\ \end{aligned} \\ B \times C = F \\ f_{i,j} = \sum_{l=1}^p b_{i,l}\cdot c_{l,j} \\ A \times F = G \\ \begin{aligned} g_{i,j} &= \sum_{k=1}^n(a_{i,k}\cdot f_{k,j}) \\ &= \sum_{k=1}^n(a_{i,k}\times (\sum_{l=1}^p b_{k,l}\cdot c_{l,j})) \\ &= \sum_{k=1}^n(\sum_{l=1}^p a_{i,k} \cdot b_{k,l}\cdot c_{l,j}) \end{aligned}$$

显然

$$\sum_{l=1}^p (\sum_{k=1}^{n} a_{i,k}\cdot b_{k,l}\cdot c_{l,j}) = \sum_{k=1}^n(\sum_{l=1}^p a_{i,k} \cdot b_{k,l}\cdot c_{l,j}) \\ e_{i,j} = g_{i,j} \\ E = G \\ A\times B\times C = A\times(B\times C)$$

¶ 快速幂

若 $2|p$ 则 $a^p = (a^2)^{p/2}$

若 $2|p+1$ 则 $a^p = a \cdot (a^2)^{(p-1)/2}$

将 $O(n)$ 变成 $O(\log_2n)$

¶ 例题

¶ 洛谷 P1939 【模板】矩阵加速（数列）

没什么好说的

由于每项只和前 $3$ 项有关，我们可以开一个 $1\times 3$ 的矩阵表示数列。

然后考虑一下转移。

$$\begin{bmatrix} a_1 & a_2 & a_3 \\ \end{bmatrix} \rightarrow \begin{bmatrix} a_2 & a_3 & a_4 \\ \end{bmatrix}= \begin{bmatrix} a_2 & a_3 & a_1+a_3 \\ \end{bmatrix}$$

用心感觉一下 转移矩阵就是（也可以解一个九元方程

$$\begin{bmatrix} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 1 \\ \end{bmatrix}$$

可以有

$$\begin{bmatrix} a_1 & a_2 & a_3 \end{bmatrix} \times \begin{bmatrix} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 1 \\ \end{bmatrix}= \begin{bmatrix} a_2 & a_3 & a_1+a_3 \\ \end{bmatrix}= \begin{bmatrix} a_2 & a_3 & a_4 \\ \end{bmatrix}$$

踩了洛谷g++的一个坑，sync_with_stdio(false)必须加在所有输入输出之前。

#include <iostream>
#include <cstring>
typedef long long LL;
using std::cin; using std::cout; using std::cerr; using std::endl; using std::min; using std::max;
const int MOD = 1000000007;
struct matrix{
    long long a[5][5], m, n;
    void init(const int &nM, const int &nN){
        std::memset(this->a, 0, sizeof(this->a));
        this->m = nM;
        this->n = nN;
    }
    matrix operator=(const matrix &s){
        this->init(s.m, s.n);
        for (int i = 0; i <= s.m; ++i)
            for (int j = 0; j <= s.n; ++j)
                this->a[i][j] = s.a[i][j];
        return *this;
    }
    matrix operator*(const matrix &s){
        // if (this->n != s.m)
            // cerr << "Matrixs Can't Time";
        matrix ans;
        ans.init(this->m, s.n);
        for (int i = 1; i <= this->m; ++i){
            for (int j = 1; j <= s.n; ++j){
                for (int k = 1; k <= this->n; ++k){
                    ans.a[i][j] += this->a[i][k]*s.a[k][j]%MOD;
                }
            }
        }
        return ans;
    }
}A, B, C;
std::ostream& operator<<(std::ostream &out, const matrix &s){
    for (int i = 1; i <= s.m; ++i){
        out << '[';
        for (int j = 1; j <= s.n-1; ++j){
            out << s.a[i][j] << ' ';
        }
        out << s.a[i][s.n];
        out << "]\n";
    }
    return out;
}
int main(){
    int T, n;
    std::ios::sync_with_stdio(false);
    cin >> T;
    while (T--){
        cin >> n;
        if (n <= 3){
            cout << 1 << endl;
            continue;
        }
        A.init(1,3);
        A.a[1][1] = A.a[1][2] = A.a[1][3] = 1;
        B.init(3,3);
        B.a[1][1] = B.a[1][2] = B.a[2][2] = B.a[2][3] = B.a[3][1] = 0;
        B.a[1][3] = B.a[2][1] = B.a[3][2] = B.a[3][3] = 1;
        n -= 4;
        C = B;
        while (n){
            if (n&1){
                C = C * B;
            }
            B = B * B;
            n = n >> 1;
        }
        A = A * C;
        cout << A.a[1][3]%MOD << endl;
    }
    return 0;
}

¶ 摆花

递归很好想

#include <iostream>
#include <cstring>
const int MOD = 1000000007;
using std::cin; using std::cout; using std::cerr; using std::endl; using std::min; using std::max;
int dp[2][101], a[101][101], n, m;
int main(){
    cin >> n >> m;
    dp[1][0] = 1;
    // cin.get();
    char c;
    for (register int i = 1; i <= m; ++i){
        dp[1][i] = 1;
        for (register int j = 1; j <= m; ++j){
            cin >> c;
            a[i][j] = c - '0';
        }
    }
    int cur = 0, prev = 1;
    for (register int i = 2; i <= n; ++i){
        for (register int j = 0; j <= m; ++j){
            dp[cur][j] = 0;
            for (register int k = 0; k <= m; ++k){
                if (a[j][k]) continue;
                dp[cur][j] += dp[prev][k];
                dp[cur][j] %= MOD;
            }
        }
        std::swap(cur, prev);
    }
    int ans = 0;
    for (int i = 0; i <= m; ++i){
        ans += dp[prev][i];
        ans %= MOD;
    }
    cout << ans;
    return 0;
}

可以用矩阵加速转移。

#include <iostream>
#include <cstring>
typedef long long LL;
const int MOD = 1000000007;
using std::cin; using std::cout; using std::cerr; using std::endl; using std::min; using std::max;
int n, m;
struct matrix{
    int a[110][110], m, n;
    void init(const int &nM, const int &nN){
        std::memset(this->a, 0, sizeof(this->a));
        this->m = nM;
        this->n = nN;
    }
    matrix operator=(const matrix &s){
        this->init(s.m, s.n);
        for (int i = 0; i <= s.m; ++i)
            for (int j = 0; j <= s.n; ++j)
                this->a[i][j] = s.a[i][j];
        return *this;
    }
    matrix operator*(const matrix &s){
        matrix ans;
        ans.init(this->m, s.n);
        for (int i = 1; i <= this->m; ++i){
            for (int j = 1; j <= s.n; ++j){
                for (int k = 1; k <= this->n; ++k){
                    ans.a[i][j] += (LL)this->a[i][k]%MOD*(LL)s.a[k][j]%MOD;
                    ans.a[i][j] %= MOD;
                }
            }
        }
        return ans;
    }
}A, B;
std::ostream& operator<<(std::ostream &out, const matrix &s){
    for (int i = 1; i <= s.m; ++i){
        out << '[';
        for (int j = 1; j <= s.n-1; ++j){
            out << s.a[i][j] << ' ';
        }
        out << s.a[i][s.n];
        out << "]\n";
    }
    return out;
}
int main(){
    cin >> n >> m;
    char c;
    A.init(m+1,m+1);
    for (register int i = 1; i <= m; ++i){
        A.a[i][m+1] = 1;
        for (register int j = 1; j <= m; ++j){
            cin >> c;
            A.a[i][j] = c - '0';
            A.a[i][j] ^= 1;
        }
    }
    for (register int i = 1; i <= m+1; ++i){
        A.a[m+1][i] = 1;
    }
    B.init(m+1,m+1);
    B = A;
    int p = n-2;
    while (p > 0){
        if (p%2!=0) B = B*A;
        A = A*A;
        p /= 2;
    }
    matrix ans;
    ans.init(1, m+1);
    for (int i = 1; i <= ans.n; ++i) ans.a[1][i] = 1;
    ans = ans * B;
    int sum = 0;
    for (int i = 1; i <= ans.n; ++i){
        sum += ans.a[1][i];
        sum %= MOD;
    }
    cout << sum;
    return 0;
}