最长公共子串问题的后缀数组解法
http://www.byvoid.com/blog/lcs-suffix-array/
http://www.cppblog.com/superKiki/archive/2010/10/29/131781.aspx?opt=admin
[最长公共子串]
最长公共子串(Longest Common Substring ,简称LCS)问题，是指求给定的一组字符串长度最大的共有的子串的问题。例如字符串”abcb”,”bca”,”acbc”的LCS就是”bc”。
求多串的LCS，显然穷举法是极端低效的算法。改进一些的算法是用一个串的每个后缀对其他所有串进行部分匹配，用KMP算法，时间复杂度为O(N*L^2)，其中N为字符串个数，L为每个串的长度。更优秀的有广义后缀树的方法，时间可以达到 O(N*L)。本文介绍一种基于后缀数组的LCS解法，利用二分查找技术，时间复杂度可以达到O(N*L*logL)。
[最长公共子串问题的后缀数组解法]
关于后缀数组的构建方法以及Height数组的性质，本文不再具体介绍，可以参阅IOI国家集训队2004年论文《后缀数组》(许智磊)和IOI国家集训队2009年论文《后缀数组——处理字符串的有力工具》(罗穗骞)。
回顾一下后缀数组，SA[i]表示排名第i的后缀的位置，Height[i]表示后缀SA[i]和SA[i-1]的最长公共前缀(Longest Common Prefix,LCP)，简记为Height[i]=LCP(SA[i],SA[i-1])。连续的一段后缀SA[i..j]的最长公共前缀，就是H[i-1..j]的最小值，即LCP(SA[i..j])=Min(H[i-1..j])。
求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)的。
于是问题就集中到了，如何验证给定的长度A是否为可行解。方法是，找出在Height数组中找出连续的一段Height[i..j]，使得i<=k<=j均满足Height[k]>=A，并且i-1<=k<=j中，SA[k]分属于原有N个串S1..SN。如果能找到这样的一段，那么A就是可行解，否则A不是可行解。
具体查找i..j时，可以先从前到后枚举i的位置，如果发现Height[i]>=A，则开始从i向后枚举j的位置，直到找到了Height[j+1]<A，判断[i..j]这个区间内SA是否分属于S1..SN。如果满足，则A为可行解，然后直接返回，否则令i=j+1继续向后枚举。S中每个字符被访问了O(1)次，S的长度为N*L+N-1，所以验证的时间复杂度为O(N*L)。
到这里，我们就可以理解为什么分隔符P1..PN-1必须是不同的N-1个不在字符集中的字符了，因为这样才能保证S的后缀的公共前缀不会跨出一个原有串的范围。
后缀数组是一种处理字符串的强大的数据结构，配合LCP函数与Height数组的性质，后缀数组更是如虎添翼。利用后缀数组，容易地求出了多个串的LCS，而且时空复杂度也相当优秀了。虽然比起后缀树的解法有所不如，但其简明的思路和容易编程的特点却在实际的应用中并不输于后缀树。
* Problem: POI2000 pow
* Author: Guo Jiabao
* Time: 2009.4.16 21:00
* State: Solved
* Memo: 多串最长公共子串 后缀数组 二分答案
*/
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cstring>
const int MAXL=10011,MAXN=6;
using namespace std;
struct SuffixArray
{
struct RadixElement
{
int id,k[2];
}RE[MAXL],RT[MAXL];
int N,A[MAXL],SA[MAXL],Rank[MAXL],Height[MAXL],C[MAXL];
void RadixSort()
{
int i,y;
for (y=1;y>=0;y--)
{
memset(C,0,sizeof(C));
for (i=1;i<=N;i++) C[RE[i].k[y]]++;
for (i=1;i<MAXL;i++) C[i]+=C[i-1];
for (i=N;i>=1;i--) RT[C[RE[i].k[y]]--]=RE[i];
for (i=1;i<=N;i++) RE[i]=RT[i];
}
for (i=1;i<=N;i++)
{
Rank[ RE[i].id ]=Rank[ RE[i-1].id ];
if (RE[i].k[0]!=RE[i-1].k[0] || RE[i].k[1]!=RE[i-1].k[1])
Rank[ RE[i].id ]++;
}
}
void CalcSA()
{
int i,k;
RE[0].k[0]=-1;
for (i=1;i<=N;i++)
RE[i].id=i,RE[i].k[0]=A[i],RE[i].k[1]=0;
RadixSort();
for (k=1;k+1<=N;k*=2)
{
for (i=1;i<=N;i++)
RE[i].id=i,RE[i].k[0]=Rank[i],RE[i].k[1]=i+k<=N?Rank[i+k]:0;
RadixSort();
}
for (i=1;i<=N;i++)
SA[ Rank[i] ]=i;
}
void CalcHeight()
{
int i,k,h=0;
for (i=1;i<=N;i++)
{
if (Rank[i]==1)
h=0;
else
{
k=SA[Rank[i]-1];
if (--h<0) h=0;
for (;A[i+h]==A[k+h];h++);
}
Height[Rank[i]]=h;
}
}
}SA;
int N,Ans,Bel[MAXL];
char S[MAXL];
void init()
{
int i;
freopen("pow.in","r",stdin);
freopen("pow.out","w",stdout);
scanf("%d",&N);
SA.N=0;
for (i=1;i<=N;i++)
{
scanf("%s",S);
for (char *p=S;*p;p++)
{
SA.A[++SA.N]=*p-'a'+1;
Bel[SA.N]=i;
}
if (i<N)
SA.A[++SA.N]=30+i;
}
}
bool check(int A)
{
int i,j,k;
bool ba[MAXN];
for (i=1;i<=SA.N;i++)
{
if (SA.Height[i]>=A)
{
for (j=i;SA.Height[j]>=A && j<=SA.N;j++);
j--;
memset(ba,0,sizeof(ba));
for (k=i-1;k<=j;k++)
ba[Bel[SA.SA[k]]]=true;
for (k=1;ba[k] && k<=N;k++);
if (k==N+1)
return true;
i=j;
}
}
return false;
}
void solve()
{
int a,b,m;
SA.CalcSA();
SA.CalcHeight();
a=0;b=SA.N;
while (a+1<b)
{
m=(a+b)/2;
if (check(m))
a=m;
else
b=m-1;
}
if (check(b))
a=b;
Ans=a;
}
int main()
{
init();
solve();
printf("%dn",Ans);
return 0;
}
