最近看了很多關於LCS(Longest common subsequence problem,最長公共子序列)的文章，大部分問題都只是求出最大公共子序列的長度，或者打印處其中的任意一個最大子序列即可，但是如何快速的打印出所有的最大長度子序列？這個問題好像看到的不多。本文給出了傳統的DP(dynamic programming，動態規划)算法進行求解的過程，並用c語言實現。另外參考一篇論文實現了其中的一種打印所有最大公共子序列的算法，這個算法比起傳統的算法而言，時間復雜度大大降低.
一：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了。這里給出實現代碼：
View Code
1 #include <stdio.h> 2 #include <stdlib.h> 3 #include <string.h> 4 #include "stack.h" 5 #define MAX_LEN 1024 6 typedef int **Matrix; 7 8 void GetLCSLen(char *str1, char *str2, Matrix pc, Matrix pb, int nrow, int ncolumn); 9 Matrix GreateMatrix(int nrow, int ncolumn); 10 void DeleteMatrix(Matrix p, int nrow, int ncolumn); 11 void TraceBack(char *str1, Matrix pb, int nrow, int ncolumn); 12 13 void GetLCSLen(char *str1, char *str2, Matrix pc, Matrix pb, int nrow, int ncolumn) 14 { 15 int i,j; 16 /************initial the edge***************/ 17 for(i=0; i<nrow; i++) 18 { 19 pc[i][0] = 0; 20 pb[i][0] = 0; 21 } 22 for(j=0; j<ncolumn; j++) 23 { 24 pc[0][j] = 0; 25 pb[0][j] = 0; 26 } 27 /************DP*****************************/ 28 for(i=1; i<nrow; i++) 29 { 30 for(j=1; j<ncolumn; j++) 31 { 32 if(str1[i-1] == str2[j-1]) 33 { 34 pc[i][j] = pc[i-1][j-1] + 1;//由左上節點轉移而來 35 pb[i][j] = 1;//標記為1 36 } 37 else if(pc[i-1][j] >= pc[i][j-1]) 38 { 39 pc[i][j] = pc[i-1][j];//由上節點轉移而來 40 pb[i][j] = 2;//標記為2 41 } 42 else 43 { 44 pc[i][j] = pc[i][j-1];//由左節點轉移而來 45 pb[i][j] = 3;//標記為2 46 } 47 } 48 } 49 } 50 void TraceBack(char *str1, Matrix pb, int nrow, int ncolumn) 51 { 52 int ntemp; 53 if(str1 == NULL || pb == NULL) 54 return; 55 if(nrow == 0 || ncolumn == 0) 56 return; 57 ntemp = pb[nrow-1][ncolumn-1]; 58 switch(ntemp) 59 { 60 case 1: 61 printf("locate:(%d,%d),%4c\n", nrow-1, ncolumn-1, str1[nrow-2]);//打印公共字符，這里下標是nrow-2，因為矩陣的坐標值(i,j)比字符串的實際下標大1 62 TraceBack(str1, pb, nrow-1, ncolumn-1);//向左上角遞歸 63 break; 64 case 2: 65 TraceBack(str1, pb, nrow-1, ncolumn);//向上方向遞歸 66 break; 67 case 3: 68 TraceBack(str1, pb, nrow, ncolumn-1);//向左方向遞歸 69 break; 70 default: 71 break; 72 } 73 }
我們給出一個測試：
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棧的狀態如圖所示。虛線框節點表示已彈出，下同。
第八步：繼續彈出store棧頂，發現又是邊界元素，繼續壓入print棧，打印print棧。清理print棧。
第九步：清理完后，繼續步驟2.
好了，下面的過程就是重復進行上面的這些步驟了，自己動手畫一下就一目了然了。
四：代碼實現
用c語言實現關鍵代碼：
View Code
1 void TraceBack2(char *str1, Matrix pc, Matrix pb, int nrow, int ncolumn); 2 void PrintStack(Stack *ps, char *str1, int len1); 3 void SearchE(Matrix pb, int curposx, int curposy, int *eposx, int *eposy, int ntype); 4 5 void PrintStack(Stack *ps, char *str1, int len1) 6 { 7 if(ps == NULL || str1 == NULL) 8 return; 9 int ntemp = ps->top; 10 int index = -1; 11 while((index = ps->data[ntemp].nrow) <= len1) 12 { 13 ntemp--; 14 printf("%2c",str1[index-1]); 15 } 16 printf("\n"); 17 } 18 19 void TraceBack2(char *str1, Matrix pc, Matrix pb, int nrow, int ncolumn) 20 { 21 if(str1 == NULL || pc == NULL || pb == NULL) 22 return; 23 24 Stack store, print;//構造兩個棧store，print 25 Element storetop;//store棧的棧頂節點 26 Element element;//臨時變量 27 Element virtualnode;//虛擬節點 28 int ntoplen;//保存store棧頂節點的LCS長度 29 int ex1,ey1,ex2,ey2;//矩形搜索的兩個節點的坐標 30 int i,j; 31 32 InitialStack(&store);//初始化 33 InitialStack(&print); 34 virtualnode = CreateElement(pc[nrow-1][ncolumn-1]+1, nrow, ncolumn); 35 Push(&store, &virtualnode);//壓入虛擬節點到store 36 37 while(!IsEmpty(&store)) 38 { 39 Pop(&store, &storetop);//從棧頂取出一個節點 40 if(storetop.nrow == 1 || storetop.ncolumn == 1)//如果是邊界節點 41 { 42 Push(&print, &storetop); 43 PrintStack(&print, str1, nrow-1);//打印print棧里面除虛擬節點之外的所有節點 44 GetTop(&store, &element); 45 ntoplen = element.nlcslen;//當前store的棧頂節點的LCS長度 46 47 /**********彈出print棧中所有LCS長度小於等於ntoplen的節點**************/ 48 while(GetTop(&print, &element) && element.nlcslen<=ntoplen) 49 { 50 Pop(&print, &element); 51 } 52 } 53 else 54 { 55 Push(&print, &storetop); 56 SearchE(pb, storetop.nrow-1, storetop.ncolumn-1, &ex1, &ey1, 0); 57 SearchE(pb, storetop.nrow-1, storetop.ncolumn-1, &ex2, &ey2, 1/*also other value is ok*/); 58 59 if(ex1 == ex2 && ey1 ==ey2) 60 { 61 element = CreateElement(pc[ex1][ey1], ex1, ey1); 62 Push(&store,&element);//壓入store棧，回到步驟2 63 } 64 else 65 { 66 for(i=ey2; i<=ey1; i++) 67 for(j=ex1; j<=ex2; j++) 68 { 69 if(pb[i][j] == 1) 70 { 71 element = CreateElement(pc[i][j], i, j); 72 Push(&store, &element); 73 } 74 } 75 } 76 } 77 78 } 79 } 80 void SearchE(Matrix pb, int curposx, int curposy, int *eposx, int *eposy, int ntype) 81 { 82 switch(pb[curposx][curposy]) 83 { 84 case 1: 85 *eposx = curposx; 86 *eposy = curposy; 87 return; 88 case 2: 89 SearchE(pb, curposx-1, curposy, eposx, eposy, ntype); 90 break; 91 case 3: 92 SearchE(pb, curposx, curposy-1, eposx, eposy, ntype); 93 break; 94 case 4: 95 if(ntype == 0) 96 SearchE(pb, curposx-1, curposy, eposx, eposy, ntype);//搜索e1點，如過碰到分叉點，向上繼續搜索 97 else 98 SearchE(pb, curposx, curposy-1, eposx, eposy, ntype);//搜索e2點，如過碰到分叉點，向左繼續搜索 99 break; 100 } 101 }
同樣的測試，這種算法能打印出全部的最長公共子序列。
五：總結
該算法能將傳統算法的指數級復雜度降低到max{O(mn),O(ck)},k為最大公共子序列的個數。詳細證明見論文。因為用兩個棧存儲了所有有效跳躍點，使得許多重復比較被忽略。棧的有序性也能巧妙的得到任意一條最大公共子序列。
本算法的全部實現代碼下載路徑：http://pan.baidu.com/share/link?shareid=125428&uk=285541510
歡迎討論。
參考論文：
《利用矩陣搜索求所有最長公共子序列的算法》，宮潔卿，安徽工程科技學院學報，vol23,No.4,Dec.,2008