后綴數組：倍增法和DC3的簡單理解

一些定義：設字符串S的長度為n,S[0~n-1]。

子串：設0<=i<=j<=n-1，那么由S的第i到第j個字符組成的串為它的子串S[i,j]。

后綴：設0<=i<=n-1，那么子串S[i,n-1]稱作它的后綴，用Suffix[i]表示。

串比較：對於兩個串S1，S2，設長度分別為n1,n2。若存在一個位置i，使得對於0<=j<i滿足S1[j]=S2[j]且S1[i]<S2[i]，那么我們稱S1<S2。如果S1是S2的一個前綴，那么也有S1<S2。兩個串相等當且僅當長度相同且所有位置的字母都相同。所以，對於S的任意兩個不同的后綴Suffix[i],Suffix[j] ,它們一定是不相等的，因為它們的長度都不同。

后綴數組：設我們用數組sa表示S的后綴數組，0<=sa[i]<=n-1，表示將S的n個后綴從小到大排序后，排名第i的后綴的位置是sa[i]。

名次數組：設我們用數組rank表示S的名次數組，0<=rank[i]<=n-1，表示將S的n個后綴從小到大排序后，后綴Suffix[i]的排名是rank[i]。很明顯，sa[rank[i]]=i。

現在我們的問題是，給出一個字符串S，長度為n，S[0~n-1]，字符集大小為m。求出S的后綴數組sa。有兩種方法計算這個sa，倍增法和DC3法。我們設m<=n（一般情況也是這樣子的吧。。）。那么倍增法的時間復雜度是O(nlog(n))，DC3的時間復雜度是O(n)。兩個方法的空間復雜度都是O(n)。

 

倍增法

倍增法的思路是:

（1）首先計算S[0],S[1],...,S[n-1]的排名(注意這個單個字符的排序)。比如，對於aabaaaab，排序后為：1,1,2,1,1,1,1,2

（2）計算子串S[0,1],S[1,2],S[2,3],...,S[n-2,n-1],S[n-1,null] 的排名（注意最后一個的第二個字符為空），由於我們知道了單個字符的排名， 那么每個子串可以用一個二元組來表示，比如S[0,1]={1,1}，S[1,2]={1,2},S[2,3]={2,1}，等等，也就是aa,ab,ba,aa,aa,aa,ab,b$（我們用$表示空）的排名，排序后為：1,2,4,1,1,1,2,3

（3）計算子串S[0,1,2,3],S[1,2,3,4],S[2,3,4,5],...,S[n-4,n-3,n-2,n-1],S[n-3,n-2,n-1,$],S[n-2,n-1,$,$],S[n-1,$,$,$]。方法與上面相同。依次類推，每次使用兩個2^(x-1)長度的子串來計算2^x次方長度的子串的排名，直到某一次排序后n個數字各不相同。最后，對於串aabaaaab，如下所示

在實現的時候，一般在串的最后加一個空字符。這個字符比其他任何字符都小。下面是一份實現的程序

class SuffixArray
{
private:
    static const int N=100005; /**字符串最大長度**/
    int wa[N],wb[N],wd[N],r[N];

    bool isSame(int *r,int a,int b,int len)
    {
        return r[a]==r[b]&&r[a+len]==r[b+len];
    }
    void da(int *r,int *sa,int n,int m)
    {
        int *x=wa,*y=wb,*t;
        for(int i=0;i<m;i++) wd[i]=0;
        for(int i=0;i<n;i++) wd[x[i]=r[i]]++;
        for(int i=1;i<m;i++) wd[i]+=wd[i-1];
        for(int i=n-1;i>=0;i--) sa[--wd[x[i]]]=i;
        /**基數排序計算長度為1的子串的排名
           相同的 越靠前 排名越小**/
        for(int j=1,p=1;p<n;j<<=1,m=p)
        {
            p=0;
            for(int i=n-j;i<=n-1;i++) y[p++]=i;
            for(int i=0;i<n;i++) if(sa[i]>=j) y[p++]=sa[i]-j;
            /**y[i]表示對於組成2^j的所有子串的二元組
               {pi,qi}來說，第二關鍵字即qi排名為i的位置為y[i] **/


            for(int i=0;i<m;i++) wd[i]=0;
            for(int i=0;i<n;i++) wd[x[i]]++;
            for(int i=1;i<m;i++) wd[i]+=wd[i-1];
            for(int i=n-1;i>=0;i--) sa[--wd[x[y[i]]]]=y[i];
             /**這里倒着枚舉 當兩個位置的y[i]和y[j]對應的x
                相同時，后面的排名大，因為它的第二關鍵字
                即y的排名大 而x在外面也決定了排名以第一
                關鍵字為主 **/
            t=x;x=y;y=t;p=1;
            x[sa[0]]=0;
            for(int i=1;i<n;i++) x[sa[i]]=isSame(y,sa[i-1],sa[i],j)?p-1:p++;
        }
    }
public:
    /**字符串S,長度n,S[0,n-1]，為方便，我們假設它只包含小寫字母
       最后的后綴數組存儲在sa[1~n]中  0<=sa[i]<=n-1
    **/
    void calSuffixArray(char *S,int n,int *sa)
    {
        for(int i=0;i<n;i++) r[i]=S[i]-'a'+1;
        r[n++]=0;
        da(r,sa,n,27);
    }
};

