MATRIX,多么美妙的一个单词！

[TOC]

定义

矩阵是一个按照长方形排列的元素集合。简单地说，矩阵可以理解为一个二维数组，其中每一个位置都存放了一个元素。

例如，以下是一个大小为 n\times m 的 矩阵。

A= \left[ \begin{matrix} a_{11} & a_{12} & \cdots & a_{1m} \\ a_{21} & a_{22} & \cdots & a_{2m} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nm} \\ \end{matrix} \right]

运算

加减法

设矩阵 A_{n\times m},B_{n\times m},C_{n\times m}。

没啥好说的，对应位置相减就完事了。

C = A\pm B,C_{i,j} = A_{i,j} \pm B_{i,j}

乘法

设矩阵 A_{n\times m},B_{m\times p},C_{n\times p} = A \times B。

矩阵相乘相对比较难理解，我们先从矩阵乘向量开始理解。

矩阵乘向量

所谓向量，本质上就是一个 m = 1 的矩阵。直接用一个例子来展示矩阵乘向量。

A= \left[ \begin{matrix} 1 &3\\ 4 &0\\ 2 &1\\ \end{matrix} \right] B = \left[ \begin{matrix} 1\\ 5\\ \end{matrix} \right] \\ A\times B = \left[ \begin{matrix} 1\times 1 + 5 \times 3 \\ 4\times 1 + 0 \times 5 \\ 2\times 2 + 1 \times 5 \\ \end{matrix} \right] = \left[ \begin{matrix} 16 \\ 4 \\ 10 \\ \end{matrix} \right]

注意，这里一个 n\times m 的矩阵乘上一个长度为 m 的向量，得到的结果是一个长度为 n 的向量。

这个过程也说明了向量的长度必须与矩阵的列数相等，答案向量长度与矩阵的行数相等。

如果不理解的话这张图或许可以帮上你(From Zhihu:爱编程的高老师)：

注意这里时间复杂度实际上是 O(n^2) 的，[某些题目]()中可能需要利用到这个性质。

矩阵乘矩阵

矩阵相乘只有在第一个矩阵的列数和第二个矩阵的行数相同时才有意义。

设 A 为 P \times M 的矩阵，B 为 M \times Q 的矩阵，设矩阵 C 为矩阵 A 与 B 的乘积，

其中矩阵 C 中的第 i 行第 j 列元素可以表示为：

C_{i,j} = \sum_{k=1}^MA_{i,k}B_{k,j}

如果没看懂上面的式子，没关系。通俗的讲，在矩阵乘法中，结果 C 矩阵的第 i 行第 j 列的数，就是由矩阵 A 第 i 行 M 个数与矩阵 B 第 j 列 M 个数分别相乘再相加得到的。

例如：

\begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{bmatrix} * \begin{bmatrix} b_{11} & b_{12} \\ b_{21} & b_{22} \end{bmatrix} = \begin{bmatrix} a_{11}b_{11}+a_{12}b_{21} & a_{11}b_{12}+a_{12}b_{22} \\ a_{21}b_{11}+a_{22}b_{21} & a_{21}b_{12}+a_{22}b_{22} \end{bmatrix}

矩阵乘法满足结合律，不满足一般的交换律。利用结合律，矩阵乘法可以利用快速幂的思想来优化。

这启示我们实际上我们可以定义广义矩阵乘法，即不用加法和乘法。那怎么保证满足结合律呢？

需要加法满足交换律结合律，乘法满足结合律，加法和乘法满足分配律，即 a(b+c)=ab+ac，例如 \min 的话就是 \min(a+b, a+c)=a+\min(b, c) (此处感谢 UOJ 群的 nealchen 大佬)

在接下来会谈到应用。

代码：

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

    friend inline matrix operator * (const matrix &A,const matrix &B) {
        matrix c;
        for (int i = 1;i <= n;++i) {
            for (int j = 1;j <= n;++j) {
                c.a[i][j] = 0;
                for (int k = 1;k <= n;++k) {
                    c.a[i][j] = (c.a[i][j] + (A.a[i][k] * B.a[k][j]) % mod) % mod;
                }
            }
        }
        return c;
    }
}A,E,F;

matrix matffpow(matrix a,ll b) {
    matrix ans = E;
    for (;b;b >>= 1) {
        if (b & 1) ans = ans * a;
        a = a * a;
    }
    return ans;
}

矩阵加速递推

一维加速递推

例：Fibonacci 数列

斐波那契数列是满足如下性质的一个数列：

