转自:
后缀数组(Suffix Array)——理论和思想
这一篇和下一篇博文我准备写一下我在参加ACM/ICPC期间曾经研究过的后缀数组。关于后缀数组,网上有很多英文资料,但是很多现在的研究结果都是受1991年Udi Manber & Gene Myers的《Suffix arrays: a new method for on-line string searches》中所提的方法的启发,采用倍增的思想。当然,现在有国外学者提到的三分+分治的线性构造方法除外。
写这两篇文章主要是为了把一种思想写下来,同时练练已经生疏了很久的算法。我曾经在老的博客上写过一个比较丑陋的后缀数组构造算法,在后一篇文章中我将结合近期我看到的资料对它进行优化,使其变得比较美观 :-)
我们定义一个字符串A,其表示:
A = a0a1a2…an-1
其中ai(0 ≤ i < n)都是字符集E中的字符,这个字符集是全序的,也就是说其中任意两个字符都是可以比较大小的。例如,用ASCII码编码的英文字符集就是这样的一个E的特例。
|A|表示字符串的长度,其值为n;Ai(0 ≤ i < n)表示一个后缀,它其实表示一个字符串:
Ai = aiai+1ai+2…an-1
我们让运算符“≤”表示两个串按照字典序比较,然后定义运算符“≤h”表示两个串的前h个字符按照字典序比较(=h、<h等同理),那么就有:
结论1 若Aj =h Ak且Aj+h ≤h Ak+h,则Aj ≤2h Ak (j+h, k+h < n,“≤”换成“=、<、>”等等依然成立)
这个是1989年Udi Manber & Gene Myers发明 nlogn 后缀数组生成算法的主要理论依据。然而他们的天才之处不是在于看到了这个结论而是将这个结论与“复用”的思想结合在了一起。
什么是后缀数组呢?
后缀数组(Suffix Array)是一个存放索引的数组,如果把这个数组命名为SA,那么有:
A[SA[x]] ≤ A[SA[y]],其中0 ≤ x ≤ y < n
要产生这样一个数组,我们可以用最普通的sort/qsort结合strcmp,这个看起来是一个不错的想法,但是考虑到strcmp其实不是常数时间复杂度而是线性时间复杂度的,所以这个算法就显得不是那么高效了。Udi Manber & Gene Myers的聪明之处就是将结论1中的h取成了“1、2、4、8、16…”这样的指数数列,只要每趟比较保证字符串A的所有后缀按≤h有序,那么相应地,每个后缀扩展成2h长度之后,只要比较其后h个字符就行了,而这后h个字符其实就是其他某个后缀的前h个字符,其实已经比较过了,直接结合结论1可以得出结果,仅仅花了常数时间。这样,经过logn趟比较,后缀数组就生成了。例如:
A = abba
A0 = abba
A1 = bba
A2 = ba
A3 = a
h = 1 时,A0 A1 A2 A3 的≤h有序序列为:A0 =h A3 <h A2 =h A1
h = 2 时,要决定A1与A2的≤2h比较结果,因为 A2 =h A1,根据结论1,只要看a与ba的≤h比较结果就行了,而我们欣喜地发现,这个结果其实就是A2与A3的≤h比较结果,在h = 1的时候早就得出了结论——A3 <h A2!
所以,我们只要经过 logn 趟比较,每趟比较花费O(n)的时间进行2个关键字的桶排序,那么就可以得到一个后缀数组了!
PS:现在有更好的线性的结果了,但是算法相对复杂,我也没有怎么看过 :-)
-----------------------------------------------------------------------------------
转自:
后缀数组(Suffix Array)——实现与应用
我在前面一篇文章中已经概要地讲了后缀数组的基本理论依据,下面结合一个 ACM/ICPC 竞赛题目来说说后缀数组的简单应用。我们先来实现后缀数组 O(nlogn) 的构造算法。我曾经在老的博客上写过一个比较丑陋的后缀数组构造算法,我在产生写这两篇文章的想法时,有去网上搜了一下,看了别人的一些实现和一些以前留下的论文,现对之前的算法进行优化,使其变得比较美观一些 :-)
我的构造算法用了O(4n)的空间复杂度,这个和Udi Manber & Gene Myers的论文中提到的O(2n)的空间复杂度还是有差距的,但是考虑如果按照他们的算法写出来,那么代码必然更长更臭(我之前那个算法就是受了他们思想的很大影响才造就了其丑陋程度),所以还是牺牲一点空间吧。此外,我还看到过几个空间为 O(2n) 而且比一般O(nlogn) 快的算法,但是代码和思想都非常复杂,不利于掌握。
定义一种类型:
1 typedef unsigned char uchar;
后缀数组构造算法:
01 void CreateSuffixArray(uchar* szText,
02 int L, int** _S, int** _R, int** _T1, int** _T2)
03 {
04 int i, h, h2, *T, *S1, *S2, *R, *B;
05
06 S1 = *_S; // h阶后缀数组
07 S2 = *_T1; // 2h阶后缀数组
08 R = *_R; // h阶Rank数组
09 B = *_T2; // 某个桶空余空间尾部的索引,兼任2h阶Rank数组
10
11 // 花O(n)的时间对h = 1进行计数排序
12 for(i = 0; i < 256; i++)
13 B = 0;
14 for(i = 0; i < L; i++)
15 B[szText]++;
16 for(i = 1; i < 256; i++)
17 B += B[i - 1];
18 for(i = 0; i < L; i++)
19 S1[--B[szText]] = i;
20
21 // 计算Rank(1),因为仅仅是1阶的Rank,所有有并列的
22 for(R[S1[0]] = 0, i = 1; i < L; i++)
23 {
24 if(szText[S1] == szText[S1[i - 1]])
25 R[S1] = R[S1[i - 1]];
26 else
27 R[S1] = R[S1[i - 1]] + 1;
28 }
29
30 // log(n)趟O(n)的倍增排序
31 // SA(h) => Rank(h) => SA(2h) => Rank(2h) => …
32
33 for(h = 1; h < L && R[S1[L - 1]] < L - 1; h <<= 1)
34 {
35 // 计算Rank(h)相同的后缀形成的h桶尾部的索引
36 // 即有多少个后缀的h前缀相同,它们被放在一个桶中
37 for(i = 0; i < L; i++)
38 B[R[S1]] = i;
39
40 // 求SA(2h)
41 // 在同一个h桶中,所有的后缀的h前缀肯定相同,
42 // 那么比较他们的2h前缀,只要比较其2h前缀后半的
43 // 长度为h的串即可,而这个串恰恰是后面某个后缀的
44 // h前缀,所以我们逆向遍历有序的SA(h),
45 // 将S1 – h号前缀放到它所在桶的{zh1}一个空位置,
46 // 同时,桶尾前进一个位置,这样即形成了2h桶排序
47 for(i = L – 1; i >= 0; i–)
48 if(h <= S1)
49 S2[B[R[S1 - h]]–] = S1 – h;
50
51 // 对于长度不超过h的{zh1}几个后缀,由于在h阶段
52 // 它们每个实际上都已经独立分桶了(长度为h的也是)
53 // 而且他们的桶中有且仅有一个元素,
54 // 所以只要直接复制他们h阶段的SA值就可以了
55 // 同时,由于采用滚动数组,所以S2中“残留”了
56 // h/2个有效的数据,所以最终我们只需复制h/2个数据
57 for(i = L – h, h2 = L – (h >> 1); i < h2; i++)
58 S2[B[R]] = i;
59
60 T = S1; S1 = S2; S2 = T;
61
62 // 计算Rank(2h)
63 // 2h阶段是否要分桶只需看相邻两个2h前缀前后两半
64 // h前缀是否全部h阶相等
65 for(B[S1[0]] = 0, i = 1; i < L; i++)
66 {
67 // 这里不用考虑S1 + h会越界
68 // 如果i达到了S1 + h越界的数值,
69 // 那么前面一个条件显然不会满足了
70 // 因为此时i前缀肯定已经独立分桶了
71 if(R[S1] != R[S1[i - 1]] ||
72 R[S1 + h] != R[S1[i - 1] + h])
73 {
74 B[S1] = B[S1[i - 1]] + 1;
75 }
76 else
77 B[S1] = B[S1[i - 1]];
78 }
79
80 T = B; B = R; R = T;
81 }
82
83 if(*_S != S1)
84 *_S = S1, *_T1 = S2;
85 if(*_R != R)
86 *_R = R, *_T2 = B;
87 }
介绍一个重要概念:LCP!LCP是Longest Common Prefix的缩写,即最长公共前缀,表示某个串从{dy}个字符开始对应位置字符相同的连续的位置数。比如,后缀abcda和后缀abcca的LCP就是3。我们将后缀数组中连续的两个后缀Ai-1和Ai的LCP称为Ai的Height,即Height(i) = LCP(j , j – 1),并规定ASA[0]的Height为0。那么很显然,后缀数组某个区间的两个区间边界元素所表示的后缀的LCP就是区间内所有元素所代表的后缀的Height的最小值。我们要求这个LCP,就相当于一个RMQ(Range Minimum Query)问题,当Height已知的时候,只要常数时间就可以求出RMQ,即所求的LCP。所以,关键是如何降低求Height数组的复杂度。不过人们发现Height数组有一个令人兴奋的性质。令 h(x) = Height(Rank(x)),即x号前缀的Height值,那么,
当 x > 0 且 Rank(x) > 0 时, h(i) ≥ h(i – 1) – 1
这个在这里就不证明了,反正证明过程相当巧妙 :-) 利用这个性质,有了下面的这个线性的求Height的算法:
01 void CalculateHeight(uchar* szText,
02 int L, int* S, int* R, int* H, int* T)
03 {
04 int i, j, k;
05
06 for(k = 0, i = 0; i < L; i++)
07 {
08 if(R == 0)
09 H = 0;
10 else
11 {
12 for(j = S[R - 1]; szText[i + k] == szText[j + k]; k++);
13
14 H[R] = k;
15
16 if(k > 0)
17 k–;
18 }
19 }
20 }
初一看,这个不是 O(n2) 的吗?其实根据上面说的性质,可以证明,它是线性的,证明也略了
下面是一个具体的ACM/ICPC竞赛题目的解法,原题你可以在这里找到:
01 char C[200002];
02 int D[4][200001];
03
04 int main()
05 {
06 int i, l1, l2, b;
07 int *S, *R, *H, *T;
08
09 gets(C);
10 l1 = (int)strlen(C);
11 C[l1] = '$';
12 gets((char*)C + l1 + 1);
13 l2 = l1 + 1 + (int)strlen(C + l1 + 1);
14
15 S = D[0]; R = D[1];
16 H = D[2]; T = D[3];
17
18 CreateSuffixArray((uchar*)C, l2, &S, &R, &H, &T);
19 CalculateHeight((uchar*)C, l2, S, R, H, T);
20
21 // 求两个串的最长公共子串,只要让两个串s1、s2
22 // 连接在一起形成一个新串,求出新串的SA、Rank和Height
23 // 很显然,最长公共子串肯定出现在后缀数组某相邻两项之中
24 // 根据Height的定义,扫描一遍Height数组,找相邻两个分别开始于
25 // s1和s2串某个位置的后缀,求出所有满足这个条件的{zd0}Height即可
26
27 for(b = 0, i = 1; i < l2; i++)
28 {
29 if(S < l1 && S[i - 1] > l1 ||
30 S > l1 && S[i - 1] < l1)
31 {
32 if(H > b)
33 b = H;
34 }
35 }
36
37 printf(“%d\n”, b);
38
39 return 0;
40 }
后缀数组的用处很大,除了上面的求两个串的最长公共字串之串之外,多模式匹配、最长回文串、全文检索等等都它的拿手好戏,可以说后缀数组是后缀树良好的替代品。
相关链接: