DP与矩阵

前言

Sky

本人DP很差，但打ACM的队友DP很强，但是在有些比较难的题里面，推出的DP式子递推次数过多容易超时，这个时候就需要用矩阵快速幂来优化

一维数组，多层状态

假设递推是
$$
f_n = a_1 f_{n-1} + a_2 f_{n-2} + \cdots + a_k f_{n-k}
$$
我们可以构造一个 $k$ 维向量
$$
F_n =
\begin{bmatrix}
f_n \\
f_{n-1} \
f_{n-2} \\
\vdots \\
f_{n-k+1}
\end{bmatrix}
$$
那么递推式可以写成
$$
F_n = M \cdot F_{n-1}
$$
其中矩阵 $M$ 是
$$
\begin{bmatrix}
a_1 & a_2 & a_3 & \cdots & a_{k-1} & a_k \\
1 & 0 & 0 & \cdots & 0 & 0 \\
0 & 1 & 0 & \cdots & 0 & 0 \\
0 & 0 & 1 & \cdots & 0 & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
0 & 0 & 0 & \cdots & 1 & 0
\end{bmatrix}
$$
即
$$
F_n =
\begin{bmatrix}
a_1 f_{n-1} + a_2 f_{n-2} + \cdots + a_k f_{n-k} \\
f_{n-1} \
f_{n-2} \\
\vdots \\
f_{n-k+1}
\end{bmatrix}
$$

通项公式

假设我们有初始向量
$$
F_n =
\begin{bmatrix}
f_{k} \\
f_{k-1} \\
f_{k-2} \\
\vdots \\
f_{1}
\end{bmatrix}
$$
那么
$$
F_n = M^{n-k} \cdot F_k
$$
算出来第一行那个元素就是 $f_n$

e.g. 斐波那契数列

我们有斐波那契数列的递推式
$$
f_n = f_{n-1} + f_{n-2}
$$
那么可以构造矩阵
$$
M =
\begin{bmatrix}
1 & 1 \\
1 & 0
\end{bmatrix}
$$
还有向量
$$
F_2 =
\begin{bmatrix}
f_2 \\
f_1
\end{bmatrix}
$$
那么递推式子可以转换成
$$
F_n = M^{n-2} \cdot F_2
$$

Example Code C++

#include <bits/stdc++.h>
#define MAXN 101
#define LL long long
using namespace std;

const LL MOD = 1e9+7;

// re = i * j (Matrix)
// size(i) = n*m
// size(j) = m*k
void matrixTimes(LL i[MAXN][MAXN], LL j[MAXN][MAXN], LL re[MAXN][MAXN], int n, int m, int k)
{
    LL a[MAXN][MAXN] = {};
    for(int x = 1; x <= n; x++)
        for(int y = 1; y <= m; y++)
            for(int z = 1; z <= k; z++)
                a[x][y] = (a[x][y] + (i[x][z] * j[z][y]) % MOD) % MOD;
    for(int x = 1; x <= n; x++)
        for(int y = 1; y <= m; y++) re[x][y] = a[x][y];
}

void matrixFastPow(LL A[MAXN][MAXN], int n, LL k)
{
    LL I[MAXN][MAXN] = {};
    for(int x = 1; x <= n; x++) I[x][x] = 1;
    while(k > 0)
    {
        if(k & 1) matrixTimes(I, A, I, n, n, n);
        matrixTimes(A, A, A, n, n, n);
        k >>= 1;
    }
    for(int x = 1; x <= n; x++)
        for(int y = 1; y <= n; y++) A[x][y] = I[x][y];
}

LL Fib[MAXN][MAXN], M[MAXN][MAXN];

int main()
{
    M[1][1] = M[1][2] = M[2][1] = 1;
    Fib[1][1] = Fib[2][1] = 1;
    int n;    
    cin >> n;
    matrixFastPow(M, 2, n-2);
    LL re[MAXN][MAXN] = {};
    matrixTimes(M, Fib, re, 2, 2, 2);
    cout << re[1][1] << endl;
    return 0;
}

二维数组，单层状态

有人问：主播主播，你这个怎么只有一维数组递推啊，要是我有二维你不炸了吗

其实，二维也可以，不过就要稍微改一下

