最長公共子串

最長公共子串(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年論文《字尾陣列》(許智磊)

回顧一下字尾陣列，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應為不同

接下來，求出字串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;
}