浅谈 DDP 与 广义矩阵乘法

更好的阅读体验戳此进入

引入

例题 #1

题目 SP1043 GSS1 - Can you answer these queries I。

题面：对于序列 $a_n$$a_n$，$q$$q$ 次询问求 $[l_i,r_i]$$\left[ l_i, r_i \right]$ 的区间最大子段和。

$1\le n\le 5\times {10}^4,a_i\le 15007$$1 \le n \le 5 \times 10^4, a_i \le 15007$。$q$$q$ 在 int 范围内。

有显然的 DP，令 $dp(i)$$dp(i)$ 表示以 $a_i$$a_i$ 结尾的区间 $[1,i]$$\left[ 1, i \right]$ 的最大子段和。

一般的区间最大子段和转移为 $dp(i)=max(dp(i-1)+a_i,a_i)$$dp(i) = \max(dp(i - 1) + a_i, a_i)$，对于本题则为：

但是每次询问都重新维护一遍显然会寄，所以考虑一些人类智慧。

广义矩阵乘法

首先对于普通的矩阵乘法，我们随便举个例子：

我们考虑把这玩意省略的东西写全一点：

我们发现这里一共有两种符号，$+$$+$ 和 $\times$$\times$，考虑替换这两个符号，假设我们存在两个二元运算符 $\oplus$$\oplus$ 和 $\otimes$$\otimes$，替换一下就是：

这就是广义矩阵乘法了。

然后为了让这东西能有点用，所以现在我们还需要让广义矩阵乘法具有结合律，也就是：

或者写成：

展开一下有：

不难发现，如果 $\otimes$$\otimes$ 对 $\oplus$$\oplus$ 存在分配律，或者说满足 $(a\oplus b)\otimes c=(a\otimes c)\oplus (b\otimes c)$$(a \oplus b) \otimes c = (a \otimes c) \oplus (b \otimes c)$，那么把上面的两个式子拆开之后，最后就会发现实际上是相同的。

所以我们便可得出这个结论——如果 $\otimes$$\otimes$ 对 $\oplus$$\oplus$ 存在分配律，那么形成的广义矩阵乘法矩阵之间就满足结合律。

DDP

为了便于理解，我们先引入一个更简单的例子。

例题 #0

题面：求斐波那契数列。

十分显然的 $dp(i)=dp(i-1)+dp(i-2)$$dp(i) = dp(i - 1) + dp(i - 2)$。

我们考虑一些让复杂度更高奇妙的方法。

我们考虑把这个转移矩阵记作 $X$$X$，转移的起始和终点分别写成一个矩阵，为了凑齐矩阵，可以把终点塞一个 $0$$0$，或者 $dp(i-1)$$dp(i - 1)$ 这种无关的值，如替换为 $0$$0$，则有：

此时的 $\otimes$$\otimes$ 就是 $\times$$\times$，$\oplus$$\oplus$ 就是 $+$$+$，也就是说我们这个矩阵之间的乘法就是经典的矩阵乘法，则不难算出我们需要的矩阵 $X$$X$ 为：

也就是说：

于是我们不如继续尝试把 $\mathtt{\text{LHS}}$$\texttt{LHS}$ 的带有 $dp$$dp$ 的式子接着拆下去：

这东西虽然可以转移的，但是自己观察一下，这东西的形式和最开始的柿子好像形式很相似，于是我们考虑如果不去选择 $0$$0$ 而去选择 $dp(i-1)$$dp(i - 1)$ 可能会更方便一些（当然选择 $0$$0$ 也是可以转移的，只是会更复杂一些）。所以我们最开始的式子就变成了：

然后这样也很简单就能求出来 $X$$X$ 了：

然后我们继续往下转移：

此时我们便能发现一些奇妙的推导：

于是这个时候我们便可以矩阵快速幂（因此我们才需要使广义矩阵之间满足结合律，否则用不了矩阵快速幂或者线段树猫树之类的）把原来的线性求，变成现在的 $\log$$\log$ 求了（虽然常数大了点）。

当然这只是最最基础的 DDP。

例题 #0.5

题目：LG-P1115 最大子段和。

