求所有最大公共子序列的算法实现-最大公共子序列

一：LCS解析

首先看下什么是子序列？定义就不写了，直接举例一目了然。如对于字符串：“student”，那么su，sud，sudt等都是它的子序列。它可以是连续的也可以不连续出现，如果是连续的出现，比如stud，一般称为子序列串，这里我们只讨论子序列。

什么是公共子序列？很简单，有两个字符串，如果包含共同的子序列，那么这个子序列就被称为公共子序列了。如“student”和“shade”的公共子序列就有“s”或者“sd”或者“sde”等。而其中最长的子序列就是所谓的最长公共子序列(LCS)。当然，最长公共子序列也许不止一个，比如：“ABCBDAB”和“BDCABA”，它们的LCS为“BCBA”，“BCAB”，“BDAB”。知道了这些概念以后就是如何求LCS的问题了。

通常的算法就是动态规划(DP)。假设现在有两个字符串序列：X={x1,x2,...xi...xm},Y={y1,y2,...yj...yn}。如果我们知道了X={x1,x2,...xi-1}和Y={y1,y2,...yj-1}的***公共子序列L，那么接下来我们可以按递推的方法进行求解：

1)如果xi==yj,那么{L,xi(或yj)}就是新的LCS了，其长度也是len(L)+1。这个好理解，即序列{Xi,Yj}的***解是由{Xi-1,Yj-1}求得的。

2)如果xiyj,那么可以转换为求两种情况下的LCS。

A: X={x1,x2,...xi}与Y={y1,y2,...yj-1}的LCS，假设为L1

B: X={x1,x2,...xi-1}与Y={y1,y2,...yj}的LCS，假设为L2

那么xiyj时的LCS=max{L1,L2},即取***值。同样，实际上序列{Xi,Yj-1}和{Xi-1,Yj}都可以由{Xi-1,Yj-1}的***解求得。

怎么样，是不是觉得这种方法很熟悉？当前问题的***解总是包含了一个同样具有***解的子问题，这就是典型的DP求解方法。好了，直接给出上面文字描述解法中求LCS长度的公式：

这里用一个二维数组存储LCS的长度信息，i，j分别表示两个字符串序列的下标值。这是求***公共子序列长度的方法，如果要打印出***公共子序列怎么办？我们还需要另外一个二维数组来保存求解过程中的路径信息，方便***进行路径回溯，找到LCS。如果看着很含糊，我下面给出其实现过程。

二：DP实现

很多博文上面都有，基本上是用两个二维数组c[m][n]和b[m][n]，一个用来存储子字符串的LCS长度，一个用来存储路径方向，然后回溯。

其中二维数组b[i][j]的取值为1或2或3，其中取值为1时说明此时xi=yj，c[i][j]=c[i-1][j-1]+1。如果将二维数组看成一个矩阵，那么此时代表了一个从左上角到右下角的路径。如果取值为2，说明xi≠yj,且c[i][j]=c[i-1][j]，代表了一个从上到下的路径，同理取值为3代表一个从左到右的路径。

***我们可以根据c[m][n]的值知道***公共子序列的长度。然后根据b[i][j]回溯，可以打印一条LCS。其中b[i][j]=1的坐标点对应的字符同时在两个序列中出现，所以依次回溯这个二维数组就可以找到LCS了。这里给出实现代码：

[[98078]]View Code

我们给出一个测试：

1     char str1[MAX_LEN] = "BADCDCBA"; 2     char str2[MAX_LEN] = "ABCDCDAB"; 3 
4     GetLCSLen(str1, str2, C, B, str1len+1, str2len+1); 5     TraceBack(str1, B, str1len+1, str2len+1);

详细的代码见文章结束处给出的链接。本测试output:BDCDB

问题：上面的方法中，我们没有单独考虑c[i-1][j]==c[i][j-1]的情况，所以在回溯的时候打印的字符只是其中一条***公共子序列，如果存在多条公共子序列的情况下。怎么解决？我们对b[i][j]二维数组的取值添加一种可能，等于4，这代表了我们说的这种多支情况，那么回溯的时候我们可以根据这个信息打印更多可能的选择。这个过程就不写代码了，其实很简单，以下面的路径回溯图举例，你从（8,8)点开始按b[i][j]的值指示的方向回溯，把所有的路径遍历一遍，如果是能达到起点（0,0)的路径，就是LCS了，有多少条打印多少条。可是，