DC3算法

（1）將所有的后綴分成兩部分，一部分是模3不等於0的，比如Suffix[1],Suffix[2],Suffix[4],Suffix[7]等，第二部分是模3等於0的后綴，Suffix[0],Suffix[3]等。首先計算第一部分每個后綴的排名（計算的時候假設沒有第二部分的這些后綴）。方法是將Suffix[1]和Suffix[2]連起來（ 連起來之前要把Suffix[1]和Suffix[2]的長度都變為3的倍數，如果不是，就在后面補上字符集中沒有且小於字符集中所有字符的字符，比如0）。對於串S=aabaaaaba，Suffix[1]=abaaaaba,Suffix[2]=baaaaba，分別補成3的倍數，Suffix[1]=abaaaaba0,Suffix[2]=baaaaba00，最后拼成的串為abaaaaba0baaaaba00 。如下所示

然后從前向后，每三個字符一組，即aba,aaa,ba0,baa,aab,a00，我們發現，他們分別是Suffix[1],Suffix[4],Suffix[7],Suffix[2],Suffix[5],Suffix[8] 的前3個字母。我們求出這六個的排名為4,2,5,6,3,1（注意，如果排序后還有相同的數字，也就是還有兩個相同的串，比如3,2,4,5,2,1，那么要繼續求，因為兩個2之后的數字4大於1，所以第二的位置的2代表的后綴大於第5個位置的2代表的后綴。其實這個問題跟剛才的問題是相同的，所以可以遞歸求）。這樣，我們最后得到了所有模3不等於0的位置的后綴的排名。

（2）計算模3等於0的位置的排名。這些位置的后綴，可以看做一個字符加上某個第一部分的一個后綴，這也很容易通過一次基數排序（就像倍增法的二元組一樣）求得。對於上面的串，模3為0的后綴的排名為Suffix[9]<Suffix[3]<Suffix[0]<Suffix[6]

（3）合並第一部分和第二部分的排名。注意，上面求出的第一部分第二部分的排名都沒有考慮另外一部分。合並的時候我們需要比較第一部分的某個后綴和第二部分的某個后綴。分兩種情況。第一種是比較Suffix[3*i] 和Suffix[3*j+1]，我們把它們看做：

Suffix[3*i]=S[3*i]+Suffix[3*i+1]

Suffix[3*j+1]=S[3*j+1]+Suffix[3*j+2]

Suffix[3*i+1]和Suffix[3*j+2]可以直接比較，因為它們都屬於第一部分，而S[3*i]和S[3*j+1]也可以直接比較；

第二種情況是Suffix[3*i] 和Suffix[3*j+2]，把它們看做是

Suffix[3*i]=S[3*i]+S[3*i+1]+Suffix[3*i+2]

Suffix[3*j+2]=S[3*j+2]+S[3*j+3]+Suffix[3*(j+1)+1]

Suffix[3*i+2]和Suffix[3*(j+1)+1]可以直接比較，它們都屬於第二部分。而  前面是兩個單個字符，可以直接比較。這樣，就可以合並所有的后綴得到答案。

class SuffixArray
{
private:
    static const int N=100005; /**字符串最大長度**/
    int wa[N*3],wb[N*3],wv[N*3],ws[N*3];
    int r[N],sa[N];

    int c0(int *r,int a,int b)
    {
        return r[a]==r[b]&&r[a+1]==r[b+1]&&r[a+2]==r[b+2];
    }

    int c12(int k,int *r,int a,int b)
    {
        if(k==2) return r[a]<r[b]||r[a]==r[b]&&c12(1,r,a+1,b+1);
        else return r[a]<r[b]||r[a]==r[b]&&wv[a+1]<wv[b+1];
    }