形如
$$
\begin{cases}
f_{i,1} = a_{1,1} f_{i-1,1} + a_{1,2} f_{i-1,2} + \cdots + a_{1,k} f_{i-1,k} \\
f_{i,2} = a_{2,1} f_{i-1,1} + a_{2,2} f_{i-1,2} + \cdots + a_{2,k} f_{i-1,k} \\
\ \ \vdots \\
f_{i,k} = a_{k,1} f_{i-1,1} + a_{k,2} f_{i-1,2} + \cdots + a_{k,k} f_{i-1,k}
\end{cases}
$$

可以尝试构造矩阵乘法
$$
\begin{bmatrix}
f_{i-1,1} & f_{i-1,2} & \cdots & f_{i-1,k}
\end{bmatrix}
\cdot
\begin{bmatrix}
a_{1,1} & a_{2,1} & \cdots & a_{k,1} \\
a_{1,2} & a_{2,2} & \cdots & a_{k,2} \\
\vdots & \vdots & \ddots & \vdots \\
a_{1,k} & a_{2,k} & \cdots & a_{k,k}
\end{bmatrix}=
\begin{bmatrix}
f_{i,1} & f_{i,2} & \cdots & f_{i,k}
\end{bmatrix}
$$

e.g. 洛谷P4910 帕秋莉的手环

像洛谷P4910这道题，设 $dp_{i,0}$ 表示第 $i$ 个选金的答案，$dp_{i,1}$ 表示第 $i$ 个选木的答案，则有dp式子：

$$
\begin{cases}
dp_{i,0} = dp_{i-1,0} + dp_{i-1,1} \\
dp_{i,1} = dp_{i-1,0}
\end{cases}
$$

当然因为是手环，所以还要分类讨论一下第一个是金还是木

这时候，我们通过观察可以发现矩阵乘法

$$
\begin{bmatrix}
dp_{i-1,0} & dp_{i-1,1}
\end{bmatrix}
\cdot
\begin{bmatrix}
1 & 1 \\
1 & 0
\end{bmatrix}=
\begin{bmatrix}
dp_{i,0} & dp_{i,1}
\end{bmatrix}
$$

那么就可以得到通项公式

$$
\begin{bmatrix}
dp_{1,0} & dp_{1,1}
\end{bmatrix}
\cdot
\begin{bmatrix}
1 & 1 \\
1 & 0
\end{bmatrix}^{n-1}=
\begin{bmatrix}
dp_{n,0} & dp_{n,1}
\end{bmatrix}
$$

Example Code C++

#include <bits/stdc++.h>
#define MAXN 101
#define INF 0x7fffffff
#define LL long long
using namespace std;

const LL MOD = 1e9+7;

// re = i * j (Matrix)
// size(i) = n*m
// size(j) = m*k
void matrixTimes(LL i[MAXN][MAXN], LL j[MAXN][MAXN], LL re[MAXN][MAXN], int n, int m, int k)
{
    LL a[MAXN][MAXN] = {};
    for(int x = 1; x <= n; x++)
        for(int y = 1; y <= m; y++)
            for(int z = 1; z <= k; z++)
                a[x][y] = (a[x][y] + (i[x][z] * j[z][y]) % MOD) % MOD;
    for(int x = 1; x <= n; x++)
        for(int y = 1; y <= m; y++) re[x][y] = a[x][y];
}

void matrixFastPow(LL A[MAXN][MAXN], int n, LL k)
{
    LL I[MAXN][MAXN] = {};
    for(int x = 1; x <= n; x++) I[x][x] = 1;
    while(k > 0)
    {
        if(k & 1) matrixTimes(I, A, I, n, n, n);
        matrixTimes(A, A, A, n, n, n);
        k >>= 1;
    }
    for(int x = 1; x <= n; x++)
        for(int y = 1; y <= n; y++) A[x][y] = I[x][y];
}

LL n, dp[MAXN][MAXN], M[MAXN][MAXN];

void cal(LL re[MAXN][MAXN])
{
    M[1][1] = M[1][2] = M[2][1] = 1;
    M[2][2] = 0;
    matrixFastPow(M, 2, n-1);
    matrixTimes(dp, M, re, 1, 2, 2);
}

int main()
{
    int T;
    cin >> T;
    while(T--)
    {
        LL ans = 0;
        cin >> n;
        dp[1][1] = 1;
        dp[1][2] = 0;
        LL re1[MAXN][MAXN] = {};
        cal(re1);
        ans = (ans + re1[1][1] + re1[1][2]) % MOD;
        dp[1][1] = 0;
        dp[1][2] = 1;
        LL re2[MAXN][MAXN] = {};
        cal(re2);
        ans = (ans + re2[1][1]) % MOD;
        cout << ans << endl;
    }
    return 0;
}