后缀数组简介
在计算机科学里, 后缀数组(英语:suffix array)是一个通过对字符串的所有后缀经过排序后得到的数组。此数据结构被运用于全文索引、数据压缩算法、以及生物信息学。
后缀数组被乌迪·曼伯尔与尤金·迈尔斯于1990年提出,作为对后缀树的一种替代,更简单以及节省空间。它们也被Gaston Gonnet 于1987年独立发现,并命名为“PAT数组”。
后缀数组(Suffix Array)本质上是对一个字符串的所有后缀进行排序,之后会得到两个信息,分别是:
- S A [ i ] SA[i] SA[i],表示排名为 i i i 的后缀的起始位置。
- R a n k [ i ] Rank[i] Rank[i],表示起始位置为 i i i 的后缀的排名。
S A SA SA数组和 R a n k Rank Rank 数组是类似反函数的关系,有: S A [ R a n k [ i ] ] = R a n k [ S A [ i ] ] = i SA[Rank[i]]=Rank[SA[i]]=i SA[Rank[i]]=Rank[SA[i]]=i,也就是说如果我求出了其中一个数组,就可以很快知道另一个数组。利用这两个信息,可以处理很多字符串相关的问题。
SA的求法
显然有 O ( n 2 log n ) O(n^2\log{n}) O(n2logn) 的暴力做法:即对每个后缀进行快速排序,由于字符串比较大小的复杂度是 O ( l e n ) O(len) O(len),因此总的复杂度是 O ( n 2 log n ) O(n^2\log{n}) O(n2logn),我们需要一个办法来优化他。
这里利用了倍增的思想(类似 S T ST ST表)。我们给 S A SA SA 数组和 R a n k Rank Rank 数组增加一维:令 S A [ i ] [ k ] SA[i][k] SA[i][k]表示,在所有长度为 2 k 2^k 2k的子串中,排名第 i i i 的子串的起始位置。类似地定义 R a n k [ i ] [ k ] Rank[i][k] Rank[i][k] 表示 起始位置为 i i i,长度为 2 k 2^k 2k 的子串,在所有长度为 2 k 2^k 2k的子串中的排名。如果子串长度不足 2 k 2^k 2k 则用 ‘\0’(字典序最小字符)补全。
所以现在可以写出算法大概的流程,以求 R a n k Rank Rank 数组为例。
- 先求出 R a n k [ 1... l e n ] [ 0 ] Rank[1...len][0] Rank[1...len][0]
- 再求出
R
a
n
k
[
1...
l
e
n
]
[
1
]
Rank[1...len][1]
Rank[1...len][1]
⋮ \vdots ⋮ - 求出 R a n k [ 1... l e n ] [ n ] Rank[1...len][n] Rank[1...len][n]
显然当 2 n > = l e n 2^n>=len 2n>=len 时,子串排序就相当于后缀排序,因为此时所有子串都延申到了字符串的末尾。算法也就完成了。
那么如何在已知 R a n k [ i ] [ k ] Rank[i][k] Rank[i][k] 的 情况下求 R a n k [ i ] [ k + 1 ] Rank[i][k+1] Rank[i][k+1] 呢。我可以把长度为 2 k + 1 2^{k+1} 2k+1 的子串(也就是现在要求排名的)分成两个长度为 2 k 2^k 2k 的部分。这两部分的起始点分别是 i i i 和 i + 2 k i+2^k i+2k。且这两部分的排名都是已知的,分别是 R a n k [ i ] [ k ] Rank[i][k] Rank[i][k] 和 R a n k [ i + 2 k ] [ k ] Rank[i+2^k][k] Rank[i+2k][k]。因此在对长度为 2 k + 1 2^{k+1} 2k+1 的子串进行排序时,可以直接利用 R a n k [ i ] [ k ] Rank[i][k] Rank[i][k] 来比较两个子串的大小。
假设要进行比较的两个长度为 2 k + 1 2^{k+1} 2k+1 的子串的起始点分别是 i i i 和 j j j,那么我们可以先比较 R a n k [ i ] [ k ] Rank[i][k] Rank[i][k] 和 R a n k [ j ] [ k ] Rank[j][k] Rank[j][k] ,如果不相等,则比较结束。否则比较 R a n k [ i + 2 k ] [ k ] Rank[i+2^k][k] Rank[i+2k][k] 和 R a n k [ j + 2 k ] [ k ] Rank[j+2^k][k] Rank[j+2k][k]。所以这个过程相当于每一轮迭代中,对 ( R a n k [ i ] [ k ] , R a n k [ i + 2 k ] [ k ] ) (Rank[i][k],Rank[i+2^k][k]) (Rank[i][k],Rank[i+2k][k]) 这个二元组进行排序。
我们可以直接建立 P a i r Pair Pair 进行快速排序,但这样做的时间复杂度为 O ( n log 2 n ) O(n\log^2n) O(nlog2n),还不是最优解。
注意到需要排序的元素为二元组,且每个元素的值域为不超过 N N N 的整数,因此可以考虑使用基数排序来将排序过程优化到 O ( n ) O(n) O(n)。
实际实现细节与代码
在求 R a n k Rank Rank 数组的实际实现中,并不需要真的去找所有长度为 2 k 2^k 2k 的子串,也不需要多开一维数组。而是对排位信息的合并。
初始时, k = 0 k=0 k=0,也就是说我们相当于在对单个字符进行排序,这时候可以直接用字符的 ASCII 码值作为关键字。
然后 k = 1 k=1 k=1,也就是对每两个相邻字符的子串进行排序,这时候可以直接用上一次排序得到的关键字,以前一个字符的排名作为第一关键字,第二个字符的排名作为第二关键字。得到长度为 2 的子串的排序。相当于合并了上一次排序的信息。排序后重新标号。
然后 k = 2 k=2 k=2,此时可以想到,对于目前的一个位置 i i i,它的第一关键字就在位置 i i i,但他的第二关键字在位置 i + 2 i+2 i+2。因为位置 i i i 上是起始点为 i i i 的长度为 2 的子串的排名,而 i + 2 i+2 i+2 上是 i i i 为起始点的长度为 2 的子串的排名,他们两个关键字合并起来的二元组,就是 第 i i i 个长度为 4 的子串的排名。
后面以此类推。
下面以字符串 s t r = " a a b a b a b " str="aababab" str="aababab",演示一下这个过程。
- k = 0,直接按字符大小给定关键字
str | a | a | b | a | b | a | b |
---|---|---|---|---|---|---|---|
rank | 1 | 1 | 2 | 1 | 2 | 1 | 2 |
- k = 1,按照 ( i , i + 1 ) (i,i+1) (i,i+1) 的二元组给定关键字,关键字分别是:(1,1),(1,2),(2,1),(1,2),(2,1),(1,2),(2,0)。超出部分补0。重新标号结果为:
str | a | a | b | a | b | a | b |
---|---|---|---|---|---|---|---|
rank | 1 | 2 | 4 | 2 | 4 | 2 | 3 |
- k = 2,按照 ( i , i + 2 ) (i,i+2) (i,i+2) 的二元组给定关键字,关键字分别是:(1,4),(2,2),(4,4),(2,2),(4,3),(2,0),(3,0)。
str | a | a | b | a | b | a | b |
---|---|---|---|---|---|---|---|
rank | 1 | 3 | 6 | 3 | 5 | 2 | 4 |
- k = 3,按照 ( i , i + 4 ) (i,i+4) (i,i+4) 的二元组给定关键字,关键字分别是:(1,5),(3,2),(6,4),(3,0),(5,0),(2,0),(4,0)。
str | a | a | b | a | b | a | b |
---|---|---|---|---|---|---|---|
rank | 1 | 4 | 7 | 3 | 6 | 2 | 5 |
2 3 = 8 > l e n 2^3=8>len 23=8>len,算法结束,最终 R a n k = { 1 , 4 , 7 , 3 , 6 , 2 , 5 } Rank=\{1,4,7,3,6,2,5\} Rank={1,4,7,3,6,2,5},请读者自行验证。
代码
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <iostream>
#define MAXN 1000010
using namespace std;
char s[MAXN];
int N, M, SA[MAXN], Rank[MAXN<<1];
int oldrk[MAXN], id[MAXN], cnt[MAXN];
void Suffix_Array()
{
for(int i=1;i<=N;i++)
{
Rank[i] = s[i];
++cnt[Rank[i]];
}
for(int i=1;i<=M;i++)
cnt[i] += cnt[i-1];
for(int i=N;i>=1;i--)
SA[cnt[Rank[i]]--] = i;
for(int k=1;k<N;k<<=1)
{
memset(cnt, 0, sizeof(cnt));
for(int i=1;i<=N;i++)
id[i] = SA[i];
for(int i=1;i<=N;i++)
++cnt[Rank[id[i]+k]];
for(int i=1;i<=M;i++)
cnt[i] += cnt[i-1];
for(int i=N;i>=1;i--)
SA[cnt[Rank[id[i]+k]]--] = id[i];
memset(cnt, 0, sizeof(cnt));
for(int i=1;i<=N;i++)
id[i] = SA[i];
for(int i=1;i<=N;i++)
++cnt[Rank[id[i]]];
for(int i=1;i<=M;i++)
cnt[i] += cnt[i-1];
for(int i=N;i>=1;i--)
SA[cnt[Rank[id[i]]]--] = id[i];
memcpy(oldrk, Rank, sizeof(Rank));
for(int p = 0, i = 1; i <= N;i++)
{
if(oldrk[SA[i]] == oldrk[SA[i - 1]] &&
oldrk[SA[i] + k] == oldrk[SA[i - 1] + k])
{
Rank[SA[i]] = p;
}
else
{
Rank[SA[i]] = ++p;
}
}
}
}
int main()
{
scanf("%s", s+1);
N = strlen(s+1);
M = max(N, 300);
Suffix_Array();
for(int i=1;i<=N;i++)
printf("%d ", SA[i]);
return 0;
}