最大字段和模板题，因为比较水所以这里就不叙述的太详细了，令 $dp(i)$$dp(i)$ 表示以 $a_i$$a_i$ 结尾的区间 $[1,i]$$\left[ 1, i \right]$ 的最大子段和，有转移 $dp(i)=max(dp(i-1)+a_i,a_i)$$dp(i) = \max(dp(i - 1) + a_i, a_i)$，但是这样求完是还需要再扫一遍取 $max$$\max$，这样我们就不太好用矩阵转移，于是考虑再令 $ans(i)=max(dp(i),ans(i-1))$$ans(i) = \max(dp(i), ans(i - 1))$。于是我们不难推出如下转移：

同时注意这里的矩阵是广义矩阵乘法，$\otimes$$\otimes$ 为 $+$$+$，$\oplus$$\oplus$ 为 $max$$\max$，显然 $+$$+$ 对 $max$$\max$ 有分配律，则矩阵有结合律。

然后像之前一样拆下去就可以得到一大串矩阵相乘，成功大幅增加常数成功提高了转移的扩展性。

例题 #1

现在让我们回到文章伊始提到的例题，这里我们再回顾一下：

题目 SP1043 GSS1 - Can you answer these queries I。

题面：对于序列 $a_n$$a_n$，$q$$q$ 次询问求 $[l_i,r_i]$$\left[ l_i, r_i \right]$ 的区间最大子段和。

$1\le n\le 5\times {10}^4,a_i\le 15007$$1 \le n \le 5 \times 10^4, a_i \le 15007$。$q$$q$ 在 int 范围内。

有显然的 DP，令 $dp(i)$$dp(i)$ 表示以 $a_i$$a_i$ 结尾的区间 $[1,i]$$\left[ 1, i \right]$ 的最大子段和。

一般的区间最大子段和转移为 $dp(i)=max(dp(i-1)+a_i,a_i)$$dp(i) = \max(dp(i - 1) + a_i, a_i)$，对于本题则为：

和上一题一样，我们还是考虑令：

$i=l$$i = l$ 就是个特判可以不用考虑，考虑把 $l<i\le r$$l \lt i \le r$ 的式子变成广义矩乘：

然后继续往下拆：

一直拆到边界，就是：

因为我们要的答案为 $ans(r)$$ans(r)$，所以把 $i$$i$ 换成 $r$$r$。并且显然有 $ans(l)=dp(l)=a_l$$ans(l) = dp(l) = a_l$，则原式化为：

可以发现我们想要得到 $ans(r)$$ans(r)$，只需要从已知具体值的 $\left[\begin{array}{ccc}{a}_l& a_l& 0\end{array}\right]$$\begin{bmatrix} a_l & a_l & 0 \end{bmatrix}$，乘上一大堆矩阵。

为了便于书写，我们记：

则：

而在实际使用中我们还有很多小技巧与转化来让我们实现起来更简单：

如对于矩阵我们一般通过新开一个结构体然后重载 * 实现，这个时候我们可以考虑把初始矩阵补齐成 $3\times 3$$3 \times 3$ 的，这样最终答案也会变成 $3\times 3$$3 \times 3$，我们便只需要维护长宽均为 $3$$3$ 的定大小的矩阵即可，至于补齐的值是什么这个实际上是任意的，因为我们不需要这些补齐的值算出来的答案，所以最终计算出什么也是无关紧要的，计算后直接输出答案矩阵的 $(0,0)$$(0, 0)$ 位置，也就是 $ans(r)$$ans(r)$ 即可。

举个例子：

然后在重载矩阵的乘法的时候，我们应该注意，是否有可能会超出 int 或 long long 等，所以在我们取 $max$$\max$ 的时候可以考虑额外加一个 $-\mathrm{\infty }$$-\infty$，也就是当有两个数相加比 $-\mathrm{\infty }$$-\infty$ 更小的时候那么其中一定至少有一个数为 $-\mathrm{\infty }$$-\infty$，这个时候我们直接将答案变成 $-\mathrm{\infty }$$-\infty$ 即可。

