后缀数组

全文转载于 http://blog.csdn.net/xymscau/article/details/8798046

后缀数组号称字符串处理神器，不过发现好多人都只会用模板，其实这不是我们学算法的本质，我们学习算法的本质应该理解其实现原理，并加以实现，特别是算法，更讲究的是一种思想。一年前的我也是只会用别人的模板，最近却静下心来，研究了一下后缀数组，自己写了一份自己的模板。

我基本上是跟着连教的ppt来学习的，当然也少不了百度，先讲一下基本概念。(这里大量引用了连教的ppt)

基本定义：

子串  注:串!=字符串

　　字符串 S 的子串r[i..j] ， i ≤ j ，表示r 串中从 i 到 j 这一段，就是顺次排列r[i],r[i+1],...,r[j] 形成的子串。

后缀

　　后缀是指从某个位置 i 开始到整个串末尾结束的一个特殊子串。字符串r 的从 第 i 个字 符 开 始 的 后 缀 表 示 为 Suffix(i) ，也 就 是Suffix(i)=r[i..len(r)] 。

后缀数组(SA[i]存放排名第i大的子串首字符下标)

　　后缀数组 SA 是一个一维数组，它保存1..n 的某个排列 SA[1] ，SA[2] ， ……， SA[n] ，并且保证Suffix(SA[i])<Suffix(SA[i+1])， 1 ≤ i<n 。也就是将 S 的 n 个后缀从小到大进行排序之后把排好序的后缀的开头位置顺次放入SA 中。

名次数组（rank[i]存放suffix(i)的优先级）

　　名次数组 Rank[i] 保存的是 Suffix(i) 在所有后缀中从小到大排列的 “ 名次 ” 。

注：这个是排序的关键字~（这句话是我们排序的重点）

算法目标：

求得串的sa数组和rank数组

易知sa和rank互为逆操作，即sa[rank[i]] = i;

