后缀数组初探

后缀数组

本文总结了后缀数组(Suffix Array,SA)的倍增算法以及如何在O(n)预处理、O(1)查询的时间复杂度内求得任意两个后缀的最长公共前缀(Longest Common Prefix,LCP)。

1 基本定义

后缀i (suffix[i])：从下标i起始的后缀。（特别地，认为字符串本身也是自己的后缀）
后缀数组 (Saffix Array,SA)：将后缀0( ightarrow)N-1按字典序从小到大排列，SA[i]为第i (0( ightarrow)N-1)小后缀的起始位置。
名次数组 (Rank)：将后缀0( ightarrow)N-1按字典序从小到大排列，Rank[ i (0( ightarrow)N-1)]为后缀i的名次。
高度数组 (Height)：Height[i (0( ightarrow)N-1)]为suffix[ SA[i] ]和suffix[ SA[i-1] ]的最长公共前缀。（Height[0]没有意义）
辅助数组(H)：H[ i (0( ightarrow)N-1) ]为Height[ Rank[i] ]。

2 基本性质

后缀数组与名次数组互逆：SA[ rank[i] ]=i, Rank[ SA[i] ]=i。
后缀i,j的LCP为min{ Height[ Rank[i]+1( ightarrow)Rank[k] ] }。

性质过于显然，证明略。
辅助数组中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( ightarrow)Rank[i] ];
4. H[i]=Height[ Rank[i] ](ge)Min Height[ Rank[k+1]+1( ightarrow)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²)）。注意，代码实现中并不需要开一个真正的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;
	}
}

参考材料：《算法竞赛入门经典——训练指南》，刘汝佳、陈锋著，清华大学出版社