然后如果严格按照上面的算法实现出来后，会发现一个问题，在某些特殊的时候答案会错误，经过 debug 我们发现当 $l=r$$l = r$ 的时候，本来应该只有最初的初始矩阵，但是按照我们的算法会多出来一个 $[l+1,r]$$\left[ l + 1, r \right]$ 的查询，但此时 $l+1>r$$l + 1 \gt r$，我们当然可以考虑特判这个情况，但是我们实际上还可以再次转换一下思路，让这个东西更泛用，更美观。

我们考虑能否将初始矩阵 $\left[\begin{array}{ccc}{a}_l& a_l& 0\end{array}\right]$$\begin{bmatrix}a_l & a_l & 0\end{bmatrix}$ 也写成 $A_l$$A_l$，然后在最前面再弄一个新的初始矩阵 $T$$T$，让它和 $A_l$$A_l$ 运算后会变成我们想要的 $\left[\begin{array}{ccc}{a}_l& a_l& 0\end{array}\right]$$\begin{bmatrix}a_l & a_l & 0\end{bmatrix}$，于是很容易写出：

不难算出：

所以我们的式子就最终变成了：

现在这个式子就很好维护了。

这个时候我们观察那一大串 $A$$A$，不难发现这玩意好像就是把 $[l,r]$$\left[ l, r \right]$ 的 $A$$A$ 乘到一起了，而且每次询问会有不同的 $l,r$$l, r$，这不就是区间查询吗，我们于是可以考虑用一些妙妙数据结构去维护这玩意。

此时我们可能会想到树状数组，但是仔细思考一下就会发现是不行的。众所周知树状数组在运算上有一个类似差分的过程，也就是说需要有逆元的存在，但是在我们定义的矩阵和矩阵之间的广义矩乘构成的群中是否存在逆元呢？显然我们都知道矩阵求逆的一个必要条件是行列式不为 $0$$0$，这东西是否不为 $0$$0$ 是无法确定的，也就是说不一定会存在逆元，所以不能用树状数组或者前缀积之类的东西区间查询。

一个比较显然的东西就是线段树，每个节点存的是一个矩阵，或者说对应区间的矩阵乘起来（看到这里大概就又能意识到让矩阵满足结合律的重要性了）。当然这里的乘仍然指的是 $\otimes$$\otimes$ 为 $+$$+$，$\oplus$$\oplus$ 为 $max$$\max$ 的广义矩乘。用线段树维护的话最后的复杂度大概是 $O({\xi }^3(n+q))$$O(\xi^3 (n + q))$，其中 $\xi$$\xi$ 为矩阵的行数或列数，这里 $\xi =3$$\xi = 3$。

对于本题来说，我们可以发现是多次查询没有修改，所以也可以考虑使用猫树，使查询复杂度变为 $O(1)$$O(1)$，这样总复杂度就变成了 $O({\xi }^3(n\log n+q))$$O(\xi^3(n \log n + q))$，具体介绍可以看咕咕日报 一种神奇的数据结构——猫树。

当然用平衡树，分块之类的也能做，这里不再过多赘述。

这里提供一个 DDP + 猫树 实现的 GSS1。（注意猫树节点数需要补齐到 $2^t$$2^t$，且要开好数组范围的大小，否则很容易 RE）

xxxxxxxxxx
117
1
#define _USE_MATH_DEFINES
2
#include <bits/extc++.h>
3

4
#define PI M_PI
5
#define E M_E
6
#define npt nullptr
7
#define SON i->to
8
#define OPNEW void* operator new(size_t)
9
#define ROPNEW(arr) void* Edge::operator new(size_t){static Edge* P = arr; return P++;}
10

11
using namespace std;
12
using namespace __gnu_pbds;
13

14
mt19937 rnd(random_device{}());
15
int rndd(int l, int r){return rnd() % (r - l + 1) + l;}
16
bool rnddd(int x){return rndd(1, 100) <= x;}
17

18
typedef unsigned int uint;
19
typedef unsigned long long unll;
20
typedef long long ll;
21
typedef long double ld;
22

23
#define MAXN (int)(5e4)
24
#define INF (0x3f3f3f3f)
25

26
template< typename T = int >
27
inline T read(void);
28

29
int N, M;
30
int lg[MAXN << 3];
31
int pos[MAXN << 1];
32
int a[MAXN];
33

