矩阵优化dp

都快csps了，还什么都不会的菜鱼（我估计着马上就可以改了这句话了，成了都快noip了）

矩阵

我们要用矩阵优化dp，首先要知道矩阵是个什么东西（感觉其实可以不用知道）。

矩阵的很多定义啥的都可以选择去oi-wiki上去进行学习。很简单的一堆定义。读者自学不难，这里就不多赘述。

矩阵加法

就是将对应 \(A\) 和 \(B\) 两个矩阵对应位置相加。然后得到 \(C\) 这个矩阵。

矩阵乘法

这里大概讲一下矩阵乘法的原理：（突然学了，这里说的好像是内积，外积没见过用（逃）

我们会由一个 \(a \times b\) 的矩阵和一个 \(b \times c\) 的矩阵得到一个 \(a \times c\) 的矩阵。这个时候我们是将左边两个黄色的矩阵每一项相乘然后相加得到了右边矩阵的黄色位置，其他位置亦然，读者自想不难。

一个重要的性质：矩阵乘法是符合结合律的，这也是为什么我们可以用快速幂去优化dp

矩阵求逆

我们有时候会用到一个矩阵的逆矩阵，这个时候我们就要去思考如何计算出来这个逆矩阵。

我们通常是在我们的矩阵右边拼上一个单位矩阵，然后用高斯消元的方法，将左边的原矩阵，变成一个单位矩阵，此时我们右边得到的就是我们的原矩阵的逆矩阵。

矩阵除法

我们将我们要除以的矩阵转换为它的逆矩阵（类比逆元），这个时候再对他进行矩阵乘法。

矩阵求最大最小（这个名字是我自己取的）（貌似已经名算有主，叫做广义矩阵乘法）

感觉矩阵里面可能并没有这个操作，但是我们在优化dp的时候难免不会少了这种求最大最小的动规形式，所以就单独拎出来说一下，也就当作矩阵的一个基本运算吧（辉算？还是carp算）找到题目了再upd（逃

广义矩阵乘法的 \(A \times B = C\)为：

\[C_{i, j} = max_{k=1}^{n}(A_{i, k} + B_{k, j}) \]

广义矩阵乘法是具有结合律的，所以也可以用矩阵快速幂来进行优化，

先放个板子

之后再根据例题看怎么应用。

code

struct Matrix
{   
    ll a[N][N];

    inline void init() { memset(a, 0, sizeof a); }

    Matrix operator +(const Matrix &x) const 
    {
        Matrix res;

        fos(i, 0, N)
            fos(j, 0, N)
                res.a[i][j] = (a[i][j] + x.a[i][j]) % P;
        return res;
    }

    Matrix operator -(const Matrix &x) const
    {
        Matrix res;

        fos(i, 0, N - 1)
            fos(j, 0, N - 1)
                res.a[i][j] = (a[i][j] - x.a[i][j]) % P;

        return res;
    }

    Matrix operator *(const Matrix &x) const
    {
        Matrix res;
        fos(i, 0, N - 1)
            fos(k, 0, N - 1)
                fos(j, 0, N - 1) 
                    res.a[i][j] = (res.a[i][j] + (a[i][k] * x.a[k][j]) % P) % P;
        
        return res;
    }

    Matrix operator ^(ll k) const
    {
        Matrix ans, base;
        fos(i, 0, N - 1) res.a[i][i] = 1;
        fos(i, 0, N - 1)
            fos(j, 0,  N - 1) base.a[i][j] = a[i][j] % P;
        while(k)
        {
            if(k & 1) res = res * base;
            base = base * base; 
            k >>= 1;
        } 

        return res;
    }
} ans, base;

矩阵优化dp的适用范围

基本上就是对于一个线性递推式，我们构造一个答案矩阵，再构造一个转移矩阵，每次将这个转移矩阵进行快速幂，然后我们用答案矩阵去乘转移矩阵，具体就是一个快速幂的实现。

例题

举个很板的栗子。
斐波那契数列

我们都知道斐波那契数列里面有这个递推式

\[f_i = f_{i -1} + f _ {i - 2} \]

这个时候我们就可以用上面矩阵乘法的思路来进行求解：
step 1 ： 构造一个答案矩阵 \(\begin{bmatrix} f_{1} & f_{2} & f_{3}\end{bmatrix}\)

step 2 ： 我们还需要构造一个转移矩阵，这个转移矩阵就是依据我们的上面的线性递推式得出来的。 \(\begin{bmatrix} 0 & 0 & 0 \\ 1 & 0 & 1 \\ 0 & 1 & 1 \end{bmatrix}\)

