子序列动态规划
flywen
·
·
算法·理论
Part.0 前言
这是一篇写给自己的动态规划入门学习笔记，但我也希望这能帮助到其他因没有学会动态规划而苦恼的选手。
鉴于我从学 oi 开始就不擅长喜欢动态规划，如今都以 NOIP 为目标，却连最基础的橙色动态规划题都没有思路，所以我决定从零开始学习动态规划。
本篇文章是通过子序列这个专题来练习动态规划的，并不是要研究子序列。子序列是线性动态规划中最基础和常见的问题，所以本篇文章介绍一些简单的子序列问题。注意，这些问题可能存在字符串解法，但本文不进行讨论，像福建三问这种东西就不在考虑范围内。
---
### 记号约定
通常不区分时间复杂度渐进上界和渐进紧确界。时间复杂度中为显示代码流程有时不省略低阶项。
通常不严格要求各类求和等记号格式。
通常所有记号都是整数，如果没有取整符号，默认取下整。
---
### 代码
本人代码风格一坨，就看个乐吧。
尽量保持和正文一样的变量名。
不重要的宏定义和快读都删了。我大多数不会给完整的主函数，会将部分不重要的输入和`int main()`删掉。
---
### 相关链接
* [题单链接](https://www.luogu.com.cn/training/1020271)
# Part.1 常见问题
*子序列是指选择原序列的一些元素，不要求相邻或连续，按照原序列中的顺序重新排列组成的序列。*
*比较简单的子序列问题通常是通过维护具有某特征的子序列信息（如最后一个元素），计算一个新元素能否接在后面解决的。*
## 最长上升子序列
### [B3637 最长上升子序列](https://www.luogu.com.cn/problem/B3637)
给定长为 $n$ 的序列 $a_i$，求其最长**严格**上升子序列。
$n\leq5\times10^3$。
---
定义 $f_i$ 为**以 $i$ 结尾**的最长上升子序列长度，所以答案为 $\max\limits_{i=1}^nf_i$。
对于所有 $j<i$ 并且 $a_j<a_i$，都可以将 $a_i$ 接在 $a_j$ 后面，故有：
$$f_i=\max\limits_{1\leq j<i,a_j<a_i}f_j + 1$$
时间复杂度 $O(n^2)$。
```
for(int i = 1;i <= n;++i){
for(int j = 1;j < i;++j)if(a[j] < a[i])f[i] = max(f[i],f[j]);
ans = max(ans,++f[i]);
}
```
---
### 最长上升子序列优化
由上题可知如何在 $O(n^2)$ 时间复杂度内解决最长上升子序列，我们来看如何优化至 $O(n\log n)$。
以非严格的为例，如果我们盯着转移看，$f_i=\max\limits_{1\leq j\leq i,a_j\leq a_i}{f_j}$，是对于之前处理的中所有 $a_j\leq a_i$ 的 $f_j$ 取 $\max$，显然可以离散化后使用树形数据结构，如树状数组快速查询，时间复杂度为 $O(n\log n)$。
```cpp
void modify(int x,int v){ for(;x <= tot;x += x&(-x))c[x] = max(c[x],v); }
int query(int x){ int res = 0;for(;x;x -= x&(-x))res = max(res,c[x]);return res; }
for(int i = 1;i <= n;++i){
f[i] = query(rk[a[i]]) + 1; //如果是严格的则query(rk[a[i]] - 1)
modify(rk[a[i]],f[i]);
ans = max(ans,f[i]);
}
```
另一方面，我们发现我们定义的状态 $f_i$ 是和原数组 $a_i$ 相关的，而不是和最长上升子序列相关。
设当前最长的上升子序列长度为 $len$，我们考虑维护 $t_i$ 表示长度为 $i+1$ 的最长上升子序列的**末尾元素最小值**。
当新的元素 $a_i$ 加入时，如果 $a_i$ 大于等于 $t_{len}$，那么 $t_{len + 1}\leftarrow a_i$，然后 $len\leftarrow len+1$；否则在 $t$ 中**二分搜索**大于等于 $a_i$ 的最小元素，将其替换为 $a_i$。
根据上述操作，易知替换操作没有影响上升子序列的长度，只是将其末尾替换，因为它越小越有可能使最长上升子序列长度增加；$t$ 是**单调递增**的，但**不是**最长上升子序列。
最后的 $len$ 就是答案，时间复杂度 $O(n\log n)$。
```cpp
for(int i = 1;i <= n;++i){
if(t[len] <= b[i])t[++len] = b[i];
else{
int x = lowerbound(t + 1,t + 1 + len,a[i]) - t;
t[x] = b[i];
}
}
```
---
### [P1020 [NOIP 1999 提高组] 导弹拦截](https://www.luogu.com.cn/problem/P1020)
一套导弹拦截系统的第一发炮弹能够到达任意的高度，以后每一发炮弹都不能高于前一发。
对于一个长为 $n$ 的需拦截导弹高度序列，求：
**Task1**：一套系统最多能拦截多少导弹。
**Task2**：拦截所有导弹最少要配备多少套导弹拦截系统。
$n\leq10^5$。
---
##### Task1
显然是最长下降子序列。
##### Task2
贪心。
定义 $f_i$ 为当前第 $i$ 套系统上一次拦截的导弹高度，即 $i$ 下一次拦截的高度不能超过 $f_i$。
对新的导弹高度 $h$，二分第一个比 $h$ 大的 $f_i$，将其设为 $h$；否则新开一套系统。
和最长上升子序列类似，$f$ 一定是单调递增的。
时间复杂的 $O(n\log n)$。
## 最长公共子序列
### 最长公共子序列朴素解法
给定两个长分别为 $n,m$ 的序列 $a_i,b_j$，求其最长公共子序列。
$n,m\leq10^3$。
---
定义 $f_{i,j}$ 表示 $a$ 前 $i$ 位，$b$ 前 $j$ 位构成的最大公共子序列长度。
当 $a_i=b_j$ 时，用 $f_{i-1,j-1}+1$ 去更新 $f_{i,j}$。
此外继承之前的状态：$f_{i,j}=\max\{f_{i-1,j},f_{i,j-1}\}$。
时间复杂度 $O(nm)$。
```cpp
for(int i = 1;i <= n;++i){
for(int j = 1;j <= m;++j){
f[i][j] = max(f[i - 1][j],f[i][j - 1]);
if(a[i] == b[j])f[i][j] = max(f[i][j],f[i - 1][j - 1] + 1);
}
}
```
---
### [P1439 两个排列的最长公共子序列](https://www.luogu.com.cn/problem/P1439)
给定两个长为 $n$ 的 $1$ 到 $n$ 的排列 $a_i,b_j$，求其最长公共子序列。
$n\leq10^5$。
---
看似很难解决，但是将一个最长公共子序列在两排列中的对应元素连接起来，会是一排不交叉的线。这个理解是从 [ABC#439 E Kite](https://www.luogu.com.cn/problem/AT_abc439_e)/[P2782 友好城市](https://www.luogu.com.cn/problem/P2782) 来的。
```
3 2 1 4 5
/ | |
/ | |
/ | |
1 2 3 4 5
```
设 $v_{a_i}=i$，即每个数在 $a$ 中的下标，那么题目在求 $v_{b_i}$ 的最长上升子序列。
这是一个非常巧妙的转化，一定要仔细理解。
时间复杂度 $O(n\log n)$。
# Part.2 子序列建模
*本板块主要讲解子序列动态规划题目，难度主要提高及以上。*
## 计数/容斥
### [P1819 公共子序列](https://www.luogu.com.cn/problem/P1819)
求 $3$ 个长为 $n$ 的小写字母序列不同的公共子序列数量，不包括空序列，对 $10^8$ 取模。
$n\leq150$。
---
定义 $f_{i,j,k}$ 表示三个序列下标分别是前 $i,j,k$ 时的方案数。
当 $a_i,b_j,c_k$ 不相等时，根据**容斥原理**，我们可以得到：
$$f_{i,j,k}=f_{i-1,j,k}+f_{i,j-1,k}+f_{i,j,k-1}-f_{i-1,j-1,k}-f_{i,j-1,k-1}-f_{i-1,j,k-1}+f_{i-1,j-1,k-1}$$
但当 $a_i=b_j=c_k$ 时，我们要考虑将它们接到后面，那么接是一种，不接是一种，故总共情况是 $2\times f_{i-1,j-1,k-1}$。然而还有重复的，因为以前也有可能出现三个相同的，这是不接和这种情况是没有区别的，故：
$$f_{i,j,k}=2\times f_{i-1,j-1,k-1}-f_{lsta_i-1,lstb_j-1,lstc_k-1}$$
其中 $lsta,lstb,lstc$ 分别表示当前字符在该串中上一次出现的位置。
最后记得去除空串。
时间复杂度 $O(n^3)$。
```cpp
f[0][0][0] = 1;
for(int i = 0;i <= n;++i){
for(int j = 0;j <= n;++j){
for(int k = 0;k <= n;++k){
if(i && j && j && a[i] == b[j] && b[j] == c[k]){
f[i][j][k] = f[i - 1][j - 1][k - 1] * 2;
if(lsta[i] && lstb[j] && lstc[k])f[i][j][k] -= f[lsta[i] - 1][lstb[j] - 1][lstc[k] - 1];
}
else{
if(i)f[i][j][k] += f[i - 1][j][k];
if(j)f[i][j][k] += f[i][j - 1][k];
if(k)f[i][j][k] += f[i][j][k - 1];
if(i && j)f[i][j][k] -= f[i - 1][j - 1][k];
if(k && j)f[i][j][k] -= f[i][j - 1][k - 1];
if(i && k)f[i][j][k] -= f[i - 1][j][k - 1];
if(i && j && k)f[i][j][k] += f[i - 1][j - 1][k - 1];
}
}
}
}
```
---
### [P2516 [HAOI2010] 最长公共子序列](https://www.luogu.com.cn/problem/P2516)
给定两个长分别为 $n,m$ 的序列 $a_i,b_j$，求其最长公共子序列及其个数，对 $10^8$ 取模。
$n,m\leq5\times10^3$。
---
考虑我们的转移：
$$f_{i,j}=\max\{f_{i-1,j},f_{i,j-1},f_{i-1,j-1}+1|a_i=b_j\}$$
那么定义 $g_{i,j}$ 为方案数，当转移过来的是结果最大值那么就加上转移过来的方案数。
但是这样还不对，因为当 $f_{i-1,j-1}=f_{i,j}$ 时，说明 $f_{i-1,j}$ 和 $f_{i,j-1}$ 都从 $f_{i-1,j-1}$ 转移，所以 $g_{i-1,j}$ 和 $g_{i,j-1}$ 都对 $g_{i,j}$ 产生贡献，根据容斥原理，$g_{i,j}$ 还要减去 $g_{i-1,j-1}$。
边界 $g_{0,j}=g_{i,0}=1$。
数组要滚动一下。
时间复杂度 $O(n^2)$。
```
f[i&1][j] = max(f[(i^1)&1][j],f[i&1][j - 1]);
g[i&1][j] = 0;
if(a[i] == b[j]){
f[i&1][j] = max(f[i&1][j],f[(i^1)&1][j - 1] + 1);
if(f[i&1][j] == f[(i^1)&1][j - 1] + 1)g[i&1][j] += g[(i^1)&1][j - 1];
}
if(f[(i^1)&1][j - 1] == f[i&1][j])g[i&1][j] -= g[(i^1)&1][j - 1];
if(f[(i^1)&1][j] == f[i&1][j])g[i&1][j] += g[(i^1)&1][j];
if(f[i&1][j - 1] == f[i&1][j])g[i&1][j] += g[i&1][j - 1];
g[i&1][j] = (g[i&1][j] + mod) % mod;
```
---
### [P3970 [TJOI2014] 上升子序列](https://www.luogu.com.cn/problem/P3970)
给定长为 $n$ 的序列 $a_i$，求有多少长度至少为 $2$ 的严格上升子序列，对 $10^9+7$ 取模。
$n\leq10^5$。
---
先离散化，定义 $f_i$ 为以**数** $i$ 结尾的子序列个数，那么 $f_i=\sum\limits_{j=1}^if_j$，可以用树状数组计算。
根据之前 [P1819 公共子序列](https://www.luogu.com.cn/problem/P1819) 的经验，我们知道肯定要去重的，重复的当然是上一次 $i$ 出现时的情况，故减去即可。
不要偷懒去用 $lst$ 判断第一次出现，否则在取模情况下 $lst$ 会为 $0$，$90$ 分。
时间复杂度 $O(n\log n)$。
```cpp
for(int i = 1;i <= n;++i){
long long sum = query(a[i] - 1);
long long cnt = (sum - lst[a[i]] + mod) % mod;
ans = (ans + cnt) % mod;
if(!vst[a[i]]){ //第一次出现{a[i]}算一个,但不算到答案中
vst[a[i]] = 1;
cnt = (cnt + 1) % mod;
}
modify(a[i],cnt);
lst[a[i]] = sum;
}
```
## 差分
### [P2215 [HAOI2007] 上升序列](https://www.luogu.com.cn/problem/P2215)
给定长为 $n$ 的序列 $a_i$，$q$ 组问询，求是否存在长为 $L$ 的严格上升子序列，如果存在输出下标字典序最小的。保证数据随机。
$n,q\leq10^4$。
---
反过来，定义 $f_i$ 为以 $i$ 开头的最长上升子序列长度。
如果最大的 $f_i<L$ 则无解；否则每次问询从 $1$ 扫到 $n$，如果 $f_i\geq L$ 并且 $a_i$ 比上一次输出 $a_{lst}$ 的要大，就输出 $a_i$ 并将 $L\leftarrow L-1,lst\leftarrow i$。原理是转移一定是从后面的 $+1$ 来的。
时间复杂度 $O(nq)$。
```cpp
int L = in();
if(L > mx)printf("Impossible\n");
else{
int lst = 0;
for(int i = 1;i <= n;++i){
if(f[i] >= L && a[i] > a[lst]){
L--;
lst = i;
printf("%d ",a[i]);
if(!L)break;
}
}
printf("\n");
}
```
---
### [P4484 [BJWC2018] 最长上升子序列](https://www.luogu.com.cn/problem/P4484)
求长为 $n$ 的随机排列最长上升子序列的期望值，对 $998244353$ 取模。
$n\leq28$。
---
设 $f_i$ 是在 $i$ 位置及之前时的最长上升子序列长度（可以不以 $i$ 结尾），那么对其进行差分为 $d_i$，则 $d_i$ 要么是 $0$，要么是 $1$，这是显然的。
既然如此，我们可以状压 $d$ 数组，定义 $dp_{i,d}$ 表示 $1$ 到 $i$ 的排列，$d$ 状态下的方案数，则答案为：
$$\frac{1}{n!}\sum\limits_{d=0}^{2^{n-1}-1}dp_{n,d}\times\mathrm{len}(d)$$
其中 $\mathrm{len}(d)$，表示 $d$ 中二进制位 $1$ 的个数，程序中可以使用`__builtin_popcount`求解。
我们可以钦定从小到大插入 $i$ 到位置 $k$，由于 $i$ 是最大的，故 $d_k=1$。而它后面第一个为 $1$ 的位置 $p$ 就应当 $d_p=0$，因为从这里开始 $i$ 对最长上升子序列的长度就没有影响了。
注意到 $d_1$ 为 $1$，可以少压一位，但时间复杂度 $O(n2^n)$ 依旧很难通过，所以——打表！
可以用[杨表（Young tableau）](https://harris.green/archives/dify-post-1774152478)高效解决，但不在本文讨论范围内，留一道要用杨表的题 [P3774 [CTSC2017] 最长上升子序列](https://www.luogu.com.cn/problem/P3774)，有兴趣可以去了解一下。
```cpp
fac = dp[1][0] = 1;
printf("1,\n");
for(int i = 2;i <= n;++i){
fac = (long long)fac * i % mod;
for(int d = 0;d < (1<<(i-2));++d){
if(!dp[(i^1)&1][d])continue;
dp[i&1][d] = ((long long)dp[i&1][d] + dp[(i^1)&1][d]) % mod; //将i插入到序列首位
for(int j = i - 2;j >= 0;--j){ //第j位
int nd = (d>>j)<<(j+1); //清空后j位
nd |= (1<<j); //第j位设置为1
int t = d&((1<<j)-1);
t = rem(t); //移除后面的最高位
nd |= t;
dp[i&1][nd] = ((long long)dp[i&1][nd] + dp[(i^1)&1][d]) % mod;
}
}
long long ans = 0;
for(int d = 0;d < (1<<(i-1));++d){
ans = (ans + (long long)dp[i&1][d] * (__builtin_popcount(d) + 1) % mod) % mod;
dp[(i^1)&1][d] = 0; //把下一次要用的清空
}
printf("%lld,\n",ans * inv(fac) % mod);
}
```
## 杂项
### [P1410 子序列](https://www.luogu.com.cn/problem/P1410)
给定一个长为偶数 $n$ 的序列 $a_i$，判断其是否能分为两个长为 $\frac{n}{2}$ 的严格上升子序列。
$n\leq2\times10^3$。
---
记以 $a_{i-1}$ 的上升子序列为 $A$，另一个为 $B$。
定义 $f_{i,j}$ 为第 $i$ 个元素，$A$ 的长度为 $j$ 时，$B$ 的最小长度。
当 $a_{i-1}<a_i$ 时，可以将 $a_i$ 接到 $A$ 后面，故 $f_{i,j}=\min\{f_{i,j},f_{i-1,j-1}\}$。
当 $f_{i-1,j-1}<a_i$ 时，可以将 $a_i$ 接到 $B$ 后面，调换 $A,B$，故 $f_{i,i-j+1}=\min\{f_{i,i-j+1},a_{i-1}\}$，因为原来的 $A$ 以 $a_{i-1}$ 结尾。
时间复杂度 $O(n^2)$。
```cpp
for(int i = 0;i <= n;++i){
for(int j = 0;j <= n;++j)f[i][j] = inf;
}
f[0][0] = 0;
for(int i = 1;i <= n;++i){
for(int j = 1;j <= i;++j){
if(a[i - 1] < a[i])f[i][j] = min(f[i][j],f[i - 1][j - 1]);
if(f[i - 1][j - 1] < a[i])f[i][i - j + 1] = min(f[i][i - j + 1],a[i - 1]);
}
}
if(f[n][n / 2] == inf)printf("No!\n");
else printf("Yes!\n");
```
---
### [P4309 [TJOI2013] 最长上升子序列](https://www.luogu.com.cn/problem/P4309)
$n$ 次操作，第 $i$ 次在序列下标 $x_i$ 后插入元素 $i$，求每次插入后的最长上升子序列长度，允许离线。
$n\leq10^5$。
---
假设我们已经拿到最后的数组 $a_i$，鉴于值域就是 $n$，那么我们直接用树状数组就能求出以 $a_i$ 结尾的最长上升子序列长度 $f_i$。令 $g_{a_i}=f_i$，由于后插入的元素一定更大，不会影响之前的最长上升子序列，故第 $i$ 次插入点答案就是 $\max\limits_{j=1}^ig_i$。
考虑怎么得到 $a_i$，如果用 vector 的话是 $O(n^2)$（但是能过），如果用平衡树的话那就有点杀鸡用牛刀了，这里介绍 [@c2020HXW](https://www.luogu.com.cn/user/89561) 的[方法](https://www.luogu.com.cn/article/ujcjc8l9)。
既然按题目顺序很困难，那么我们可以倒着来。首先 $n$ 的位置就是 $x_n+1$，那么推广一下 $i$ 的位置就是第 $x_i+1$ 个**没有被确定**的位置。
那现在就很简单了，这个位置是具有单调性的，直接用树状数组维护，在树状数组上倍增即可。
时间复杂度 $O(n\log n)$。
```cpp
int ask(int idx){
int res = 0,s = 0;
for(int i = 18;i >= 0;--i){
if(s + (1<<i) - t[res + (1<<i)] < idx){
res += (1<<i);
s += (1<<i) - t[res];
}
}
return res + 1;
}
for(int i = n;i;--i){
int x = ask(X[i] + 1);
add(x,1);
a[x] = i;
}
```
---
### [P4310 绝世好题](https://www.luogu.com.cn/problem/P4310)
对长为 $n$ 的序列 $a_i$，求最长的子序列满足相邻项的按位与结果不为 $0$。
$n\leq10^5,a\leq10^9$。
---
定义 $f_b$ 为当前最后一项第 $b$ 位为 $1$ 的最大长度。对于 $a_i$，更新所有 $a_i$ 能接的即可。
时间复杂度 $O(n\log a)$。
```cpp
for(int i = 1;i <= n;++i){
int mx = 0;
for(int b = 0;b <= 30;++b)if((1<<b)&a[i])mx = max(mx,f[b]);
for(int b = 0;b <= 30;++b)if((1<<b)&a[i])f[b] = max(mx + 1,f[b]);
}
```