题解 SP12076 【LCS0 - Longest Common Subsequence】
神之右大臣
2020-03-03 09:07:28
# 最长公共子序列的优化
## 0.前言：
或许大家不知道除了$n^{2}$的LCS做法，如果不知道，那看这篇文章就对了；
需要读者的前置技能：
1. **$n^{2}$的dp求LCS的做法**；
2. **二进制运算的一些小技巧**；
看完本篇文章，你会学会的技能：
**给定两个字符串$a$,$b$(没有字符集的限制),它们的长度均小于等于100000，求它们的最长公共子序列**
请注意：要做到这一点需要毒瘤的代码实现技巧和恶魔般的内存优化，如果你是手残+脑残，只能解决$n<=50000$的数据范围。
------------
## 1. 首先，我们来看看一种可能被卡的玄学优化：
我们先不看dp方程，想一想能否有**不是二维dp**的做法（时间复杂度可以大于$n^{2}$，只要不再是$n^{2}$的dp）：
显然有，因为洛谷的LCS模板就是这样的，虽然它有字符集的限制。
假设现在有两个字符串a，b；它们的长度分别是$n$,$m$;
对于a中的每个字符$a[k]$，找到$a[k]$在b中的每个出现的位置(下标从1开始)，并把这些位置倒序排序，存在vector $c[k]$中；
然后从$1$到$k$把$c[]$连起来，求出它们的最长**严格上升**子序列。
------------
### 1.1 让我们快乐的举个例子：
设a:$acabcabdb$ b: $abaabb$
$c[1]={4,3,1}$
$c[2]={}$
$c[3]={4,3,1}$
$c[4]={6,5,2}$
$c[5]={}$
$c[6]={4,3,1}$
$......$
$c[9]={6,5,2}$
然后将$c[]$拼起来,得到：{4,3,1, ,4,3,1,6 5 2, ,4,3,1,......,6,5,2}
我们将拼起来的数组中所有空元素删除，对其求最长**严格上升**子序列，就是$a$和$b$的最长公共子序列；
------------
### 1.2 简述该方法的正确性：
因为是严格上升子序列，每个数最多只会出现一次，所以最长上升子序列中出现的数两两不同；
由定义可知：b[c[x][y]]=a[x]；
为了得到最长的公共子序列，我们需要找到最多的$pair(x,c[x][y])$，使得满足所有选出的数对，第一维的数两两不同，第二维的数也是两两不同的，并且当第一维是上升的时候，第二维同样上升；
这不就是最长上升子序列吗，而至于最长上升子序列LIS，我们可以通过树状数组做到$O(nlogn)$求出；
------------
### 1.3 时间复杂度论证：
至于时间复杂度和空间复杂度嘛，是个玄学；
1.上限：时间复杂度上限是$O(nmlog(nm))$,空间复杂度上限是$O(nm)$;
比如说a串是$aaaaaaaaa......$,b串也是$aaaaaa.........$;那么c数组会很大很大，时间复杂度甚至没有n^n的dp好；
但时间复杂度下限是$O(n*log(n))$,空间复杂度下限是$O(n)$;
比如说a串是一个排列，b串也是一个排列，那么c数组会很小很小，时间复杂度产生了巨大的优化；
对于随机数据，设它的字典集(即所有出现过的字符的种类数目)大小为$w$,那么时间复杂度从$O(n^{2})$变成了$O(\frac{n^{2}} {w} *logn)$
因此求两个排列的LCS便是严格的$O(nlogn)$复杂度。
重要的事情说三遍：
**1.如果$w>logn$,那么就是正优化，否则就是负优化。**
**2.如果$w>logn$,那么就是正优化，否则就是负优化。**
**3.如果$w>logn$,那么就是正优化，否则就是负优化。**
所以如果数据是~~用脚造的~~**随机出的(这很重要)**，并且字典集比较大，这种方法就会跑得比较快；
### ~~好了，完结撒花~~~
Q：欸？不对不对，这糊弄三岁小孩纸呢？这这不是在字符集有限制的情况下才能产生正优化吗？举报举报！出现标题党！与标题严重不符！
A：别着急，接着向下看，因为简单的优化应该放在前面不是吗？
------------
## 2. 二进制优化：
刚才的算法仅仅是在随机数据下字典集比较大的情况下使用的方法，那么如果字典集比较小（更准确的说是不那么大）呢？怎么办？难道只能是严格$O(n^{2})$的dp了吗；
众所周知，二进制的位运算是很快的，可以考虑从这方面结合dp方程入手。
让我们先来看一看dp方程：
$\left\{\begin{matrix}f[i][j]=max(f[i-1][j],f[i][j-1])&\\ f[i][j]=max(f[i][j],f[i-1][j-1]+1)&,b[i]=a[j]\end{matrix}\right.$
**特别注意：我的方程的定义：是$b[i]=a[j]$而不是$a[i]=b[j]$,并且矩阵的下标是从0开始，字符串$a，b$的下标也是从0开始。**
我们发现，这个$f$所表示的矩阵对于任意的位置$f[i][j]$,它最多比$f[i-1][j-1]$大1；(比较显然吧~)
那么我们可以构造一个01矩阵$M$，对于任意的$f[i][j]$,$f[i][j]=\sum_{k=0}^{j} M[i][k]$;
显然的，答案就是$f[m-1][n-1]=\sum_{k=0}^{n-1} M[m-1][k]$
------------
### 2.1 得到$M$数组的方法：
#### 2.1.1.首先，我们将a串反转(下标也反转过来)。然后我们建立一个bitset类型的数组$t$,对于每个在$b$中出现过的字符，利用二进制表示(可以使用bitset)出在a中所有的出现位置；
比如说：对于a串:$acabcabdb$ b串: $abaabb$
将a串反转得到$bdbacbaca$;
$'a'=000100101$
$'b'=101001000$
接下来我们观察一下得到的$M$数组：
其中：红框框圈起来的部分便是$M$数组:
#### 2.1.2.对与每一行，我们可以将其看成一个二进制bitset类型。
设第i行的二进制表示是s[i]，我们来说一下s[i]的求法。
用第s[1]和s[2]的转化为例：
$s[1]=000001001$。**我们将每一位1与前面的部分断开（这就是'块'定义）**，也就是说：$s[1]=00000,100,1$
然后对于第2行的字符b[2]:$a$,将$t['a']$按照$s[1]$的断点断开，也就是说：$t['a']=00010,010,1$
我们将$s$和$t$进行$or$运算，得到：$00010,110,1$
对于每块，我们选取最右侧的1的位置，把$s[2]$的这一位赋值成1。
也就是说：$s[2]=00010,010,1$
#### 2.1.3.对于每块(块的定义在上面的例子中提到了)，我们选取最右侧的1的位置，把$s[2]$的这一位赋值成1；
------------
### 2.2 这种做法的含义：
我们这样做的意义便是尽可能的找到一个1的个数**大于等于**$s[1]$中1的个数（也就是最长公共子序列的部分匹配），且尽可能的对后面的影响小(所以选择每一块内靠右的)；
我们发现，只有当$s[i]$中所有分出的块中存在**全是0**的块(也就是最靠左的块全是0)才有可能使得匹配数+1
------------
### 2.3 利用二进制运算加速递推$M$数组：
然而怎么找每一块中最靠右的1呢？暴力的话复杂度还不如$n^{2}$。
这时候想到了我们的好朋友：二进制运算；
首先我们通过$xor$运算得到$x[i]=s[i-1]$ $or$ $t['b[i]']$。
然后我们可以想办法把每一块内**最右侧的1(包含1)以及后面的0都变成1，块内其余的数都变成0**，得到一个bitset类型的tmp。接着将tmp和$x[i]$做$and$运算，就是取到了每一块中最右侧的1的位置的新的bitset；
那么怎样把每一块内最右侧的1及其后面的0都变成1，块内其余的数都变成0呢？
**我们将$x-((s[i-1]<<1)|1)$的结果与$x[i]$进行异或运算就好了；**
对于$x-(s[i]<<1)|1$的解释：将每一块内最右侧一个1变成0，并将这个1右侧的所有0变成1，其余的位置不变；
对于与$x[i]$进x-行异或的解释：将每一块内最右侧的1到块的结束都变成1，其余的都变成0；
综上所述，得到式子：$s[i]=(s[i-1] or t['b[i]'])and((s[i-1]| or ['b[i]']) xor ((s[i-1] or t['b[i]'])-((s[i-1]<<1)|1))))$;
时间复杂度是$O(nm)$(但是基本都是位运算，常数很小)，空间复杂度$O(nm)$(用bitset类型，实际用到的内存很小)；
那么我们可以轻松地写出一个求LCS的题解(没有滚动数组)：
```cpp
#include <bits/stdc++.h>
#define inc(i,a,b) for(register int i=a;i<=b;i++)
#define dec(i,a,b) for(register int i=a;i>=b;i--)
using namespace std;
bitset<20011> s[20001];
bitset<20011> t[26];
char a[20010],b[20010];
int n,m;
void operator-=(bitset<20011> &a,bitset<20011> b) {
unsigned int last = 0;
for (int i=0;i<7000;i++)
if (a[i]>=b[i]+last)
a[i]=a[i]-(b[i]+last),last=0;
else
a[i]=a[i]-(b[i]+last),last=1;
}
int main()
{
scanf("%s %s",a,b);
n=strlen(a); m=strlen(b);
inc(i,0,25){
t[i].reset();
inc(j,0,min(n-1,20000)){
if(a[j]-'a'==i) t[i].set(j,1);
}
}
//cout<<t['G'-'A']
bitset<20011> tmp=t[b[0]-'a'],tmp2,tmp3=tmp,tmp4;tmp2.set(0,1); tmp4.set(0,1);
tmp3-=tmp2;
s[0]=tmp&(tmp3^tmp);
inc(i,1,m-1){
tmp=s[i-1]|t[b[i]-'a'];
tmp3=tmp;
tmp2=(s[i-1]<<1)|tmp4;
tmp3-=tmp2;
s[i]=tmp&(tmp3^tmp);
}
dec(i,m-1,0){
dec(j,n-1,0){
cout<<s[i][j]<<" ";
}
cout<<b[i]<<" "<<i;
cout<<endl;
}
dec(i,n-1,0) cout<<a[i]<<" ";
cout<<endl;
dec(i,n-1,0) cout<<i<<" ";
// cout<<(s[2]|t[1])<<endl;
cout<<(int)(s[m-1].count());
}
/*
acabcabdb
abaabb
*/
```
(早期码风若毒瘤请见谅)；
------------
### 2.4利用指令集 _**思想**_ 来优化二进制运算
但是我们发现它不香，**因为bitset类型不存在减法**(如果你直接写减法编译会错的)，所以需要自己写减法；
而且用bitset常数巨大，这比普通的$O(n^{2})$算法快不了多少(随机情况下是普通写法的$\frac{1}{10}$)；
那怎么办呢？我们想到了指令集(指令集是不让用的啦~，但是可以借鉴它的思想，就是那个把几个数封装成一块的思想)；
首先自己定义bitset类型(反正都要自己定义减法再多定义几个运算也无所谓的吧)；
我们**把若干位封装在一起**（不足补0），得到一个不超过long long的二进制数，将其看作十进制正整数进行二进制运算（and，or，xor，<< 等）;
### 让我们来举个例子：
假设bitset类型是：10001011001010。按照3位封装的话，就是010,001,011,001,010，那么新数组就是：2，1，3，1，2。完全可以$O(1)$的进行整数long long类型的运算。
------------
然后优化就结束啦~我们成功的优化掉了一个巨大的常数(64);
有人或许觉得二进制常数优化优化不了什么，但这个优化与直接用bitset运算的快慢就如同FFT迭代版与FFT递归版的区别，如同吸了氧气的毒瘤二进制版莫队与普通版莫队的区别；
经实践可知：这种做法完全可以通过$a$,$b$串的长度均在$70000$数据强度时且时限是$1s$的数据点。如果写法更优美一点更毒瘤一点，我们可以通过$n\leqslant 100000$的数据点； (好了，看到这就证明我不是标题党了吧。)
------------
A：又在逗小孩纸了，时间是降下来了，但内存咋办啊？
Q：内存限制不成问题，观察式子发现，$s[i]$只和$s[i-1]$和$t$有关，在字符集较小的情况下完全可以滚动$s$通过,字符集大的情况下一边滚动$s$,一边在$t$中只记录1的位置有哪些(内存均摊复杂度$O(m)$)，在使用$t$的时候还原成二进制就好了)；
下面放上本人蒟蒻的代码(没有滚动等内存优化，在$1s$内只能通过$n<=70000$的数据。因为这篇代码主要是给大家提供一个**自定义带有指令集思想的bitset类型**的模板)：
```cpp
#include <bits/stdc++.h>
using namespace std;
const int N = 50010, M = 2200;
int n, m;
char a[N], b[N];
struct BitSet {
unsigned int a[M];
} Ans, X, Y, A[N];
BitSet operator&(BitSet a, BitSet b) {
for (int i = 0; i < M; i++) a.a[i] &= b.a[i];
return a;
}
BitSet operator^(BitSet a, BitSet b) {
for (int i = 0; i < M; i++) a.a[i] ^= b.a[i];
return a;
}
BitSet operator|(BitSet a, BitSet b) {
for (int i = 0; i < M; i++) a.a[i] |= b.a[i];
return a;
}
void operator-=(BitSet &a, BitSet b) {
unsigned int last = 0;
for (int i = 0; i < M; i++)
if (a.a[i] >= b.a[i] + last)
a.a[i] -= b.a[i] + last, last = 0;
else
a.a[i] -= b.a[i] + last, last = 1;
}
void operator<<=(BitSet &a, BitSet b) // a=b*2+1
{
unsigned int last = 1;
for (int i = 0; i < M; i++) {
unsigned int x = b.a[i] >> 31;
a.a[i] = (b.a[i] << 1 | last);
last = x;
}
}
void Insert(BitSet &A, int x) { A.a[x >> 5] |= 1u << (x & 31); }
int count(BitSet A) {
int ans = 0;
for (int i = 0; i < M; i++) ans += __builtin_popcount(A.a[i]);
return ans;
}
int main() {
scanf("%s %s",a+1,b+1);
n=strlen(a+1); m=strlen(b+1);
for (int i = 1; i <= n; i++) {
Insert(A[a[i]-'a'], i);
}
for (int i = 1; i <= m; i++) {
Y <<= Ans;
Ans = Ans | A[b[i]-'a'];
X = Ans;
X -= Y;
Ans = Ans & (Ans ^ X);
for(int j=0;j<=n;j++){
cout<<Ans.a[j]<<" ";
}
cout<<endl;
}
printf("%d\n", count(Ans));
return 0;
}
```
## 总结：
或许在省选中乃至NOI、IOI中不会考察LCS这种基础算法，但谁又说得准呢。之前有一年的省选不是还考四边形不等式优化都过不去的石子合并了吗？需要用什么加西亚瓦克斯算法(这是什么鬼？？？)，详见**洛谷P5569 [SDOI2008]石子合并**
其实更主要的还是这种二进制优化的思想，它可以被应用到很多意想不到的地方。