step 3 ： 接下来就是进行一个快速幂的计算。进行多少次，取决于你希望让这个答案的位置位于我们这个答案矩阵的第几个位子，若让我们的答案为答案矩阵的第二项，我们将进行 \(n - 2\) 次变换，然后得到 \(\begin{bmatrix} f_{n - 1} & f_{n} & f_{n + 1} \end{bmatrix}\)。

code

const ll N = 3, P =  1e9 + 7;

ll n;

struct Matrix
{
    ll a[N][N];

    Matrix operator *(const  Matrix &x) const
    {
        Matrix res;
		memset(res.a, 0, sizeof res.a);
        fos(i, 0, N - 1)  
            fos(k, 0, N - 1)
                fos(j, 0, N - 1)
                    res.a[i][j] = (res.a[i][j] + (a[i][k] * x.a[k][j]) % P) % P;
	
        return res;
    }
} ans, base;

inline void init()
{
    base.a[1][0] = base.a[1][2] = base.a[2][1] = base.a[2][2] = 1;
    ans.a[0][0] = ans.a[0][1] = 1, ans.a[0][2] = 2;
}

inline void qmi(ll k)
{
    while(k)
    {
        if(k & 1) ans = ans * base;
        base = base * base;
        k >>= 1;
    }
}

int main()
{   
    read(n); init();
    if(n <= 2) {puts("1"); return 0;} 
    ll k = n - 2;
	qmi(k);
	
	cout << ans.a[0][2] << endl;
    return 0;
}

但是我们发现，我们构造出来的转移矩阵中，第一横行都为 \(0\)，所以我们就可以考虑，省去一维。

新的答案矩阵 \(\begin{bmatrix}f_1 & f_2\end{bmatrix}\)，新的转移矩阵为 \(\begin{bmatrix} 0 & 1 \\ 1 & 1 \end{bmatrix}\)。剩下的代码读者自己实现不难。

---

T1

广义斐波那契数列

这个可以作为一个很好的练习题，既不难，又可以检验之前学会了没有。

这个时候我们知道了广义斐波那契数列的递推式：

\[f_{i} \gets p \times f_{i - 1} + q \times f_{i - 2} \]

还是和上面的步骤一样，构造矩阵，然后转移。

不难得出我们的答案矩阵为 \(\begin{bmatrix} f_{1} f_{2}\end{bmatrix}\) ，但是由于我们的递推式的不同，我们的转移矩阵将变成 \(\begin{bmatrix} 0 & q \\ 1 & p \end{bmatrix}\)

之后即是相同的矩阵快速幂。

---

数学作业

这个题目用来做一个提升感觉并不为过。

这个题目我们将会用到一种技巧，因为我们在转移的时候不是仅仅只有 dp 数组的转移，还要考虑带个常数，所以说我们会采取一个扩维的方法，单独扩出来一维，放上 1，然后转移我们那个拼接的数字。

之后就是我们的转移方程了（不会有人不会推吧）：

\[f_{i} = f_{i-1} * 10 ^ {\lfloor {\log_{10} i} \rfloor + 1} + i \]

\[\begin{bmatrix} f_i \\ i \\ 1 \end{bmatrix} = \begin{bmatrix} f_{i - 1} \\ i - 1 \\ 1\end{bmatrix} \times \begin{bmatrix} 10 ^ {\lfloor\log_{10} i \rfloor} + 1 & 0 & 0 \\ 1 & 1 & 0 \\ 0 & 1 & 1\end{bmatrix} \]

---

Power of Matrix

可以考虑平移到我写的题解

我们这个题目可以用来学习到一个解决问题的方法。

对于我们要求的矩阵的幂次和，我们可以采取二次分治的思想。

不断递归然后求解。具体详细内容见我题解。

---

迷路

这个题目又会给我们引申出来一个知识点：拆点。

对于这个拆点这个概念我们可以对应分层图的思想来看。

这是对于要运用矩阵乘法但是边权不仅仅只为一的时候，我们可以在两个边权不为一的点建立边权个点。中间的虚点仅做转移运算用，而不存在任何实际意义，所以最后统计答案不可以加上虚点的答案。

对于这个题目：我们可以将每个点拆成九个，然后每一个点 \((i, j)\) 表示的是距离 \(i\) 节点为 \(j + 1\) 的一个虚点。为什么是 \(j + 1\) 是因为题目中：

不可以在某一个点逗留。

我们矩阵 \(f_{t, i, j} = k\) 表示的是 \(i\) 到 \(j\) 的距离为 \(t\) 的路径条数为 \(k\) 条。

然后经过 \(t\) 次变化就可以计算出来我们终点的值。

附上代码：

code

