后缀数组

摘要

后缀数组由 Manber & Myers 在 1990 年首先提出《Suffix arrays: a new method for on-line string searches》,用以替代后缀树,并且改进了存储空间的需求。

简介

后缀数组(Suffix Array)是某一字符串的所有后缀按照字典序的一个排列。

数组的索引从 0 开始。

后缀数组

考虑以下字符串:

i    : 1 2 3 4 5 6
s[i] : b a n a n a

对应的后缀数组如下:

suf[1]:banana
suf[2]:anana
suf[3]:nana
suf[4]:ana
suf[5]:na
suf[6]:a

按字典序排序

suf[6]:a
suf[4]:ana
suf[2]:anana
suf[1]:banana
suf[5]:na
suf[3]:nana

则后缀数组为

i   : 1 2 3 4 5 6
sa[i]:6 4 2 1 5 3

表示字典序排第 i 个的后缀的首字母下标

则每一个后缀对应排序后的位置为:

i    :1 2 3 4 5 6
rk[i]:4 3 6 2 5 1

表示下标为 i 的后缀的排名。

构造后缀数组

后缀数组的关键在于sa[]rk[].

直接使用字符串排序,复杂度 .

使用桶排序 + 倍增可将复杂度降低到 .

原理

如果我们已知每一个后缀,长度为 的前缀的sa[]rk[],那么我们可以在 O(n) 时间内推导出长度为 2j 的sa[]rk[]

每次 倍增,即可在 O(nlogn) 时间内求出后缀数组。

算法剖析

考虑 s[]=”banana”.
SuffixArray.png

对于长度为 1 的子串,直接预处理rk[](有相同的存在)

当 j=1 时,我们已知长度为 1 的子串的rk[];将两个长度为 1 的子串拼在一起,他们的rk[]组合成一个二元组(最后一位后面没有字符,默认为 0)。

怎么推导出长度为 2 的子串的后缀数组呢?

对二元组桶排序。第一关键字为第一元,第二关键字为第二元。我们将它理解为,先按前 个排序;如果相等,再按后 个排序。

最后更新rk[]s[],j 倍增。

循环上述算法,直到 大于了序列长度,则算法结束。

背板要素

后缀数组的板子,笔者认为十分难背。这里笔者整理一下后缀数组模板的要素帮助背板

变量: ,数组 .

桶排序:排序对象是 ,排序关键字是 。(相对保持不变的关键字为 ).

初始化: .

第一次:当做步长为0,那么一二关键字重合,所以 则由字符决定, 求得

倍增:

  1. ,按第二关键字排列 在过程中可能重复,但 中的元素不会重复);
  2. ,按第一关键字给 排序(桶排序);
  3. ,用旧的 配合 求新的 (交换的技巧)

代码实现

#include<bits/stdc++.h>
using namespace std;
const int N=1e6+6;

namespace SA{
int l;// 字符串长度,1 起点
int sa[N],rk[N];
int t[N],bin[N],sz;//sz: 桶的大小
void qsort(){
for(int i=0;i<=sz;i++)bin[i]=0;
for(int i=1;i<=l;i++)bin[rk[t[i]]]++;// 可以用 bin[rk[i]],不过用 bin[rk[t[i]] 可读性更强
for(int i=1;i<=sz;i++)bin[i]+=bin[i-1];
for(int i=l;i>=1;i--)sa[bin[rk[t[i]]]--]=t[i];
}
void make(char *s){
l=strlen(s+1),sz=max(l,'z'-'0'+1);// 初始化 l,sz
for(int i=1;i<=l;i++)t[i]=i,rk[i]=s[i]-'0'+1;// 初始化 t,rk(1 为最高排名)
qsort();// 对 t 排序,结果转移到sa中
for(int j=1;j<=l;j<<=1){// 开始倍增
int tot=0;
for(int i=l-j+1;i<=l;i++)t[++tot]=i;// 这些后缀没有第二关键字,因此排在最前面
for(int i=1;i<=l;i++)if(sa[i]-j>0)t[++tot]=sa[i]-j;//{sa[i]-j}是按照第二关键字排列的
qsort();
swap(t,rk);//t 变成上一次的 rk
rk[sa[1]]=tot=1;
for(int i=2;i<=l;i++)
rk[sa[i]]=t[sa[i-1]]==t[sa[i]]&&t[sa[i-1]+j]==t[sa[i]+j]?tot:++tot;
}
}
}
char s[N];
int main(){
scanf("%s",s+1);
SA::make(s);
for(int i=1;i<=SA::l;i++)printf("%d ",SA::sa[i]);
return 0;
}

首先要知道, 表示后缀的下标,并且取遍 中的元素。{tmp[i]}数列满足后缀序列 是按字典序排列的,其中当 即为空串

也就是说, 是已经以第二关键字排好序的。而我们要完成的工作是这样的:把 以第一关键字(第一关键字就是前 个字符的字典序)排序,并保留第二关键字(第二关键字就是第 的字符串的字典序)的相对顺序不变。把排序的结果存储在

而我们的桶排序是一种稳定排序,因此使用桶排序时,同一个桶内的元素的相对顺序不会被打乱,也就保证了第二关键字的顺序

再看看我们的桶排序函数在干什么。把相同 的放在一个桶里面,这里的 ,不是 ,即我们以第一关键字为桶。我们没有模拟出二元组,那么我们是怎么实现二元组排序的呢?那就是,把第一关键字放在桶里面,然后按第二关键字的顺序拿出来!

