最长公共子串算法
本文详细介绍了求解最长公共子串问题的三种方法：基本算法、动态规划（DP方案）及后缀数组方法。通过实例讲解了每种方法的具体实现过程，并对各种方法的时间和空间复杂度进行了对比分析。
本文转自： http://zhedahht.blog.163.com/blog/static/254111742007376431815/
题目：如果字符串一的所有字符按其在字符串中的顺序出现在另外一个字符串二中，则字符串一称之为字符串二的子串。注意，并不要求子串（字符串一）的字符必须连续出现在字符串二中。请编写一个函数，输入两个字符串，求它们的最长公共子串，并打印出最长公共子串。
例如：输入两个字符串BDCABA和ABCBDAB，字符串BCBA和BDAB都是是它们的最长公共子串，则输出它们的长度4，并打印任意一个子串。
分析：求最长公共子串（Longest Common Subsequence, LCS）是一道非常经典的动态规划题，因此一些重视算法的公司像MicroStrategy都把它当作面试题。
完整介绍动态规划将需要很长的篇幅，因此我不打算在此全面讨论动态规划相关的概念，只集中对LCS直接相关内容作讨论。如果对动态规划不是很熟悉，请参考相关算法书比如算法讨论。
先介绍LCS问题的性质：记Xm={x0, x1,…xm-1}和Yn={y0,y1,…,yn-1}为两个字符串，而Zk={z0,z1,…zk-1}是它们的LCS，则：
1. 如果xm-1=yn-1，那么zk-1=xm-1=yn-1，并且Zk-1是Xm-1和Yn-1的LCS；
2. 如果xm-1≠yn-1，那么当zk-1≠xm-1时Z是Xm-1和Y的LCS；
3. 如果xm-1≠yn-1，那么当zk-1≠yn-1时Z是Yn-1和X的LCS；
下面简单证明一下这些性质：
1. 如果zk-1≠xm-1，那么我们可以把xm-1（yn-1）加到Z中得到Z’，这样就得到X和Y的一个长度为k+1的公共子串Z’。这就与长度为k的Z是X和Y的LCS相矛盾了。因此一定有zk-1=xm-1=yn-1。
既然zk-1=xm-1=yn-1，那如果我们删除zk-1（xm-1、yn-1）得到的Zk-1，Xm-1和Yn-1，显然Zk-1是Xm-1和Yn-1的一个公共子串，现在我们证明Zk-1是Xm-1和Yn-1的LCS。用反证法不难证明。假设有Xm-1和Yn-1有一个长度超过k-1的公共子串W，那么我们把加到W中得到W’，那W’就是X和Y的公共子串，并且长度超过k，这就和已知条件相矛盾了。
2. 还是用反证法证明。假设Z不是Xm-1和Y的LCS，则存在一个长度超过k的W是Xm-1和Y的LCS，那W肯定也X和Y的公共子串，而已知条件中X和Y的公共子串的最大长度为k。矛盾。
3. 证明同2。
有了上面的性质，我们可以得出如下的思路：求两字符串Xm={x0, x1,…xm-1}和Yn={y0,y1,…,yn-1}的LCS，如果xm-1=yn-1，那么只需求得Xm-1和Yn-1的LCS，并在其后添加xm-1（yn-1）即可；如果xm-1≠yn-1，我们分别求得Xm-1和Y的LCS和Yn-1和X的LCS，并且这两个LCS中较长的一个为X和Y的LCS。
如果我们记字符串Xi和Yj的LCS的长度为c[i,j]，我们可以递归地求c[i,j]：
/ 0 if i<0 or j<0
c[i,j]= c[i-1,j-1]+1 if i,j>=0 and xi=xj
\ max(c[i,j-1],c[i-1,j] if i,j>=0 and xi≠xj
上面的公式用递归函数不难求得。但从前面求Fibonacci第n项(本面试题系列第16题）的分析中我们知道直接递归会有很多重复计算，我们用从底向上循环求解的思路效率更高。
为了能够采用循环求解的思路，我们用一个矩阵（参考代码中的LCS_length）保存下来当前已经计算好了的c[i,j]，当后面的计算需要这些数据时就可以直接从矩阵读取。另外，求取c[i,j]可以从c[i-1,j-1] 、c[i,j-1]或者c[i-1,j]三个方向计算得到，相当于在矩阵LCS_length中是从c[i-1,j-1]，c[i,j-1]或者c[i-1,j]的某一个各自移动到c[i,j]，因此在矩阵中有三种不同的移动方向：向左、向上和向左上方，其中只有向左上方移动时才表明找到LCS中的一个字符。于是我们需要用另外一个矩阵（参考代码中的LCS_direction）保存移动的方向。
参考代码如下：
#include "string.h"
// directions of LCS generation
enum decreaseDir {kInit = 0, kLeft, kUp, kLeftUp};
/////////////////////////////////////////////////////////////////////////////
// Get the length of two strings' LCSs, and print one of the LCSs
// Input: pStr1 - the first string
// pStr2 - the second string
// Output: the length of two strings' LCSs
/////////////////////////////////////////////////////////////////////////////
int LCS(char* pStr1, char* pStr2)
{
if(!pStr1 || !pStr2)
return 0;
size_t length1 = strlen(pStr1);
size_t length2 = strlen(pStr2);
if(!length1 || !length2)
return 0;
size_t i, j;
// initiate the length matrix
int **LCS_length;
LCS_length = (int**)(new int[length1]);
for(i = 0; i < length1; ++ i)
LCS_length[i] = (int*)new int[length2];
for(i = 0; i < length1; ++ i)
for(j = 0; j < length2; ++ j)
LCS_length[i][j] = 0;
// initiate the direction matrix
int **LCS_direction;
LCS_direction = (int**)(new int[length1]);
for( i = 0; i < length1; ++ i)
LCS_direction[i] = (int*)new int[length2];
for(i = 0; i < length1; ++ i)
for(j = 0; j < length2; ++ j)
LCS_direction[i][j] = kInit;
for(i = 0; i < length1; ++ i)
{
for(j = 0; j < length2; ++ j)
{
if(i == 0 || j == 0)
{
if(pStr1[i] == pStr2[j])
{
LCS_length[i][j] = 1;
LCS_direction[i][j] = kLeftUp;
}
else
LCS_length[i][j] = 0;
}
// a char of LCS is found,
// it comes from the left up entry in the direction matrix
else if(pStr1[i] == pStr2[j])
{
LCS_length[i][j] = LCS_length[i - 1][j - 1] + 1;
LCS_direction[i][j] = kLeftUp;
}
// it comes from the up entry in the direction matrix
else if(LCS_length[i - 1][j] > LCS_length[i][j - 1])
{
LCS_length[i][j] = LCS_length[i - 1][j];
LCS_direction[i][j] = kUp;
}
// it comes from the left entry in the direction matrix
else
{
LCS_length[i][j] = LCS_length[i][j - 1];
LCS_direction[i][j] = kLeft;
}
}
}
LCS_Print(LCS_direction, pStr1, pStr2, length1 - 1, length2 - 1);
return LCS_length[length1 - 1][length2 - 1];
}
/////////////////////////////////////////////////////////////////////////////
// Print a LCS for two strings
// Input: LCS_direction - a 2d matrix which records the direction of
// LCS generation
// pStr1 - the first string
// pStr2 - the second string
// row - the row index in the matrix LCS_direction
// col - the column index in the matrix LCS_direction
/////////////////////////////////////////////////////////////////////////////
void LCS_Print(int **LCS_direction,
char* pStr1, char* pStr2,
size_t row, size_t col)
{
if(pStr1 == NULL || pStr2 == NULL)
return;
size_t length1 = strlen(pStr1);
size_t length2 = strlen(pStr2);
if(length1 == 0 || length2 == 0 || !(row < length1 && col < length2))
return;
// kLeftUp implies a char in the LCS is found
if(LCS_direction[row][col] == kLeftUp)
{
if(row > 0 && col > 0)
LCS_Print(LCS_direction, pStr1, pStr2, row - 1, col - 1);
// print the char
printf("%c", pStr1[row]);
}
else if(LCS_direction[row][col] == kLeft)
{
// move to the left entry in the direction matrix
if(col > 0)
LCS_Print(LCS_direction, pStr1, pStr2, row, col - 1);
}
else if(LCS_direction[row][col] == kUp)
{
// move to the up entry in the direction matrix
if(row > 0)
LCS_Print(LCS_direction, pStr1, pStr2, row - 1, col);
}
}
扩展：如果题目改成求两个字符串的最长公共子字符串，应该怎么求？子字符串的定义和子串的定义类似，但要求是连续分布在其他字符串中。比如输入两个字符串BDCABA和ABCBDAB的最长公共字符串有BD和AB，它们的长度都是2。
扩展题目答案：
这个LCS跟前面说的最长公共子序列的LCS不一样，不过也算是LCS的一个变体，在LCS中，子序列是不必要求连续的，而子串则是“连续”的。即：
题：给定两个字符串X，Y，求二者最长的公共子串，例如X=[aaaba]，Y=[abaa]。二者的最长公共子串为[aba]，长度为3。
本节给出三种不同的实现方式，并对比分析每种方法的复杂度，内容如下：
==基本算法==
==DP方案==
==后缀数组==
==各方法复杂度分析==
==================================
基本算法
其实对于最长公共子串，还是比较简单易想的，因为子串是连续的，这就方便了很多。最直接的方法就是用X每个子串与Y的每个子串做对比，求出最长的公共子串。代码如下：
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
|
/* 最长公共子串 Longest Common Substring */
int maxlen; /* 记录最大公共子串长度 */
int maxindex; /* 记录最大公共子串在串1的起始位置 */
void outputLCS(char * X); /* 输出LCS */
/* 最长公共子串 基本算法 */
int comlen(char * p, char * q)
{
int len = 0;
while(*p && *q && *p++ == *q++)
{
++len;
}
return len;
}
void LCS_base(char * X, int xlen, char * Y, int ylen)
{
for(int i = 0; i < xlen; ++i)
{
for(int j = 0; j < ylen; ++j)
{
int len = comlen(&X[i],&Y[j]);
if(len > maxlen)
{
maxlen = len;
maxindex = i;
}
}
}
outputLCS(X);
}
|
==================================
DP方案
既然最长公共子串是最长公共子序列的变体，那么最长公共子串是不是也可以用动态规划来求解呢？
我们还是像之前一样“从后向前”考虑是否能分解这个问题，在最大子数组和中，我们也说过，对于数组问题，可以考虑“如何将arr[0,...i]的问题转为求解arr[0,...i-1]的问题”，类似最长公共子序列的分析，这里，我们使用dp[i][j]表示 以x[i]和y[j]结尾的最长公共子串的长度，因为要求子串连续，所以对于X[i]与Y[j]来讲，它们要么与之前的公共子串构成新的公共子串；要么就是不构成公共子串。故状态转移方程
- X[i] == Y[j]，dp[i][j] = dp[i-1][j-1] + 1
- X[i] != Y[j]，dp[i][j] = 0
对于初始化，i==0或者j==0，如果X[i] == Y[j]，dp[i][j] = 1；否则dp[i][j] = 0。
代码如下：
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
|
/* 最长公共子串 DP */
int dp[30][30];
void LCS_dp(char * X, int xlen, char * Y, int ylen)
{
maxlen = maxindex = 0;
for(int i = 0; i < xlen; ++i)
{
for(int j = 0; j < ylen; ++j)
{
if(X[i] == Y[j])
{
if(i && j)
{
dp[i][j] = dp[i-1][j-1] + 1;
}
if(i == 0 || j == 0)
{
dp[i][j] = 1;
}
if(dp[i][j] > maxlen)
{
maxlen = dp[i][j];
maxindex = i + 1 - maxlen;
}
}
}
}
outputLCS(X);
}
|
==================================
后缀数组
前面提过后缀数组的基本定义，与子串有关，可以尝试这方面思路。由于后缀数组最典型的是寻找一个字符串的重复子串，所以，对于两个字符串，我们可以将其连接到一起，如果某一个子串s是它们的公共子串，则s一定会在连接后字符串后缀数组中出现两次，这样就将最长公共子串转成最长重复子串的问题了，这里的后缀数组我们使用基本的实现方式。
值得一提的是，在找到两个重复子串时，不一定就是X与Y的公共子串，也可能是X或Y的自身重复子串，故在连接时候我们在X后面插入一个特殊字符‘#’，即连接后为X#Y。这样一来，只有找到的两个重复子串恰好有一个在#的前面，这两个重复子串才是X与Y的公共子串。
代码如下：
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
|
/* 最长公共子串 后缀数组 */
char * suff[100];
int pstrcmp(const void *p, const void *q)
{
return strcmp(*(char**)p,*(char**)q);
}
int comlen_suff(char * p, char * q)
{
int len = 0;
while(*p && *q && *p++ == *q++)
{
++len;
if(*p == '#' || *q == '#')
{
return len;
}
}
return 0;
}
void LCS_suffix(char * X, int xlen, char * Y, int ylen)
{
int suf_index = maxlen = maxindex = 0;
int len_suff = xlen + ylen + 1;
char * arr = new char [len_suff + 1]; /* 将X和Y连接到一起 */
strcpy(arr,X);
arr[xlen] = '#';
strcpy(arr + xlen + 1, Y);
for(int i = 0; i < len_suff; ++i) /* 初始化后缀数组 */
{
suff[i] = & arr[i];
}
qsort(suff, len_suff, sizeof(char *), pstrcmp);
for(int i = 0; i < len_suff-1; ++i)
{
int len = comlen_suff(suff[i],suff[i+1]);
if(len > maxlen)
{
maxlen = len;
suf_index = i;
}
}
outputLCS(suff[suf_index]);
}
|
下面给出三种实现方案所用到的打印输出LCS程序以及测试用例程序。
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
|
/* 输出LCS
* 在后缀数组方法中，maxindex=0
* 因为传进来的就是后缀数组suff[]，从0打印maxlen个字符
*/
void outputLCS(char * X)
{
if(maxlen == 0)
{
printf("NULL LCS\n");
return;
}
printf("The len of LCS is %d\n",maxlen);
int i = maxindex;
while(maxlen--)
{
printf("%c",X[i++]);
}
printf("\n");
}
void main()
{
char X[] = "aaaba";
char Y[] = "abaa";
/* 基本算法 */
LCS_base(X,strlen(X),Y,strlen(Y));
/* DP算法 */
LCS_base(X,strlen(X),Y,strlen(Y));
/* 后缀数组方法 */
LCS_suffix(X,strlen(X),Y,strlen(Y));
}
|
==================================
各方案复杂度对比
设字符串X的长度为m，Y的长度为n，最长公共子串长度为l。
对于基本算法，X的子串（m个）和Y的子串（n个）一一对比，最坏情况下，复杂度为O(m*n*l)，空间复杂度为O(1)。
对于DP算法，由于自底向上构建最优子问题的解，时间复杂度为O(m*n)；空间复杂度为O(m*n)，当然这里是可以使用滚动数组来优化空间的，滚动数组在动态规划基础回顾中多次提到。
对于后缀数组方法，连接到一起并初始化后缀数组的时间复杂度为O(m+n)，对后缀数组的字符串排序，由于后缀数组有m+n个后缀子串，子串间比较，故复杂度为O((m+n)*l*lg(m+n))，求得最长子串遍历后缀数组，复杂度为O(m+n)，所以总的时间复杂度为O((m+n)*l*lg(m+n))，空间复杂度为O(m+n)。
总的来说使用后缀数组对数据做一些“预处理”，在效率上还是能提升不少的。
本节相关代码可以到这里下载。
（全文完）
勘误：（感谢 @jiawei）
后缀数组方法comlen_suff的实现逻辑上有错误，要保证找到公共子串后二者只有一个#号，这样才表示来自不同的字符串，勘误代码如下：
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
|
int comlen_suff(char * p, char * q)
{
int len = 0;
while(*p && *q && *p++ == *q++)
{
++len;
if(*p == '#' || *q == '#')
{
break;
}
}
int count = 0;
while(*p)
{
if(*p++ == '#')
{
++count;
break;
}
}
while(*q)
{
if(*q++ == '#')
{
++count;
break;
}
}
if(count == 1)
return len;
return 0;
}
|
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