#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define rl register ll
#define foa(i, a, b) for(rl i=a; i < b; ++ i)
#define fos(i, a, b) for(rl i=a; i <= b; ++ i)
#define fop(i, a, b) for(rl i=a; i >= b; -- i)
#define ws putchar(' ')
#define wl putchar('\n')

template <class T> inline void waw(T x)
{
    if(x > 9) waw(x / 10);
    putchar(x % 10 ^ 48);
}

template <class T> inline void ww(T x)
{
    if(x < 0) x = ~x + 1, putchar('-');
    waw(x);
}

template <class T> inline void read(T &res)
{
    char ch = getchar(); bool f = 0; res = 0;
    for(; !isdigit(ch); ch = getchar()) f |= ch == '-';
    for(; isdigit(ch); ch = getchar()) res = (res << 1) + (res << 3) + (ch ^ 48);
    res = f ? ~res + 1 : res;
}

const ll N = 12, M = 128, P = 2009;

ll n, t, np;

struct Matrix
{
    ll a[M][M];

    Matrix() {memset(a, 0, sizeof a); }

    Matrix operator *(const Matrix &x) const
    {
        Matrix res;

        foa(i, 1, M) foa(k, 1, M) foa(j, 1, M)
            res.a[i][j] = (res.a[i][j] + (a[i][k] * x.a[k][j]) % P) % P;
        
        return res;
    }
} a;

inline ll id(ll x, ll y) { return x + y * n; }

inline Matrix qmi(ll k)
{
    Matrix res = a; -- k; 

    while(k)
    {
        if(k & 1) res = res * a;
        a = a * a;
        k >>= 1;
    }

    return res;
}

int main()
{
    read(n), read(t);

    fos(i, 1, n) 
    {
        fos(j, 1, 8) a.a[id(i, j)][id(i, j - 1)] = 1; // 沟通上下每一层

        fos(j, 1, n)
        {
            char ch = getchar(); ch ^= 48;
            if(ch) a.a[i][id(j, ch - 1)] = 1; // 建立我们的拆点，类似于分层图的构造。
        }
        getchar();
    }

    Matrix ans = qmi(t);

    waw(ans.a[1][n]);
    return 0;
}

---

CF222E

我们看到 \(n\) 的范围极其大，所以考虑倍增思想。

我们不难可以考虑到我们这个计数的式子是这样：

\[f_{i, j} \gets \sum _ {k = 1} ^ {52} f_{i - 1, k} \times [k \in \text{j前面的合法字母}] \]

我们得到这个式子之后不难想到后面这个维度判断我们可以构造一个矩阵 \(base\)，里面的 \(base_{i, j}\) 表示 \(i\) 是 \(j\) 这个字母前面的合法字母。

再构造一个 \(ans\) 矩阵，\(1, i\) 表示当前是以 \(i\) 结尾的方案数。我们这个当前表示我们填这个合法 \(DNA\) 到达的位置。这个显然是可以用快速幂来解决的。

code

快读快写部分略去。

const ll N = 55, P = 1e9 + 7;

ll n, m, k;

unordered_map<char, ll> mp;

struct Matrix
{
    ll a[N][N];

    Matrix () { memset(a, 0, sizeof a); }

    Matrix operator *(const Matrix &x) const
    {
        Matrix res;

        foa(i, 1, N) foa(k, 1, N) foa(j, 1, N)
            res.a[i][j] = (res.a[i][j] + (a[i][k] * x.a[k][j]) % P) % P;
        
        return res;
    }
} ans, base;

inline void qmi(ll k)
{
    while(k)
    {
        if(k & 1) ans = ans * base;
        base = base * base; k >>= 1;
    }
}

int main()
{
    read(n), read(m), read(k);

    fos(i, 0, 25) mp[(char)('a' + i)] = i + 1;
    fos(i, 0, 25) mp[(char)('A' + i)] = i + 27;

    fos(i, 1, m) ans.a[1][i] = 1;

    fos(i, 1, m) fos(j, 1, m) base.a[i][j] = 1;

    fos(i, 1, k)
    {
        char str[3]; cin >> str;
        ll a, b;
        a = mp[str[0]], b = mp[str[1]];
        base.a[a][b] = 0;
    }

    qmi(n - 1);

    ll res = 0;

    fos(i, 1, m) res = (res + ans.a[1][i]) % P;

    ww(res), wl;
    return 0;
}

---

CF576D

这个题目我们要引申出来一个 \(trick\) (?)

对于一个邻接矩阵 \(G\)，\(G ^ k\) 就是经过 \(k\) 后的状态。

运用这个知识我们就能想到，\(G\) 是我们的地图 \(k\) 是我们走过了多少条边，然后

---

P6190