34
struct Matrix3{
35
    int val[3][3];
36
    Matrix3(int v00, int v01, int v02, int v10, int v11, int v12, int v20, int v21, int v22):
37
        val{
38
            {v00, v01, v02},
39
            {v10, v11, v12},
40
            {v20, v21, v22}
41
        }{;}
42
    Matrix3(int val[][3]){for(int i = 0; i <= 2; ++i)for(int j = 0; j <= 2; ++j)this->val[i][j] = val[i][j];}
43
    Matrix3(void) = default;
44
    friend const Matrix3 operator * (const Matrix3 &x, const Matrix3 &y){
45
        int val[3][3];
46
        for(int i = 0; i <= 2; ++i)for(int j = 0; j <= 2; ++j)val[i][j] = -INF;
47
        for(int i = 0; i <= 2; ++i)for(int j = 0; j <= 2; ++j)for(int p = 0; p <= 2; ++p)
48
            val[i][j] = max({val[i][j], x.val[i][p] + y.val[p][j], -INF});
49
        return Matrix3(val);
50
    }
51
    void Print(void){
52
        for(int i = 0; i <= 2; ++i)for(int j = 0; j <= 2; ++j)
53
            printf("%d%c", val[i][j], j == 2 ? '\n' : ' ');
54
    }
55
}mt[MAXN];
56

57
class CatTree{
58
private:
59
    Matrix3 tr[30][MAXN << 1];
60
    #define MID ((gl + gr) >> 1)
61
    #define LS (p << 1)
62
    #define RS (LS | 1)
63
public:
64
    void Build(int p = 1, int dep = 1, int gl = 1, int gr = N){
65
        // printf("Building l = %d, r = %d, p = %d, dep = %d\n", gl, gr, p, dep);
66
        if(gl == gr)return pos[gl = gr] = p, void();
67
        tr[dep][MID] = mt[MID];
68
        tr[dep][MID + 1] = mt[MID + 1];
69
        //Tips: Matrix operation does not have the commutative law.
70
        for(int i = MID - 1; i >= gl; --i)tr[dep][i] = mt[i] * tr[dep][i + 1];
71
        for(int i = MID + 1 + 1; i <= gr; ++i)tr[dep][i] = tr[dep][i - 1] * mt[i];
72
        Build(LS, dep + 1, gl, MID);
73
        Build(RS, dep + 1, MID + 1, gr);
74
    }
75
    Matrix3 Query(int l, int r){
76
        if(l == r)return mt[l = r];
77
        int dep = lg[pos[l]] - lg[pos[l] ^ pos[r]];
78
        // tr[dep][l].Print(), tr[dep][r].Print();
79
        return tr[dep][l] * tr[dep][r];
80
    }
81
}ct;
82

83
int main(){
84
    N = read();
85
    int rN(1); while(rN < N)rN <<= 1;
86
    for(int i = 1; i <= N; ++i)
87
        a[i] = read(),
88
        mt[i] = Matrix3(0, -INF, -INF, a[i], a[i], -INF, a[i], a[i], 0);
89
    N = rN;
90
    lg[0] = lg[1] = 1;
91
    for(int i = 2; i <= (N << 2) + 10; ++i)lg[i] = lg[i >> 1] + 1;
92
    ct.Build();
93
    M = read();
94
    for(int m = 1; m <= M; ++m){
95
        int l = read(), r = read();
96
        auto ans = Matrix3(-INF, 0, 0, -INF, -INF, -INF, -INF, -INF, -INF) * ct.Query(l, r);
97
        printf("%d\n", ans.val[0][0]);
98
    }
99
    fprintf(stderr, "Time: %.6lf\n", (double)clock() / CLOCKS_PER_SEC);
100
    return 0;
101
}
102

103
template < typename T >
104
inline T read(void){
105
    T ret(0);
106
    short flag(1);
107
    char c = getchar();
108
    while(c != '-' && !isdigit(c))c = getchar();
109
    if(c == '-')flag = -1, c = getchar();
110
    while(isdigit(c)){
111
        ret *= 10;
112
        ret += int(c - '0');
113
        c = getchar();
114
    }
115
    ret *= flag;
116
    return ret;
117
}

例题 #2

题目：[ABC246Ex] 01? Queries。