Rank[sa[i]] = i;（所以我们只要求得其一，就能O(n）算出另一个）

注：这个结论只在最后完成排序的时候符合。

但sa和rank的定义一直都是适用的。

原因是最后的时候不会存在相同（rank相等）的两个子串。

算法基本流程

•设排序的的当前长度是h。Suffix(i,h) 表示suffix(i)前h个字符(大于length会截断)

•先按H=1，对suffix(i,H)(0<i<s.length)排序

•倍增长度H,利用之前排序H/2长度后得到的rank数组作为关键字,把后H/2部分作为第二关键字,把前H/2部分作为第一关键字,对H长度的子串作排序.

•由于是倍增长度,所以最多作logn次排序

那么复杂要做到nlogn，显然排序要o(n)，O(n)一般都选计数排序。

计数排序：

http://baike.baidu.com/view/1209480.htm 不会的自己看看这个，是代码的主要部分。

这里选计数排序还有一个重要的原因，它是一个稳定排序，这就保证了数组的下标识第二关键字，我们前面说了，对于倍增长度H，利用之前排序H/2长度后得到的rank数组作为关键字,把后H/2部分作为第二关键字，嗯，就是这里，所以我们要先排后H/2的序，然后得到新的数组序列，下标就是第二关键字了，数组里面就是前H/2 rank的值，这是第一关键字，那么直接排序就相当于先对前H/2排序，如果这里相等，那么就会按下标排序，既第二关键字排序

以下内容请按代码手动模拟一个串abab的构造过程，求sa数组

具体见代码实现：

•按H=1进行计数排序 

[cpp] view plaincopy

1. //cnt是计数排序的辅助数组，k是第一关键字，id是第二关键字下标数组,r是以下标为第二关键字的新构数组，w存放的是字符串信息。sa保存的是排第i的是谁  
2. int *k = rk,*id = height,*r = res, *cnt = wa;//计数排序  
3. rep(i,up) cnt[i] = 0;  
4. rep(i,len) cnt[k[i] = w[i]]++;  
5. rep(i,up) cnt[i+1] += cnt[i];  
6. for(int i = len - 1; i >= 0; i--) {  
7.     sa[--cnt[k[i]]] = i;  
8. }  

求第二关键字（想想为什么构造w数组的时候末尾要加个0）

[cpp] view plaincopy

1. //cnt是计数排序的辅助数组，k是第一关键字，id是第二关键字下标数组,r是以下标为第二关键字的新构数组，w存放的是字符串信息，<span style="font-family: Arial, Helvetica, sans-serif;">sa保存的是排第i的是谁</span>  
2. for(int i = len - d; i < len; i++) id[p++] = i;  
3. rep(i,len) if(sa[i] >= d) id[p++] = sa[i] - d;//id保存了按后h/2排序的的序列，即排第i的后h/2的是原数组中的那个  
4. rep(i,len) r[i] = k[id[i]];//构造新的排序数组  

对新数组排序

[cpp] view plaincopy

1. //cnt是计数排序的辅助数组，k是第一关键字，id是第二关键字下标数组,r是以下标为第二关键字的新构数组，w存放的是字符串信息，sa保存的是排第i的是谁  
2. rep(i,up) cnt[i] = 0;  
3. rep(i,len) cnt[r[i]]++;  
4. rep(i,up) cnt[i+1] += cnt[i];  
5. for(int i = len - 1; i >= 0; i--) {  
6.     sa[--cnt[r[i]]] = id[i];  
7. }   

得到新的关键字(即按H长度排序后的离散序列)

[cpp] view plaincopy

1. //cnt是计数排序的辅助数组，k是第一关键字，id是第二关键字下标数组,r是以下标为第二关键字的新构数组，w存放的是字符串信息，sa保存的是排第i的是谁  
2. swap(k,r);  
3. p = 0;  
4. k[sa[0]] = p++;  
5. rep(i,len-1) {  
6.     if(sa[i]+d < len && sa[i+1]+d <len &&r[sa[i]] == r[sa[i+1]]&& r[sa[i]+d] == r[sa[i+1]+d])  
7.     k[sa[i+1]] = p - 1;  
8.     else k[sa[i+1]] = p++;  
9. }  

重复以上最后可以得到sa数组

Sa和rank有什么用？

求height数组！！

height[i] 表示sa[i]和sa[i-1]的最长前缀，height的构造看代码手推一定能弄懂，自己看看吧

下面给出全模板代码

[cpp] view plaincopy

1. #define rep(i,n) for(int i = 0;i < n; i++)  
2. using namespace std;  
3. const int size  = 200005,INF = 1<<30;  
4. int rk[size],sa[size],height[size],w[size],wa[size],res[size];  
5. void getSa (int len,int up) {  
6.     int *k = rk,*id = height,*r = res, *cnt = wa;  
7.     rep(i,up) cnt[i] = 0;  
8.     rep(i,len) cnt[k[i] = w[i]]++;  
9.     rep(i,up) cnt[i+1] += cnt[i];  
10.     for(int i = len - 1; i >= 0; i--) {  
11.         sa[--cnt[k[i]]] = i;  
12.     }  
13.     int d = 1,p = 0;  
14.     while(p < len){  
15.         for(int i = len - d; i < len; i++) id[p++] = i;  
16.         rep(i,len)  if(sa[i] >= d) id[p++] = sa[i] - d;  
17.         rep(i,len) r[i] = k[id[i]];  
18.         rep(i,up) cnt[i] = 0;  
19.         rep(i,len) cnt[r[i]]++;  
20.         rep(i,up) cnt[i+1] += cnt[i];  
21.         for(int i = len - 1; i >= 0; i--) {  
22.             sa[--cnt[r[i]]] = id[i];  
23.         }   
24.         swap(k,r);  
25.         p = 0;  
26.         k[sa[0]] = p++;  
27.         rep(i,len-1) {  
28.             if(sa[i]+d < len && sa[i+1]+d <len &&r[sa[i]] == r[sa[i+1]]&& r[sa[i]+d] == r[sa[i+1]+d])  
29.                 k[sa[i+1]] = p - 1;  
30.             else k[sa[i+1]] = p++;  
31.         }  
32.         if(p >= len) return ;  
33.         d *= 2,up = p, p = 0;  
34.     }  
35. }  
36. void getHeight(int len) {  
37.     rep(i,len) rk[sa[i]] = i;  
38.     height[0] =  0;  
39.     for(int i = 0,p = 0; i < len - 1; i++) {  
40.         int j = sa[rk[i]-1];  
41.         while(i+p < len&& j+p < len&& w[i+p] == w[j+p]) {  
42.             p++;  
43.         }  
44.         height[rk[i]] = p;  
45.         p = max(0,p - 1);  
46.     }  
47. }  
48. int getSuffix(char s[]) {  
49.     int len = strlen(s),up = 0;   
50.     for(int i = 0; i < len; i++) {  
51.         w[i] = s[i];  
52.         up = max(up,w[i]);  
53.     }  
54.     w[len++] = 0;  
55.     getSa(len,up+1);  
56.     getHeight(len);  
57.     return len;  
58. }  

最后给出几题习题

•Poj 2774 –最长公共连续子串，入门题目

•Poj1743—最长不重叠重复子串

•Hint：二分的判定要小心点，这题有点特别。

•Poj3294—出现次数超过一半的最长子串

•Hint：判断组中不同串出现次数的技巧很关键

•Poj3261—重复k次可重叠子串。

•Hint：会了上面两题，这题应该很简单，可以试试用单调栈。

•Poj2758—后缀数组+rmq

•Hint：这题难度不在rmq，而在于写代码的能力和查询的算法实现。