最长公共子串(Longest Common Substring)是一个非常经典的问题,它的基本描述为“给定两个字符串,求出它们之间最长的相同子字符串(要求连续)的长度”。例如以下两个字符串S和T的的最长公共子串“howmuchiloveyoumydearmother”,长度为27.
对于求多个字符串的LCS,byvoid 在“ 最长公共子串问题的后缀数组解法”一文中给出了解法,将求N个串的最长公共子串转化为 求一些后缀的最长公共前缀的最大值,这些后缀应分属于N个串。设N个串分别为S1,S2,S3,...,SN,具体方法如下:
S="yeshowmuchiloveyoumydearmotherreallyicannotbelieveit" T="yeaphowmuchiloveyoumydearmother"求N个最长为L的字符串的的LCS的方法大致可分为以下几类
- 枚举法显然是简单但极端低效的算法,改进一些的算法是用一个串的每个后缀对其他所有串进行部分匹配,用KMP算法,时间复杂度为O(NL2).
- 动态规划解法,时间复杂度为O(N2)
- 后缀数组与高度数组解法,利用二分查找技术,时间复杂度为O(NLlogL)。
- 广义后缀树方法,时间复杂度为可以达到 O(NL)。
动态规划
作为最长公共子序列问题的简化,常规解法为动态规划(DP),参考代码如下。
int LCS(string S,string T){
if(S==""||T=="")return 0;
if(S.length()>T.length())swap(S,T);
int ls=S.length(),MAX=-1;
int now=0,former=1;
vector<vector<int>> memo(2,vector<int>(ls));
for(auto t:T){
memo[now][0]=(t==S[0]);
for(int i=1;i<ls;++i){
if(S[i]!=t)memo[now][i]=0;
else MAX=max(MAX,memo[now][i]=memo[former][i-1]+1);
}
swap(now,former);
}
return MAX;
}
后缀数组算法
后缀数组SA[i]表示排名第i的后缀的位置,高度数组(LCP数组)Height[i]表示后缀SA[i]和SA[i-1]的最长公共前缀(Longest Common Prefix, LCP),简记为Height[i]=LCP(SA[i],SA[i-1]),详细介绍可参见后缀数组(Suffix Arrary) .求两个字符串S和T的最长公共子串,可以构造一个新串S+T,然后建立后缀数组,求出LCP数组。可以知道,最长公共前缀一定是出现在后缀数组中相邻的两个子串里的(否则越离越远,只可能更小)。利用这一点,我们找出满足sa[i]和sa[i+1]分属于S1,S2的串,然后对应的LCP数组中的最大值就是答案。下面给出了求两个字符串的LCS的代码,利用了后缀数组(Suffix Arrary)一文中给出的后缀树类(SA),SA有两个成员变量sa和height,分别为构造字符串S的后缀数组和高度数组。
int LCS(const string &S,const string &T) {
string R=S+T;
int len=R.size(),ls=S.size(),MAX=0;
//构造后缀数组,第二个参数为构造SA的算法(如Prefix Doubling, DC3, SA-IS等)
SA mySA(R,&SA::prefixDoubling);
for(int i=0;i<len;++i) {
if((mySA.sa[i]<ls)!=(mySA.sa[i+1]<ls))MAX=max(MAX,mySA.height[i]);
}
return MAX;
}
对于求多个字符串的LCS,byvoid 在“ 最长公共子串问题的后缀数组解法”一文中给出了解法,将求N个串的最长公共子串转化为 求一些后缀的最长公共前缀的最大值,这些后缀应分属于N个串。设N个串分别为S1,S2,S3,...,SN,具体方法如下:
- 首先建立一个串S,把这N个串用不同的分隔符连接起来。S=S1[P1]S2[P2]S3...SN-1[PN-1]SN,P1,P2,...PN-1应为不同的N-1个不在字符集中的字符,作为分隔符(后面会解释为什么)。
- 接下来,求出字符串S的后缀数组和Height数组,可以用倍增算法,或DC3算法。
- 然后二分枚举答案A,假设N个串可以有长度为A的公共字串,并对A的可行性进行验证。如果验证A可行,A'(A'<A)也一定可行,尝试增大A,反之尝试缩小A。
- 最终可以取得A的最大可行值,就是这N个串的最长公共子串的长度,可以证明,尝试次数是O(logL)的。
后缀树算法
先将串 S 构造为 SAM ,然后用 T 按如下规则去跑自动机。- 用一个变量 lcs 记录当前的最长公共子串,初始化为0。
- 设当前状态结点为 p,要匹配的字符为 c,若 go[c] 中有边,说明能够转移状态,则转移并 lcs++;
- 若不能转移则将状态移动到 p 的 par ,如果仍然不能转移则重复该过程直到 p 回到根节点,并将 lcs 置为 0;
- 如果在上一个过程中进入了能够转移的状态,则设 lcs 为当前状态的 val。
为什么失配后要移向 par 呢?因为在状态 p 上失配说明该状态的 [min,max] 所表示字符串都不是 B 中的子串,但是比它们短的后缀仍有可能是 B 的子串,而 par 指针恰好指向了该状态的后缀。
#define N 200005
struct Node { int f,next[26],l; }G[N];
class SAM {
char *st;
int n,len,last,cnt;
public:
int MAX;
SAM(string S):len(1),cnt(0),last(1),MAX(0) {
st=(char*)malloc(S.size());
for(int i=0;i<S.size();++i)add(S[i]-'a');
last=1;
}
void add(int ch) {
int i,x,p;
G[++len].l=G[last].l+1;
for(i=last;i&&!G[i].next[ch];i=G[i].f)G[i].next[ch]=len;
p=G[i].next[ch],last=len,x=i;
if(!x){G[1].next[ch]=len,G[len].f=1; return;}
if(G[p].l==G[x].l+1) G[len].f=p;
else {
G[++len].l=G[x].l+1,G[len].f=G[p].f,G[p].f=G[len-1].f=len;
for(i=0;i<26;i++)G[len].next[i]=G[p].next[i];
for(i=x;i&&G[i].next[ch]==p;i=G[i].f)G[i].next[ch]=len;
}
}
void explore(int ch) {
if(G[last].next[ch])++cnt;
else {
for(;last&&!G[last].next[ch];last=G[last].f);
if(!last){last=1,cnt=0; return;}
cnt=G[last].l+1;
}
last=G[last].next[ch],MAX=max(MAX,cnt);
}
};
int LCS(string S,string T) {
SAM mySAM(S);
for(int i=0;i<T.size();++i)mySAM.explore(T[i]-'a');
return mySAM.MAX;
}