题面：给定长度为 $N$$N$ 的仅包含 0，1，? 的字符串 $S$$S$，给定 $Q$$Q$ 组询问 $(x_1,c_1),(x_2,c_2),\cdots ,(x_q,c_q)$$(x_1, c_1), (x_2, c_2), \cdots, (x_q, c_q)$，每次将原字符串中 $x_i$$x_i$ 位置的字符改为 $c_i$$c_i$，然后输出 $S$$S$ 有多少种非空子串，? 需任意替换为 0 或 1。

$1\le N,Q\le {10}^5,1\le x_i\le N$$1 \le N, Q \le 10^5, 1 \le x_i \le N$。

我们先不考虑修改，思考对于这样一个字符串能有多少种子串，显然这个东西是个 DP。

令 $dp(i,0/1)$$dp(i, 0 / 1)$ 表示考虑前 $i$$i$ 位以 $0$$0$ 或 $1$$1$ 结尾的方案数，显然有如下转移：

若 $S_i=1$$S_i = 1$，有：

若 $S_i=0$$S_i = 0$，有：

若 $S_i=\mathtt{\text{?}}$$S_i = \texttt{?}$，有：

应该不难理解吧？如果状态里是当前这一位，那么可以接到上一个状态任意的结尾，接上这一位之后都会是一个符合要求的新串，或者丢弃以前的直接让这一位成为一个新串。反之就直接由上一次的转移而来，把这一位丢弃，而 $\mathtt{\text{?}}$$\texttt{?}$ 可以认为是任意一个，所以可以接在任意一个的后面。

这里有一个细节可以解释一下，在这一位和状态相同的时候我们会 $+1$$+1$，但是如果之前就有一个孤立的（这时很显然会有的），比如 $0$$0$，那么我们这一位又一个新的 $0$$0$ 难道不会重吗？这显然是不会的，因为上一个状态的 $0$$0$ 在这一次的讨论中已经变成 $00$$00$ 了，只有在讨论到上界的时候才会存在真正的孤立的 $0$$0$ 或 $1$$1$，于是显然可以保证不重不漏。

这么一大坨分类讨论显然很难维护，那么我们尝试把这些缩到一起：

关于式子中的方括号括起来的表达式，这个一般用来表示如果里面的表达式为 $\mathtt{\text{true}}$$\texttt{true}$，那么这个东西值为 $1$$1$，反之为 $0$$0$。

现在，我们就可以对这个简洁的式子搞事情了，因为我们后面要有很多次修改，这样的话显然可以考虑 DDP。

我们再次尝试设计出转移的矩阵：

不难算出：不难算出我们算不出来。

于是考虑再加一维

然后尝试推一下 $T$$T$，这东西个人认为实际上就是个蒙猜凑，总之最后弄一下就能出来这个（注意这里我们用的是普通矩乘规则，非广义）：

（不要问我为什么这题用不到广义矩乘还要在前面说一大堆广义矩乘的内容，我记错了，技多不压身）

然后类比之前的过程，令：

那么一直拆下去，一定有：

依然尝试找到一个这样的等式：

所以可以最终化为：

终于推完了，不难发现我们每次的修改就是 $A_i$$A_i$ 的值，然后查询整个区间。

这玩意有修改显然用不了猫树了，老老实实写线段树吧。。

复杂度大概是 $O({\xi }^3(n+q\log n))$$O(\xi^3(n + q \log n))$，其中 $\xi =3$$\xi = 3$。

xxxxxxxxxx
117
1
#define _USE_MATH_DEFINES
2
#include <bits/extc++.h>
3

4
#define PI M_PI
5
#define E M_E
6
#define npt nullptr
7
#define SON i->to
8
#define OPNEW void* operator new(size_t)
9
#define ROPNEW(arr) void* Edge::operator new(size_t){static Edge* P = arr; return P++;}
10

11
using namespace std;
12
using namespace __gnu_pbds;
13

14
mt19937 rnd(random_device{}());
15
int rndd(int l, int r){return rnd() % (r - l + 1) + l;}
16
bool rnddd(int x){return rndd(1, 100) <= x;}
17

18
typedef unsigned int uint;
19
typedef unsigned long long unll;
20
typedef long long ll;
21
typedef long double ld;
22

