后缀数组初探
后缀数组
本文总结了后缀数组(Suffix Array,SA)的倍增算法以及如何在O(n)预处理、O(1)查询的时间复杂度内求得任意两个后缀的最长公共前缀(Longest Common Prefix,LCP)。
1 基本定义
- 后缀i (suffix[i])：从下标i起始的后缀。（特别地，认为字符串本身也是自己的后缀）
- 后缀数组 (Saffix Array,SA)：将后缀0\(\rightarrow\)N-1按字典序从小到大排列，SA[i]为第i (0\(\rightarrow\)N-1)小后缀的起始位置。
- 名次数组 (Rank)：将后缀0\(\rightarrow\)N-1按字典序从小到大排列，Rank[ i (0\(\rightarrow\)N-1)]为后缀i的名次。
- 高度数组 (Height)：Height[i (0\(\rightarrow\)N-1)]为suffix[ SA[i] ]和suffix[ SA[i-1] ]的最长公共前缀。（Height[0]没有意义）
- 辅助数组(H)：H[ i (0\(\rightarrow\)N-1) ]为Height[ Rank[i] ]。
2 基本性质
-
后缀数组与名次数组互逆：SA[ rank[i] ]=i, Rank[ SA[i] ]=i。
-
后缀i,j的LCP为min{ Height[ Rank[i]+1\(\rightarrow\)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的情形)
-
suffix[k+1]必然排在suffix[i]前面;
-
LCP ( suffix[i],suffix[k+1] ) = LCP ( suffix[i-1],suffix[k])-1=Height[ Rank[i-1] ]-1=H[i-1]-1;
-
LCP ( suffix[i],suffix[k+1] ) = Min Height[ Rank[k+1]+1\(\rightarrow\)Rank[i] ];
-
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(n2)）。注意，代码实现中并不需要开一个真正的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;
}
}
参考材料：《算法竞赛入门经典——训练指南》，刘汝佳、陈锋著，清华大学出版社