    void sort(int *r,int *a,int *b,int n,int m)
    {
        for(int i=0;i<n;i++) wv[i]=r[a[i]];
        for(int i=0;i<m;i++) ws[i]=0;
        for(int i=0;i<n;i++) ws[wv[i]]++;
        for(int i=1;i<m;i++) ws[i]+=ws[i-1];
        for(int i=n-1;i>=0;i--) b[--ws[wv[i]]]=a[i];
    }

    void dc3(int *r,int *sa,int n,int m)
    {
        #define F(x) ((x)/3+((x)%3==1?0:tb))
        #define G(x) ((x)<tb?(x)*3+1:((x)-tb)*3+2)

        int *rn=r+n,*san=sa+n,ta=0,tb=(n+1)/3,tbc=0,p;
        r[n]=r[n+1]=0;
        for(int i=0;i<n;i++) if(i%3!=0) wa[tbc++]=i;

        sort(r+2,wa,wb,tbc,m);
        sort(r+1,wb,wa,tbc,m);
        sort(r,wa,wb,tbc,m);

        rn[F(wb[0])]=0; p=1;
        for(int i=1;i<tbc;i++) rn[F(wb[i])]=c0(r,wb[i-1],wb[i])?p-1:p++;

        if(p<tbc) dc3(rn,san,tbc,p);
        else for(int i=0;i<tbc;i++) san[rn[i]]=i;   
        /**以上是第一部分計算完畢*/

        for(int i=0;i<tbc;i++) if(san[i]<tb) wb[ta++]=san[i]*3;
        if(n%3==1) wb[ta++]=n-1;
        sort(r,wb,wa,ta,m);  
        /**以上是第二部分計算完畢*/

        /**合並*/
        for(int i=0;i<tbc;i++) wv[wb[i]=G(san[i])]=i;
        int i=0,j=0;
        for(p=0;i<ta&&j<tbc;p++)
        {
            sa[p]=c12(wb[j]%3,r,wa[i],wb[j])?wa[i++]:wb[j++];
        }
        while(i<ta) sa[p++]=wa[i++];
        while(j<tbc) sa[p++]=wb[j++];

        #undef F(x)
        #undef G(x)
    }

public:
    /**字符串S,長度n,S[0,n-1]，為方便，我們假設它只包含小寫字母
       最后的后綴數組存儲在SA[1~n]中  0<=SA[i]<=n-1
    **/
    void calSuffixArray(char *S,int n,int *SA)
    {
        for(int i=0;i<n;i++) r[i]=S[i]-'a'+1;
        r[n+1]=0;
        dc3(r,sa,n+1,27);
        for(int i=1;i<=n;i++) SA[i]=sa[i];
    }
};

到此位置，這就是求后綴數組sa的兩種方式。 下面，來介紹一個新的數組height。

在求出sa數組之后，我們定義一個新的數組height,height[i]表示Suffix[sa[i-1]]與Suffix[sa[i]]的最長公共前綴，也就是排名為i和排名為i-1的兩個后綴的最長公共前綴。如果我們求出了height數組，那么對於任意兩個位置i,j,我們不妨設rank[i]小於rank[j]，它們的最長公共前綴就是height[rank[i]+1]，height[rank[i]+2],...,height[rank[j]]的最小值。比如字符串為aabaaaab，我們求后綴Suffix[1]=abaaaab和Suffix[4]=aaab的最長公共前綴，如下圖所示

 

那么如何計算height？我們定義h[i]=height[rank[i]]，也就是Suffix[i]和它前一名的最長公共前綴，那么很明顯有h[i]>=h[i-1]-1。因為h[i-1]是Suffix[i-1]和它前一名的最長公共前綴，設為Suffix[k]，那么Suffix[i]和Suffix[k+1] 的最長公共前綴為h[i-1]-1，所以h[i]至少是h[i-1]-1。所以我們可以按照求h[1],h[2],h[3] 順序計算所有的height。代碼如下

const int N=100005;
int Rank[N];

void calHeight(int *r,int *sa,int n,int *height)
{
    int i,j,k=0;
    for(int i=1;i<=n;i++) Rank[sa[i]]=i;
    for(int i=0;i<n;i++)
    {
        if(k) k--;
        j=sa[Rank[i]-1];
        while(i+k<n&&j+k<n&&r[i+k]==r[j+k]) k++;
        height[Rank[i]]=k;
    }
}