最长公共子序列

或许大家不知道除了$n^{2}$的LCS做法，如果不知道，那看这篇文章就对了；

看完本篇文章，你会学会的技能：

给定两个字符串$a$,$b$,它们的长度均小于等于**100000**的最长公共子序列(但要做到这一点需要毒瘤的代码实现技巧和恶魔般的内存优化，如果你是手残+脑残，只能解决$n<=50000$的数据范围)

需要读者的前置技能：$n^{2}$的dp求LCS的做法；

这篇文章并没有深入讲解普通的LCS的dp解法，所以食用本篇文章前请先自行学习$n^{2}$的LCSdp算法；

我们先不看dp方程，想一想能否有不是二维dp的做法（时间复杂度可以大于$n^{2}$，只要不再是$n^{2}$的dp）：

有!

假设现在有两个字符串a，b；它们的长度分别是$n$,$m$;

对于a中的每个字符$a[k]$，找到$a[k]$在b中的每个出现的位置(下标从1开始)，并把这些位置倒序排序，存在vector $c[k]$中；

然后从$1~k$把$c[]$连起来，具体理解可以参照下面的例子；

设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$的最长公共子序列；

下面是感性理解时间：

因为是严格上升子序列，每个数最多只会出现一次，所以最长上升子序列中出现的数两两不同；

由定义可知：b[c[x][y]]=a[x]；

为了得到最长的公共子序列，我们需要找到最多的$pair(x,c[x][y])$，使得满足所有选出的数对，第一维的数两两不同，第二维的数也是两两不同的，并且当第一维是上升的时候，第二维同样上升；

这不就是最长上升子序列吗，而至于最长上升子序列LIS，我们可以通过树状数组做到$O(nlogn)$求出；

至于时间复杂度和空间复杂度嘛，是个玄学；

时间复杂度上限是$O(n*m*log(n*m))$,空间复杂度上限是$O(n*m)$;

比如说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)$

所以如果数据是用脚造的随机出的，并且字典集(即所有出现过的字符的种类数目)比较大，这种方法就会跑得比较快；

但是这仅仅是在随机数据下字典集比较大的情况下使用的方法，那么如果字典集比较小（更准确的说是不那么大）呢？怎么办？难道这只能是严格$O(n^{2})$的dp了吗；

这回我们不从定义入手，而是从dp方程入手；

众所周知，二进制的位运算是很快的，可以考虑从这方面入手；

$left{egin{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} ight.$

(注意我的方程的定义：是$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]$;

显然的，答案就是$sum_{k=0}^{n-1} M[m][k]$

首先，我们将a串反转。然后我们建立一个bitset类型的数组$t[]$,对于每个在$b$中出现过的字符，利用二进制表示出在a中所有的出现位置；

比如说：对于a:$acabcabdb$ b: $abaabb$

将a串反转得到bdbacbaca;

$'a'=000100101$

$'b'=101001000$

接下来我们观察一下得到的$M$数组：

 其中：红框框圈起来的部分便是$M$数组

对与每一行，我们可以将其看成一个二进制bitset类型；

设第i行的二进制表示是s[i]，我们来说一下s[i]的求法；

用第s[1]和s[2]的转化为例：

$s[1]=00000,100,1$。我们将每一位1与前面的部分断开，也就是说：$s[3]=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$

我们这样做的意义便是尽可能的找到一个1的个数大于等于$s[1]$中1的个数（也就是最长公共子序列的部分匹配），且尽可能的对后面的影响小(所以选择每一块内靠右的)；

我们发现，只有$s[i]$中分出的块中存在全是0的块才有可能使得匹配数+1

然而怎么找呢？暴力的话复杂度还不如$n^{2}$,这时候想到了我们的好朋友：二进制运算；

首先我们通过$xor$运算得到$x[i]=s[i-1] or t['b[i]']$。

然后我们可以想办法把每一块内最右侧的1(包含1)以及后面的0都变成1，块内其余的数都变成0，得到一个bitset类型的tmp。接着将tmp和x做$and$运算，就是取到了每一块中最右侧的1的位置的新的bitset； 

那么怎样把每一块内最右侧的1及其后面的0都变成1，块内其余的数都变成0呢？

首先，我们将$x-((s[i-1]<<1)|1)$的结果与$x[i]$进行异或运算就好了；

对于$(s[i]<<1)|1$的解释：将每一块内最右侧一个1变成0，并将这个1右侧的所有0变成1，其余的位置不变；

对于与$x[i]$进行异或的解释：将每一块内最右侧的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(n*m)$(但是基本都是位运算，常数很小)，空间复杂度$O(n*m)$(用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<7011> s[50001];
bitset<7011> t[26];
char a[50010],b[50010];
int n,m;
void operator-=(bitset<7011> &a,bitset<7011> 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,7000)){
            if(a[j]-'a'==i) t[i].set(j,1);
        }        
    }
    //cout<<t['G'-'A']
    bitset<7011> 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
*/

 (早期码风若毒瘤请见谅)；

但是我们发现它不香，因为bitset类型不存在减法，所以需要自己写减法；

而且用bitset常数巨大，这比普通的$O(n^{2})$算法快不了多少(随机情况下是普通写法的$frac{1}{10}$)；

那怎么办呢？我们想到了指令集(指令集是不让用的啦~，但是可以借鉴它的思想)；

首先自己定义bitset类型(反正都要自己定义减法再多定义几个运算也无所谓的吧)；

我们把若干位封装在一起，得到一个不超过二进制数，将其看作十进制正整数进行二进制运算（and，or，xor，<< 等）;

然后优化就结束啦~我们成功的优化掉了一个巨大的常数(64);

有人或许觉得二进制常数优化优化不了什么，但这个优化与直接用bitset运算的快慢就如同FFT迭代版与FFT递归版的区别，如同吸了氧气的毒瘤二进制版莫队与普通版莫队的区别；

经实践可知：这种做法完全可以通过$a$,$b$串的长度均在$70000$数据强度时且时限是$1s$的数据点；如果写法更优美一点更毒瘤一点，我们应该可以通过$nleqslant 100000$的数据点；

(内存限制不成问题，观察式子发现，$s[i]$只和$s[i-1]$和t[]有关，在字符集较小的情况下完全可以滚动$s[]$通过,字符集大的情况下一边滚动$s[]$,一边在$t[]$中只记录1的位置有哪些(内存均摊复杂度$O(m)$)，在使用$t[]$的时候还原成二进制就好了)；

下面放上本人蒟蒻的代码(没有滚动等内存优化，因为这篇代码只是给大家提供一个自定义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
", count(Ans));
    return 0;
}