什么叫按第二关键字的顺序?观察第三个循环的赋值,我们倒序枚举 i,即倒序遍历 ,然后减掉当前的桶的前缀和,这个前缀和就是当前 在排序后的的新的位置。也就是说,我们每一个桶内的元素也是倒着取出来的

别忘了桶内的元素是第一关键字,既然同一个桶内的第一关键字是一样的。桶的位置由其前缀和的差分决定(下标区间

同一个桶内的元素的顺序由第二关键字决定而我们是按第二关键字倒序的顺序来枚举对应的后缀下标,用这个下标的 找到其对应的第一关键字(即 ),也就是桶,然后在这个桶里面拿走末尾的下标。即我们倒着枚举倒着拿走,这等价于正着枚举正着拿走,那么就能够维护第二关键字的相对顺序

综上,采用这样的排序策略是正确的

排完序后就要统计新的 对应的 啦。怎么统计?这里当然是直接二元组比较啦。sa[i],sa[i-1] 是字典序相邻的两个后缀的下标,其对应的 就是上一次的 做了交换)。而 则是这两个下标后第 个下标,即第二关键字,对应的 就是第二关键字的

那么这就是 二元组啦,直接比较是否相同就好啦。另外, 不可能变小,只可能从相等变得更大(因为第一关键字排了序)

模板题 LuoguP3809

Height&h 数组

利用后缀数组维护信息,需要一些辅助数组。

定义h[i]表示suf[i]suf[sa[rk[i]-1]]的最长公共前缀(LCP),也就是字典序排序后比 suf[i] 字典序小 1 的后缀与 suf[i] 的 LCP

定义height[i]=h[sa[i]],即字典序排名第 i 个的后缀与他前一个后缀的 LCP。

构造 h 数组

结论一则:h[i]>=h[i-1]-1 (0<i<n). 证明很显然,直接用 h[i-1] 的⽅案去掉头上的元素即可

height[] 可直接索引 h[] 得到,而构造 h[]:

int calc(int x,int y,int len){
while(x+len<=l&&y+len<=l&&s[x+len]==s[y+len])++len;
return len;
}
void calc_h(){
for(int i=1;i<=l;i++)h[i]=(rk[i]==1)?0:calc(i,sa[rk[i]-1],max(h[i-1]-1,0));
}

均摊复杂度O(n).

h 数组的模板见 SP609 和 SP705(双倍经验)

一个完整的板子

#include<bits/stdc++.h>
using namespace std;
const int N=3e4+5;

struct SA{
char s[N];
int l;
int sa[N],rk[N];
int t[N],bin[N],sz;
int h[N],he[N];//h,height
void qsort(){
for(int i=0;i<=sz;i++)bin[i]=0;
for(int i=1;i<=l;i++)bin[rk[t[i]]]++;
for(int i=1;i<=sz;i++)bin[i]+=bin[i-1];
for(int i=l;i>=1;i--)sa[bin[rk[t[i]]]--]=t[i];
}
void make(){
l=strlen(s+1),sz=max(l,26);
for(int i=1;i<=l;i++)t[i]=i,rk[i]=s[i]-'a'+1;
qsort();
for(int j=1;j<=l;j<<=1){
int tot=0;
for(int i=l-j+1;i<=l;i++)t[++tot]=i;
for(int i=1;i<=l;i++)if(sa[i]-j>0)t[++tot]=sa[i]-j;
qsort();
swap(t,rk);
rk[sa[1]]=tot=1;
for(int i=2;i<=l;i++)
rk[sa[i]]=t[sa[i-1]]==t[sa[i]]&&t[sa[i-1]+j]==t[sa[i]+j]?tot:++tot;
}
}
int move(int x,int y,int len){
while(x+len<=l&&y+len<=l&&s[x+len]==s[y+len])++len;
return len;
}
void calc_h(){
for(int i=1;i<=l;i++)
h[i]=rk[i]==1?0:move(i,sa[rk[i]-1],max(h[i-1]-1,0));
}
int st[N][16];//h[sa[i]]~h[sa[i+2^j]]中的最小值
void make_st(){
for(int i=1;i<=l;i++)st[i][0]=h[sa[i]],printf("st[%d,0]:%d\n",i,st[i][0]);
for(int j=1;(1<<j)<=l;j++){
int step=1<<(j-1);
for(int i=1;i+step<=l;i++){
st[i][j]=min(st[i][j-1],st[i+step][j-1]);
printf("st[%d,%d]:%d\n",i,j,st[i][j]);
}
}
}
int lcp(int x,int y){//返回长度
if(x==y)return l-x+1;
x=rk[x],y=rk[y];
if(x>y)swap(x,y);
x++;//取不到x
int step=log(y-x+1)/log(2);
return min(st[x][step],st[y-(1<<step)+1][step]);
}
};
SA sa;
int main(){
scanf("%s",sa.s+1);
sa.make();
sa.calc_h();
sa.make_st();
for(int i=1;i<=sa.l;i++){
for(int j=1;j<=sa.l;j++){
if(i!=j&&sa.lcp(i,j))printf("%d %d %d\n",i,j,sa.lcp(i,j));
}
}
return 0;
}
/*
* BUG#1: bin[i-1]写成bin[i+1]
* 求后缀(i,j)(rk[i]<rk[j])的LCP,就是求 j in (rk[i],rk[j]] 的height[j] 的最小值,即RMQ.
*/