23
#define MAXN (int)(1e5 + 100)
24
#define MOD (int)(998244353)
25

26
template< typename T = int >
27
inline T read(void);
28

29
int N, Q;
30
int S[MAXN];
31

32
struct Matrix3{
33
    int val[3][3];
34
    Matrix3(int v00, int v01, int v02, int v10, int v11, int v12, int v20, int v21, int v22):
35
        val{
36
            {v00, v01, v02},
37
            {v10, v11, v12},
38
            {v20, v21, v22}
39
        }{;}
40
    Matrix3(int S):
41
        val{
42
            {1, S != 0, 0},
43
            {S != 1, 1, 0},
44
            {S != 1, S != 0, 1}
45
        }{;}
46
    Matrix3(int val[][3]){for(int i = 0; i <= 2; ++i)for(int j = 0; j <= 2; ++j)this->val[i][j] = val[i][j];}
47
    Matrix3(void) = default;
48
    friend const Matrix3 operator * (const Matrix3 &x, const Matrix3 &y){
49
        int val[3][3]; memset(val, 0, sizeof val);
50
        for(int i = 0; i <= 2; ++i)for(int j = 0; j <= 2; ++j)for(int p = 0; p <= 2; ++p)
51
            val[i][j] = ((ll)val[i][j] + (ll)x.val[i][p] * y.val[p][j] % MOD) % MOD;
52
        return Matrix3(val);
53
    }
54
    void Print(void){
55
        for(int i = 0; i <= 2; ++i)for(int j = 0; j <= 2; ++j)
56
            printf("%d%c", val[i][j], j == 2 ? '\n' : ' ');
57
    }
58
}mt[MAXN];
59

60
class SegTree{
61
private:
62
    Matrix3 tr[MAXN << 2];
63
    #define LS (p << 1)
64
    #define RS (LS | 1)
65
    #define MID ((gl + gr) >> 1)
66
public:
67
    void Pushup(int p){tr[p] = tr[LS] * tr[RS];}
68
    void Build(int p = 1, int gl = 1, int gr = N){
69
        if(gl == gr)return tr[p] = mt[gl = gr], void();
70
        Build(LS, gl, MID);
71
        Build(RS, MID + 1, gr);
72
        Pushup(p);
73
    }
74
    void Modify(int idx, Matrix3 v, int p = 1, int gl = 1, int gr = N){
75
        if(gl == gr)return tr[p] = v, void();
76
        if(idx <= MID)Modify(idx, v, LS, gl, MID);
77
        else Modify(idx, v, RS, MID + 1, gr);
78
        Pushup(p);
79
    }
80
    Matrix3 Query(void){return tr[1];}
81
}st;
82

83
int main(){
84
    N = read(), Q = read();
85
    string s; cin >> s;
86
    for(int i = 1; i <= (int)s.size(); ++i)
87
        S[i] = s.at(i - 1) == '?' ? -1 : int(s.at(i - 1) - '0'),
88
        mt[i] = Matrix3(S[i]);
89
    st.Build();
90
    Matrix3 origin(0, 0, 1, 0, 0, 0, 0, 0, 0);
91
    while(Q--){
92
        int p = read();
93
        char c = getchar(); while(c != '0' && c != '1' && c != '?')c = getchar();
94
        int flag = c == '?' ? -1 : int(c - '0');
95
        st.Modify(p, Matrix3(flag));
96
        auto ans = origin * st.Query();
97
        printf("%d\n", (int)((ll)(ans.val[0][0] + ans.val[0][1]) % MOD));
98
    }
99
    fprintf(stderr, "Time: %.6lf\n", (double)clock() / CLOCKS_PER_SEC);
100
    return 0;
101
}
102

103
template < typename T >
104
inline T read(void){
105
    T ret(0);
106
    short flag(1);
107
    char c = getchar();
108
    while(c != '-' && !isdigit(c))c = getchar();
109
    if(c == '-')flag = -1, c = getchar();
110
    while(isdigit(c)){
111
        ret *= 10;
112
        ret += int(c - '0');
113
        c = getchar();
114
    }
115
    ret *= flag;
116
    return ret;
117
}

例题 #3

正在写~

UPD

update-2022_10_22 初稿