又出现问题了：你发现没有，在回溯路径的时候，如果采用一般的全搜索，会进行了很多无用功。即重复了很多，且会遍历了一些无效路径，因为这些路径最终不会到达终点(0,0)，比如节点(6,3),(7,2),(8,1)。因此加大算法复杂度和时间消耗。那么如何解决？看下面的这个方法，正式进入本文正题。

路径回溯图：

加入状态4后的状态图：

三：算法改进

上面提到路径回溯过程中，一般的方法就是遍历所有可能的路径，但是一些不可能构成***公共子序列的跳跃点我们也会去计算。这里先解释下什么叫跳跃点，就是导致公共子序列长度发生变化的节点，即b[i][j]=1对应的节点(i,j)。Ok，接下来的问题是，如何不去考虑这些无效跳跃点，降低算法复杂度？参考论文里提出了这样一种方法：矩形搜索。

首先构造两个栈数据结构store和print。故名思议，一个用来储存节点，一个用来打印节点。栈的定义为：

 1 #define MAX_STACK_SIZE 1024
 2 typedef struct _Element  3 {  4     int nlcslen;  5     int nrow;  6     int ncolumn;  7 }Element;  8 typedef struct _Stack  9 { 10     int top; 11     Element data[MAX_STACK_SIZE]; 12 }Stack;

栈使用数组实现，并有一个指向顶点的下标值top。为了初始化需要，先构造了一个虚拟节点virtualnode，指向节点(m,n)的右下角，即(m+1,n+1),这个节点的LCS的长度假设为***公共子序列长度+1。将虚拟节点压入栈store，然后然后执行下面的算法：

1)栈store为空吗？是的话退出。

2)否则从store栈顶弹出节点。

3)如果这个节点为边界元素(行或列的小标为1),则将此节点压入栈print，打印栈print里面的所有节点(除virtualnode外)。查看此时store栈顶节点的LCS长度，并以这个长度为参考，弹出print栈里面所有LCS长度小于等于这个参考值的节点，跳转到第1步。

4)如果不是边界元素，那么以该节点的左上节点(i-1,j-1)为出发点，沿着该出发点指示的方向，找到***个跳跃点e1(即b[i][j]==1的点)。途中碰到分支节点(b[i][j]==4的点)时，沿着其上节点继续探索。

5)找到***个跳跃点以后，重新回到第4步的出发点，沿着该节点指示的方向，找到第二个跳跃点e2。途中碰到分支节点(b[i][j]==4的点)时，沿着其左节点继续探索

6)如果e1和e2节点为同一节点，将该节点压入store栈，回到步骤1)。

7）如果不为同一节点，在e1和e2构成的矩形范围内，搜索出所有的跳跃点，并全部压入store栈，回到步骤1)。

不明白？不要紧，我们结合上面的矩阵图一步步按照算法来，看看到底是如何计算的：

***步：压入虚拟节点6(9,9)到栈store，这里6表示这个节点的LCS长度，（9,9)表示坐标值。

第二步：store栈不为空，则弹出store栈顶，压入print栈，这时候的两个栈的状态如下面的左图。沿出发点(8,8)出发，这是个分支节点，因为b[8][8]==4,所以选择向上走，搜索到e1跳跃点(7,8)，搜索路径为：(8,8)->(7,8)。然后回到(8,8)找e2点，这时选择向左走，找到e2跳跃点(8,7)。这两个跳跃点不同，所以以e1，e2为对角线方向围成的矩形内搜索所有跳跃点，这里只有e1,和e2本身两个节点，然后将它们压入栈store。此时两个栈的状态见下面的有图。蓝色底的节点表示有store栈弹出然后压到print栈，绿色底表示新压入到store栈的跳跃点，下面所有的图都这样表示。

第三步：弹出5(7,8)到print栈，搜索到新的两跳跃节点。

第四步：

第五步：

第六步：

第七步：关键步骤来了，因为此时从store栈弹出的节点是边界元素1(1,2)，所以我们打印print栈的所有元素(红色字体节点)，而这些元素恰好构成了一个最长公共子序列(自己揣摩一下)。打印完了以后，我们要对print栈进行清理，去除不需要的节点，按照步骤2，此时store栈顶节点的LCS为1，所以print栈中弹出的节点只有1(1,2)。弹完以后，print栈的状态如图所示。虚线框节点表示已弹出，下同。