后缀数组
本文总结了后缀数组(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的情形)
-
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( ightarrow)Rank[i] ];
-
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(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;
}
}
参考材料:《算法竞赛入门经典——训练指南》,刘汝佳、陈锋著,清华大学出版社