【算法】LCS最长公共子序列

LCS—Longest Common Subsequence

最长公共子序列。一个序列，如果是两个或多个已知序列的子序列，且是所有子序列中最长的，则为最长公共子序列。

#include <stdio.h>#include <string.h>#include <stdlib.h>#define SIZE 100int main(){char A[SIZE],B[SIZE];int i,j,n;printf("输入字符串A:");char *p1 = &(A[1]);gets(p1);A[0] = '0';printf("输入字符串B:");char *p2 = &(B[1]);gets(p2);B[0] = '0';const int len_A = strlen(A)-1;const int len_B = strlen(B)-1;/*数组c用来存放LCS的长度*/int **c;c = (int **)malloc(sizeof(int*)*(len_A+1));for(i=0;i<len_A+1;i++)c[i]=(int *)malloc(sizeof(int)*(len_B+1));/*b用来指明LCS的路径*/char **b;b = (char **)malloc(sizeof(char*)*(len_A+1));for(i=0;i<len_A+1;i++)b[i]=(char *)malloc(sizeof(char)*(len_B+1));char *result;result = (char *)malloc(sizeof(char *)*(len_A+1));for(i=0;i<=len_A;i++){c[i][0] = 0;b[i][0] = '0';}for(j=0;j<=len_B;j++){c[0][j] = 0;b[0][j] = '0';}for(i=1;i<=len_A;i++)for(j=1;j<=len_B;j++){if(A[i] == B[j]){c[i][j] = c[i-1][j-1] + 1;b[i][j] = 'b';}else if(c[i-1][j] >= c[i][j-1]){c[i][j] = c[i-1][j];b[i][j] = 'u';}else{c[i][j] = c[i][j-1];b[i][j] = 'l';}}for(i=0;i<=len_A;i++){for(j=0;j<=len_B;j++)printf("%c  ",b[i][j]);printf("\n");}i = len_A;j = len_B;n = 0;while(b[i][j] != '0'){switch (b[i][j]){case 'u':i--;break;case 'l':j--;break;case 'b':result[n++] = A[i];i--;j--;break;}}free(c);free(b);result[n] = '\0';for(i=n-1;i>=0;i--)printf("%c",result[i]);}

算法分析：

最优子结构：设 X = {x1,x2,......,xm}和Y = {y1,y2,......,yn}为两个序列，并设Z = {z1,z2,......,zk}为X和Y的任意一个LCS。

1)如果xm=yn，那么zk = xm = yn而且Zk-1是Xm-1和Yn-1的一个LCS。

2)如果xm≠yn，那么zk≠xm蕴含Z是Xm-1和Y的一个LCS。

3)如果xm≠yn，那么zk≠yn蕴含Z是X和Yn-1的一个LCS。

用c[i,j]表示序列Xi和Yj的一个LCS的长度。那么有：

c[i,j] = 

case i = 0,j = 0 : c[i,j] = 0;

case i,j > 0;xi = yj:  c[i,j] = c[i-1,j-1] + 1;

case i,j > 0;xi ≠ yj:  c[i,j] = max(c[i,j-1],c[i-1,j])

伪代码：

LCS_LENGTH(X,Y)m ← length[X]n ← length[Y]for i ← 1 to m        do c[i,0] ← 0for j ← 1 to n        do c[0,j] ← 0for i ← 1 to m        do for j ← 1 to n                     do if xi = yj                                then c[i,j] ← c[i-1,j-1] + 1                                         b[i,j] ← "↖"                                else if c[i-1,j] ≥ c[i,j-1]                                            then c[i,j] ← c[i-1,j]                                                     b[i,j] ← "↑"                                           else c[i,j] ← c[i,j-1]                                                    b[i,j] ← "←"return c and b

流程图：

其中数组c中，c[i,j]表示Xi和Yj的一个LCS的长度，刚刚已经说过了，数组b则表示找到最优解的路径。下面举例，

X = {A,B,C,B,D,A,B}

Y = {B,D,C,A,B,A}

m ← length[X]n ← length[Y]for i ← 1 to m        do c[i,0] ← 0for j ← 1 to n        do c[0,j] ← 0

对c表初始化，填0

for i ← 1 to m        do for j ← 1 to n                     do if xi = yj                                then c[i,j] ← c[i-1,j-1] + 1                                         b[i,j] ← "↖"                                else if c[i-1,j] ≥ c[i,j-1]                                            then c[i,j] ← c[i-1,j]                                                     b[i,j] ← "↑"                                           else c[i,j] ← c[i,j-1]                                                    b[i,j] ← "←"

这里有两层for循环，外层循环length_A次，内层length_B次。
逐步分析for循环：

此时，表更新为：

i=2再循环

表c变化：

通过两个j的循环可以看出，每次填入c[i,j]的值，都是查询之前填入的，也就是查表，同时不修改原先的值。通过最优子结构来获得最优值，也就是说，当子结构是最优时，那么通过cut-and-paste，可以得出，最优的解。这就是动态规划的思想。

通过上面两步，我们可以验证一下，虽然j只取了2次，但是我们可以验证X={A,B} Y={B,D,C,A,B,A}的LCS。

如上图所示，LCS为AB，虽然这个结果很特别，X就是Y 的子集，但是还是可以看出动态规划得到最优解的思路。