题解 SP12076 【LCS0 - Longest Common Subsequence】

· · 题解

最长公共子序列的优化

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的题解(没有滚动数组)：

#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类型的模板)：

#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]石子合并

其实更主要的还是这种二进制优化的思想，它可以被应用到很多意想不到的地方。