F_n = \left\{\begin{aligned} 1 \space (n \le 2) \\ F_{n-1}+F_{n-2} \space (n\ge 3) \end{aligned}\right.

请你求出 F_n \bmod 10^9 + 7 的值。

n \leq 2 ^{63}

我们如果直接递推的话时间复杂度为 O(n)，无法通过本题。

考虑把递推式里的东西丢进一个矩阵里优化递推过程，我们设

F_i = \begin{bmatrix} f_i \\ f_{i-1} \end{bmatrix}

那么我们怎么从 F_{i-1} 到 F_{i} 呢？考虑构造一个单位矩阵，因为我们有 f_i = f_{i-1} + f_{i-2}，故

F_i = \begin{bmatrix} f_{i-1} + f_{i-2} \\ f_{i-1} \end{bmatrix} = \begin{bmatrix} 1 \times f_{i-1} + 1\times f_{i-2} \\ 1 \times f_{i-1} + 0 \times f_{i-2} \end{bmatrix}

反推出单位矩阵

A = \begin{bmatrix} 1 &1 \\ 1 &0 \end{bmatrix}

不难验证成立。

所以我们就有 F_i = A \times F_{i-1} = A^{i-1} \times F_1

矩阵快速幂处理即可。

这类题目一般有一种比较特征的东西，就是递推的各个数都是 f_i 的形式，我将其称为一维递推。对于这种递推方式，一般的方法是：

1. 观察递推式中有几项是相关的，比如 Fibonacci 数列中就是 f_i 与 f_{i-1},f_{i-2} 相关。
2. 将 f_i 和相关项丢进矩阵中，按照上面提到的方法展开。
3. 如果这里构造了是一列的矩阵，按项数从大到小直接将系数提取出来就可以构造出单位矩阵了。

注意如果有常数项，也可以一并将常数项丢进去处理，如下题（From OI-wiki）

f_{1} = f_{2} = 0\\ f_{n} = 7f_{n-1}+6f_{n-2}+5n+4\times 3^n

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

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

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

\begin{bmatrix}f_n\\ f_{n-1}\\ n\\ 3^n \\ 1\end{bmatrix}

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

\begin{bmatrix} f_{n+1}\\ f_{n}\\ n+1 \\ 3^{n+1} \\ 1 \end{bmatrix}

将矩阵展开有

\begin{bmatrix} 7\times f_{n} + 6 \times f_{n-1} + 5 \times n + 4 \times 3 \times 3^n + 5\\ 1\times f_{n} + 0 \times f_{n-1} + 0 \times n + 0 \times 3^n + 0\\ 0\times f_{n} + 0 \times f_{n-1} + 1 \times n + 0 \times 3^n + 1 \\ 0\times f_{n} + 0 \times f_{n-1} + 0 \times n + 3 \times 3^n + 0 \\ 0\times f_{n} + 0 \times f_{n-1} + 0 \times n + 0 \times 3^n + 1 \end{bmatrix}

转移矩阵即为

\begin{bmatrix} 7 & 6 & 5 & 12 & 5\\ 1 & 0 & 0 & 0 & 0\\ 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 3 & 0\\ 0 & 0 & 0 & 0 & 1 \end{bmatrix}

请尝试用这种方法解决习题： P1939 【模板】矩阵加速（数列）

二维加速递推

这类递推一般就结合 DP 等算法来加速递推了，详细的说，如果存在一系列形如：

f_{i,0} = f_{i-1,i} + \dots \\ f_{i,1} = f_{i-1,j} + \dots \\ \dots \\ f_{i,k} = f_{i-1,p} + \dots \\

的 常系数齐次线性递推式 都可以表示成矩阵的形式。

也就是说，对于序列的每一项 a_i 我们都可以找到一个矩阵 A_i 使得

\begin{bmatrix} f_{i,0}\\ f_{i,1}\\ \dots\\ f_{i,j} \end{bmatrix} = A_i \times \begin{bmatrix} f_{i-1,0}\\ f_{i-1,1}\\ \dots\\ f_{i-1,j} \end{bmatrix}

例：CF750E New Year and Old Subsequence

定义一个数字串满足性质nice当且仅当：该串包含子序列2017，且不包含子序列2016。

定义一个数字串的函数ugliness为：该串至少删去几个字符，可以使得剩余串满足性质nice；如果该串没有满足性质nice的子序列，则该串的ugliness是-1。

给定一个长度为n的字符串t，和q次询问，每次询问用(l,r)表示。对于每次询问，回答ugliness(t[l,r])。

4 \leq n \leq 200000,1 \leq q \leq 200000

首先，如果这个 l,r 是给定的，我们可以写出如下的 DP 方程：

dp_{i,0}=\begin{cases}dp_{i-1,0} &(s[i]\neq'2')\\dp_{i-1,0}+1 &(s[i]='2')\end{cases}

dp_{i,1}=\begin{cases}dp_{i-1,0} &(s[i]\neq'0'\land s[i]\neq'2')\\\min(dp_{i-1,1},dp_{i-1,0}) &(s[i]='2')\\dp_{i-1,0}+1 &(s[i]='0')\end{cases}

dp_{i,2}=\begin{cases}dp_{i-1,1} &(s[i]\neq'1'\land s[i]\neq'0')\\\min(dp_{i-1,2},dp_{i-1,0}) &(s[i]='0')\\dp_{i-1,1}+1 &(s[i]='1')\end{cases}

dp_{i,3}=\begin{cases}dp_{i-1,3} &(s[i]\neq'7'\land s[i]\neq'1'\land s[i]\neq'6')\\\min(dp_{i-1,3},dp_{i-1,2}) &(s[i]='1')\\dp_{i-1,2}+1 &(s[i]='6'\lor s[i]='7')\end{cases}

dp_{i,4}=\begin{cases}dp_{i-1,4} &(s[i]\neq'7'\land s[i]\neq'6')\\\min(dp_{i-1,4},dp_{i-1,3}) &(s[i]='7')\\dp_{i-1,4}+1 &(s[i]='6')\end{cases}

这个方程转移的含义如下:

1. 若 a_i=2 。
   1. 不删去 a_i 。此时 a_i 的作用在于将一个暂未匹配的状态( \emptyset )变为一个匹配 1 位的状态(2)。在这种决策下有 dp_{i,1}=dp_{i-1,0} 。
   2. 删去 a_i 。可以发现删除 a_i 起到了保留状态的作用，可以不改变一个匹配 0 位的状态。在这种决策下有 dp_{i,0}=dp_{i-1,0}+1 。
2. 若a_i=0，与①同理:
   1. dp_{i,2} 的决策点包含 dp_{i-1,1} ；
   2. dp_{i,1} 的决策点包含 dp_{i-1,1}+1 。
3. 若 a_i=1/7 ，与 1,2 同理，请自行推导。
4. 若 a_i=6 ：
   1. 不删去 a_i 。这种转移合法当且仅当当前匹配到的位数不超过 2 ，否则转移不合法。
      为什么这里有一个状态转移的合法限制呢？假设匹配了3位，那么此时就含有子序列 201 ，添上一个6后就存在了子序列 2016 而这时题目严格禁止出现的。匹配了 4 位的情况同理。
      所以， dp_{i,0},dp_{i,1},dp_{i,2} 可以直接转移过来，但 dp_{i,3},dp_{i,4} 必须分别从 dp_{i-1,3},dp_{i-1,4} 加 1 而转移过来。
5. 若 a_i 为其他的数，那么不会有任何影响，直接转移过来即可。

注意到这里，我们需要定义一个广义的矩阵乘法，也就是把 \prod 换成 \sum , 把\sum 换成 \min,即 c_{i,j} = \sum_{1 \leq k \leq p} \min\{ a_{i,k} + b_{k,j} \}

通过上面的方法，我们可以证明这个是满足结合律的，所以我们就可以直接使用矩阵乘法处理。

由转移方程我们可以得出当 s_i = 2/0/1/6/7 的不同时候的转移矩阵。如果直接暴力乘，时间复杂度为 O(qn) ，因为矩阵乘法所以会带上一个 125 的常数，无法通过。

考虑矩阵乘法是没有除法的，所以我们可以想到使用线段树优化这个过程，每一次查询等价于查询区间矩阵乘积，时间复杂度 O(n \log n)，还是带一个 125 的常数。

Code:

int n,q,a[N];

struct matrix {
    int a[5][5];

    void init() {
        for (int i = 0;i < 5;++i) {
            for (int j = 0;j < 5;++j) {
                a[i][j] = 0x3f3f3f3f;
            }
        }
    }

    friend inline matrix operator * (const matrix &A,const matrix &B) {
        matrix C;C.init();
        for (int i = 0;i < 5;++i) {
            for (int j = 0;j < 5;++j) {
                for (int k = 0;k < 5;++k) {
                    C.a[i][j] = min(C.a[i][j],A.a[i][k] + B.a[k][j]);
                }
            }
        }
        return C;
    }
};

struct node {
    int l,r;matrix val;
}tree[N<<2];

void build(int p,int l,int r) {
    tree[p].l = l,tree[p].r = r;
    if (l == r) {
        tree[p].val.init();
        for (int i = 0;i < 5;++i) tree[p].val.a[i][i] = 0;
        if (a[l] == 2) tree[p].val.a[0][0] = 1,tree[p].val.a[0][1] = 0;
        if (a[l] == 0) tree[p].val.a[1][1] = 1,tree[p].val.a[1][2] = 0;
        if (a[l] == 1) tree[p].val.a[2][2] = 1,tree[p].val.a[2][3] = 0;
        if (a[l] == 6) tree[p].val.a[3][3] = 1,tree[p].val.a[4][4] = 1;
        if (a[l] == 7) tree[p].val.a[3][3] = 1,tree[p].val.a[3][4] = 0;
        return ;
    }
    int mid = (l + r) >> 1;
    build(p<<1,l,mid);build(p<<1|1,mid+1,r);
    tree[p].val = tree[p<<1].val * tree[p<<1|1].val;
}

matrix query(int p,int x,int y) {
    if (tree[p].l >= x && tree[p].r <= y) {
        return tree[p].val;
    }
    int mid = (tree[p].l + tree[p].r) >> 1;
    if (y <= mid) return query(p<<1,x,y);
    else if (x > mid) return query(p<<1|1,x,y);
    else return query(p<<1,x,mid) * query(p<<1|1,mid+1,y);
}

signed main() {
    read(n,q);
    for (int i = 1;i <= n;++i) {a[i] = getchar() - '0';}
    build(1,1,n);
    while (q--) {
        int x,y;read(x,y);
        int ans = query(1,x,y).a[0][4];
        if (ans >= 0x3f3f3f3f) puts("-1");
        else printf("%d\n",ans);
    }
}

习题： [NOI2020] 美食家