后綴數組初探

后綴數組

​ 本文總結了后綴數組(Suffix Array,SA)的倍增算法以及如何在O(n)預處理、O(1)查詢的時間復雜度內求得任意兩個后綴的最長公共前綴(Longest Common Prefix,LCP)。

1 基本定義

1. 后綴i (suffix[i])：從下標i起始的后綴。（特別地，認為字符串本身也是自己的后綴）
2. 后綴數組 (Saffix Array,SA)：將后綴0\(\rightarrow\)N-1按字典序從小到大排列，SA[i]為第i (0\(\rightarrow\)N-1)小后綴的起始位置。
3. 名次數組 (Rank)：將后綴0\(\rightarrow\)N-1按字典序從小到大排列，Rank[ i (0\(\rightarrow\)N-1)]為后綴i的名次。
4. 高度數組 (Height)：Height[i (0\(\rightarrow\)N-1)]為suffix[ SA[i] ]和suffix[ SA[i-1] ]的最長公共前綴。（Height[0]沒有意義）
5. 輔助數組(H)：H[ i (0\(\rightarrow\)N-1) ]為Height[ Rank[i] ]。

2 基本性質

1. 后綴數組與名次數組互逆：SA[ rank[i] ]=i, Rank[ SA[i] ]=i。

2. 后綴i,j的LCP為min{ Height[ Rank[i]+1\(\rightarrow\)Rank[k] ] }。

   性質過於顯然，證明略。

3. 輔助數組中H[i]\(\ge\)H[i-1]-1。

   將所有后綴排序后，假設排在suffix[i-1]的前一個是suffix[k]，將兩個后綴分別刪除首字符，可得到suffix[i]和suffix[k+1]，結合（2）有 : (忽略suffix[i-1]本身就是排在第一個的后綴和suffix[k]長度為1的情形)

   1. suffix[k+1]必然排在suffix[i]前面;

   2. LCP ( suffix[i],suffix[k+1] ) = LCP ( suffix[i-1],suffix[k])-1=Height[ Rank[i-1] ]-1=H[i-1]-1;

   3. LCP ( suffix[i],suffix[k+1] ) = Min Height[ Rank[k+1]+1\(\rightarrow\)Rank[i] ];

   4. H[i]=Height[ Rank[i] ]\(\ge\)Min Height[ Rank[k+1]+1\(\rightarrow\)Rank[i] ];

   綜合以上，證畢。

3 后綴數組的倍增算法

​ 首先算出每個字母的Rank，然后利用Rank給所有后綴的前兩個字符（不存在的字符認為它是無窮小）排序得到以每個二元組（字符+字符）的Rank，如此再給所有后綴的前四個字符排序得到以每個二元組（2*字符+2*字符）的Rank……迭代至每個二元組的Rank各不相同，這就是SA的倍增算法。

void suffixArray() {
    for (int i=0; i<n; i++) c[s[i] ]++;
    for (int i=1; i<128; i++) c[i]+=c[i-1];
    for (int i=n-1; ~i; i--) rank[i]=c[s[i] ]--;
    for (int k=1,p=0; p!=n && k<=n; k<<=1) {
        for (int i=0; i<n; i++) b[i]=make_pair(make_pair(rank[i],rank[i+k]),i);
        sort(b,b+n), p=0; //利用sort排序二元組
        for (int i=0; i<n; i++) {
            if (i && b[i].first == b[i-1].first) rank[b[i].second]=p; //計算每個位置的rank
            else rank[b[i].second]=++p;
        }
    }
    for (int i=0; i<n; i++) SA[rank[i]-1]=i+1;
    for (int i=0; i<n; i++) printf("%d ",SA[i]);
    return 0;
}

​ 容易看出，利用快排排序二元組的倍增算法為O(NlogNlogN)。

3.1 倍增算法的基數排序優化

​ 注意到每輪對二元組的排序中，第二關鍵字的排名可以直接由上一次排序的得到的Rank推出，利用基數排序里LSD的做法，第二關鍵字求得名次后，直接對第一關鍵字開（穩定的）桶排序即可。正是同函數開頭對單個字符求名次一樣的做法。

void suffixArray(char*s,int*x,int*y,int*sa) {
    int i,k,p,n=strlen(s),m=128;
    for(i=0; i<n; ++i) ++c[x[i]=s[i] ];
    for(i=1; i<m; ++i) c[i]+=c[i-1];
    for(i=n-1; ~i; --i) sa[--c[x[i] ]]=i;
    for(k=1; k<=n; k<<=1) {
        for(i=n-k,p=0; i<n; ++i) y[p++]=i;
        for(i=0; i<n; ++i) if(sa[i]>=k) y[p++]=sa[i]-k;//y[i]第二關鍵字排名i的第一關鍵字位置
        for(i=0; i<m; ++i) c[i]=0;
        for(i=0; i<n; ++i) ++c[x[y[i] ]];
        for(i=1; i<m; ++i) c[i]+=c[i-1];
        for(i=n-1; ~i; --i) sa[--c[x[y[i] ]] ]=y[i];//基排求得二元組名次為[?]的第一關鍵字位置
        swap(x,y), p=1, x[sa[0] ]=0; //y上次排序后各后綴的前綴的名次；x本次排序后后綴的前綴的名次
        for(i=1; i<n; ++i) x[sa[i] ]=//計算本次排序后二元組的名次
            (y[sa[i] ]==y[sa[i-1] ]&&y[sa[i]+k]==y[sa[i-1]+k])?p-1:p++;
        if((m=p)>=n) break;
    } 
}

​ 顯然，這樣優化后復雜度降為O(NlogN)

4 利用輔助數組求高度

​ 由基本性質3，利用輔助數組H計算可以減少字符比較次數，實現O(n)的做法（暴力時間復雜度O(n^2^)）。注意，代碼實現中並不需要開一個真正的H數組。

void heightArray() {
    int i,j,k=0;
    for(i=0; i<n; ++i) rank[sa[i] ]=i;
    for(i=0; i<n; ++i) {
        if(k) --k;
        if(rank[i]) p=sa[rank[i]-1];
        else {height[0]=0; continue;} //已改正原書上的數組越界的錯誤
        while(s[i+k]==s[j+k]) ++k;
        height[rank[i] ]=k;
	}
}

參考材料：《算法競賽入門經典——訓練指南》，劉汝佳、陳鋒著，清華大學出版社