【后缀数组之SA数组】【真难懂啊】

Posted on 2015-08-28 11:24  LLGemini  阅读(1562)  评论(0编辑  收藏  举报

基本上一搜后缀数组网上的模板都是《后缀数组——处理字符串的有力工具》这一篇的注释,O(nlogn)的复杂度确实很强大,但对于初次接触(比如窝)的人来说理解起来也着实有些困难(比如窝就活活好了两天的光阴。。),看了那么多材料感觉《挑战程序设计》的后缀数组解释理解起来会相对容易很多,然而它的复杂度是O(nlog2n)的,主要区别是对子串排序的时候前者用了计数排序--O(n),而后者用了快排--O(nlogn),这就导致了最终的复杂度后者比前者多了一个O(logn)

 

O(nlog2n)算法

先附清爽版求SA(Suffix_Array)模板,主要思想当然仍是倍增法;

 1 /*0(nlog(n)^2)*/
 2 #include <iostream>
 3 #include <cstring>
 4 #include <cstddef>
 5 #include <cstdio>
 6 #include <string>
 7 #include <algorithm>
 8 using namespace std;
 9 const int MAXN = 10001;
10 int n,k;
11 int rank[MAXN+1],tmp[MAXN+1];
12 
13 bool comp_sa(int i, int j)
14 {
15     if(rank[i] != rank[j])
16         return rank[i] < rank[j];
17     int ri = i+k <= n? rank[i+k] : -1;
18     int rj = j+k <= n? rank[j+k] : -1;
19     return ri < rj;
20 }
21 
22 void calc_sa(string &S, int *sa) //计算字符串S的后缀数组
23 {
24     n = S.size();
25     //初始长度为1
26     for(int i = 0; i <= n; i++)
27     {
28         sa[i] = i;
29         rank[i] = i < n ? S[i] : -1;
30     }
31 
32     for( k = 1; k <= n; k *= 2)
33     {
34         sort(sa,sa+n+1,comp_sa); //双关键字快排
35 
36         //先在tmp中临时存储新计算的rank,再转存回rank中
37         tmp[sa[0]] = 0;
38         for(int i = 1; i <= n; i++)
39         {
40             tmp[sa[i]] = tmp[sa[i-1]] + (comp_sa(sa[i-1],sa[i]) ? 1: 0);
41         }
42         for(int i = 0; i <= n; i++)
43         {
44             rank[i] = tmp[i];
45         }
46     }
47 }
48 
49 int main()
50 {
51     string S = "abracadabra";
52     int *sa = new int[S.size()+1];
53      SuffixArrayMatch(S,sa,T);
54      delete [] sa;
55      sa = NULL;
56 }
View Code

 

当然初次接触的人看这代码必然不知这是什么鬼,建议可以多找几篇blog的代码注释(有很多大神的注释很细致很不错),加上手动模拟理解一下,多看几遍肯定会开窍的~;

定义(理解!):

后缀数组Suffix_Array:SA[]

  将某个字符串的所有后缀按字典序排序后得到的数组;

  SA[i] = k即表示排序后第i小的子串为S[k ... n](为好理解暂用字符串下标1~n);

  而最终我们要达到的SA数组状态如下例所示(摘自罗神的PPT):

名次数组:Rank[]

  保存后缀S[i ... n]在排序中的名次;

  rank[i] = k即表示子串S[i ... n]在所有后缀的字典序排序中为第k小;

通过图1可以看出SA[1] = 4, Rank[4] = 1; 即SA数组与Rank数组为互逆关系

 

 

后缀数组的计算

算法的基本思想 -- 倍增

什么叫倍增?

  首先计算每个位置开始的长度为1的子串的顺序,利用该结果计算长度为2的子串的顺序,再利用长度为2的子串的顺序结果计算长度为4的子串的顺序 ...... 不断倍增,知道长度大于等于原字符串长度,就得到了后缀数组;

  要计算长度为2的子串的顺序,只要排序两个字符组成的数对即可。

  比如原字符串 aaba (下面各字母的下标代表该字母在原字符串的位置)

  a1=a2=a4 < b3,那么a1a2一定小于a2b3

  要求长度为2k的子串的顺序,只要知道长度为k的子串的顺序即可。

  比如原字符串 aabac

  a1a2 < a2b3 < a4c5< b3a4,那么 a1a2b3a4 一定小于 a2b3a4c5 

  记rankk(i)为S[i, k](从i开始的长度为k的子串)在所有排好序的长度为k的子串中是第几小的;

  要计算长度为2k的子串的顺序,就只要对两个rank组成的数对进行排序即可。通过rankk(i)与rankk(i+k)的数对和rankk(j)与rankk(j+k)的数对比较(双元素比较)来代替对S[i, 2k]和S[j, 2k]的直接比较。因为比较rankk(i)rankk(j)就相当于比较S[i, k]与S[j, k],比较rankk(i+k)与rankk(j+k)就相当于比较S[i+k, k]与S[j+k, k]。

 

举个例子: abracadabra

初始化:SA[i] = i;

        rank[i] = S[i]; //初始的时为对字符串中的单个字符排序,故可以直接将rank初始为字符的ASCII码,注意此时的rank并非实际意义上的排序,仅仅是相对排序,即S[i]>S[j]则rank[i] > rank[j]而已;

    还有需要注意的一点是此处有一个处理字符串的小技巧是将SA[n] 定义为 -1;这样以后的rank值就可以从1开始排了。

k = 0,初始化,得到S[i, 1]的排序;

  sa[0] : 0         rank[0] : 97
  sa[1] : 1         rank[1] : 98
  sa[2] : 2         rank[2] : 114
  sa[3] : 3         rank[3] : 97
  sa[4] : 4         rank[4] : 99
  sa[5] : 5         rank[5] : 97
  sa[6] : 6         rank[6] : 100
  sa[7] : 7         rank[7] : 97
  sa[8] : 8         rank[8] : 98
  sa[9] : 9         rank[9] : 114
  sa[10] : 10      rank[10] : 97
  sa[11] : 11  rank[11] : -1

k = 1;得到S[i, 2]的排序;

  sa[0] : 11               rank[11] : 0
  sa[1] : 10       a      rank[10] : 1
  sa[2] : 0         ab     rank[0] : 2
  sa[3] : 7         ab     rank[7] : 2
  sa[4] : 3         ac     rank[3] : 3
  sa[5] : 5         ad     rank[5] : 4
  sa[6] : 1         br     rank[1] : 5
  sa[7] : 8         br     rank[8] : 5
  sa[8] : 4         ca     rank[4] : 6
  sa[9] : 6         da     rank[6] : 7
  sa[10] : 2        ra     rank[2] : 8
  sa[11] : 9        ra     rank[9] : 8

k = 2;得到S[i, 4]的排序;

  sa[0] : 11                rank[11] : 0
  sa[1] : 10        a        rank[10] : 1
  sa[2] : 0         abra      rank[0] : 2
  sa[3] : 7         abra      rank[7] : 2
  sa[4] : 3         acad      rank[3] : 3
  sa[5] : 5         adab      rank[5] : 4
  sa[6] : 8         bra        rank[8] : 5
  sa[7] : 1         brac      rank[1] : 6
  sa[8] : 4         cada      rank[4] : 7
  sa[9] : 6         dabr      rank[6] : 8
  sa[10] : 9        ra         rank[9] : 9
  sa[11] : 2        raca      rank[2] : 10

k = 4;得到S[i, 8]的排序;

  sa[0] : 11                  rank[11] : 0
  sa[1] : 10        a          rank[10] : 1
  sa[2] : 7         abra        rank[7] : 2
  sa[3] : 0         abracada     rank[0] : 3
  sa[4] : 3         acadabra     rank[3] : 4
  sa[5] : 5         adabra        rank[5] : 5
  sa[6] : 8         bra          rank[8] : 6
  sa[7] : 1         bracadab     rank[1] : 7
  sa[8] : 4         cadabra       rank[4] : 8
  sa[9] : 6         dabra          rank[6] : 9
  sa[10] : 9         ra              rank[9] : 10
  sa[11] : 2         racadabr      rank[2] : 11

k = 8;得到S[i, n]的排序;

      sa[0] : 11                   rank[11] : 0
  sa[1] : 10    a            rank[10] : 1
  sa[2] : 7      abra       rank[7] : 2
  sa[3] : 0      abracadabra   rank[0] : 3
  sa[4] : 3      acadabra      rank[3] : 4
  sa[5] : 5      adabra       rank[5] : 5
  sa[6] : 8      bra         rank[8] : 6
  sa[7] : 1      bracadabra     rank[1] : 7
  sa[8] : 4      cadabra      rank[4] : 8
  sa[9] : 6      dabra       rank[6] : 9
  sa[10] : 9    ra           rank[9] : 10
  sa[11] : 2    racadabra       rank[2] : 11

这样就得到了SA数组;

对于字符串的排序采用双关键字快排,复杂度为O(nlogn);每次计算后缀S[i ... n] 的 rank值时,先与该后缀排序后的前一个后缀比较,如果相等则rank值相同,否则rank值加1;

1 //先在tmp中临时存储新计算的rank,再转存回rank中
2 tmp[sa[0]] = 0;
3 for(int i = 1; i <= n; i++)
4     tmp[sa[i]] = tmp[sa[i-1]] + (comp_sa(sa[i-1],sa[i]) ? 1: 0);
5 
6 for(int i = 0; i <= n; i++)
7     rank[i] = tmp[i];

 

 

O(nlogn)算法

声明:本部分大部分内容摘自staginner博文,写得很详细很明白~推荐!

还是先附清爽版代码

 1 #include <iostream>
 2 #include <cstring>
 3 #include <cstddef>
 4 #include <cstdio>
 5 #include <string>
 6 #include <algorithm>
 7 
 8 int wa[maxn],wb[maxn],wv[maxn],ws[maxn];
 9 int cmp(int *rank, int a,int b,int l)
10 {
11     return rank[a]==rank[b] && rank[a+l]==rank[b+l];
12 }
13 void DA(int *r,int *sa,int n,int m)
14 {
15     int i, k, p, *x=wa, *y=wb, *t;
16 
17     for(i=0;i<m;i++) ws[i] = 0;
18     for(i=0;i<n;i++) ws[x[i] = r[i]]++;
19     for(i=1;i<m;i++) ws[i] += ws[i-1];
20     for(i=n-1;i>=0;i--) sa[--ws[x[i]]] = i;
21   
22     for(k=1, p=1; p<n; k*=2, m=p)
23     {
24         p=0; for(i=n-k;i<n;i++) y[p++]=i;
25         for(i=0;i<n;i++) if(sa[i]>=k) y[p++]=sa[i]-k;
26   
27         for(i=0;i<n;i++) wv[i]=x[y[i]];
28 
29         for(i=0;i<m;i++) ws[i]=0;
30         for(i=0;i<n;i++) ws[wv[i]]++;
31         for(i=1;i<m;i++) ws[i]+=ws[i-1];
32         for(i=n-1;i>=0;i--) sa[--ws[wv[i]]]=y[i];
33 
34         for(t=x,x=y,y=t,p=1,x[sa[0]]=0,i=1;i<n;i++)
35             x[sa[i]]=cmp(y,sa[i-1],sa[i],k)?p-1:p++;
36     }
37     return;
38 }
39 
40 DA(r,sa,n+1,128);
View Code

总觉得把整个代码和注释混一起看好烦躁,分解开来一部分一部分看吧~

定义:

int wa[maxn],wb[maxn],wv[maxn],ws[maxn];
/*
wa[]: 本意是保存各个后缀的rank值的,但是这里并没有去存储rank值,因为后续只是涉及wa[]的比较工作,
      因而这一步可以不用存储真实的rank值,能够反映相对的大小即可。
wb[]: 存放的是按第二关键字排序的子串首字符下标
wv[]: 存放每个子串的第一关键字
ws[]: 存放每个rank值的数目
*/

函数:

void da(int *r,int *sa,int n,int m)
{
        ...
}
/*
*r:  字符串(数组)
*sa: 后缀数组
n:   字符串中字符的个数,注意这里的n里面是包括人为在字符串末尾添加的那个0的
m:   字符串中字符的取值范围,是基数排序的一个参数,如果原序列都是字母可以直接取128,如果原序列本身都是整数的话,则m可以取比最大的整数大1的值。
*/
int i, k, p, *x=wa, *y=wb, *t;
/*
*x 代替wa数组,*y 代替wb数组
*t 作交换指针
*/

 

以下四行代码是把长度为1的子串进行基数排序 Tips: 如果不理解为什么这样可以达到计数排序的效果,建议自己实际用纸笔模拟一下!

for(i=0;i<m;i++) ws[i] = 0;
for(i=0;i<n;i++) ws[x[i] = r[i]]++;
for(i=1;i<m;i++) ws[i] += ws[i-1];
for(i=n-1;i>=0;i--) sa[--ws[x[i]]] = i;
/*
ws[]: 存放每个rank值的数目; 第1行,清零; 第2行,上面已经提到x[]保存的是后缀的相对rank值,x[i]=r[i]的意思是将x[i]初始化为各字符的ASCII值,字符的ASCII值也就可以代表长度为1的子串的相对顺序; 第3行,求出最后一个子串i的rank是多少,供第4行使用 第4行,相当于从后向前得到各子串的sa[]数组,i之所以从n-1开始循环,是为了保证在当字符串中有相等的字符串时,默认靠前的字符串更小一些。 */

注意理解上述“相对顺序”的含义,上一部分也提到过这个概念,我们并不需要求出子串的排序到底是1还是2,我们只需要达到若r[i]<r[j],则x[i]<x[j]的要求即可。这也解释了利用字符的ASCII值初始x[]数组的合理性;

主要循环:

for(k=1, p=1; p<n; k*=2, m=p)
{
    ...    
}
/*
这层循环中p的值表示的是此时关键字不同的子串的数量,也可以这么理解,所有子串排序后,相等的子串rank值相同,则关键字的范围是[1,p];
如果p达到n,即各后缀的rank与sa 已全部求出,因为后缀长度不一,所以不可能出现相等的情况;
k代表当前待合并的字符串的长度,每次将两个长度为k的字符串合并成一个长度为2k的字符串——倍增思想;
m同样代表基数排序的元素的取值范围
*/

以下两行代码实现对第二关键字的排序

此处我们需要借助罗神的插图来理解了

我们可以这么理解所谓第二关键字 (为避免混淆,原字符串r换为字母s来表示)

首字符为i长度为1的子串没有第二关键字

首字符为i长度为2的子串s[i,2]的第一关键字为s[i,1]的rank值,第二关键字为s[i+1,1]的rank值;

首字符为i长度为4的子串s[i,4]的第一关键字为s[i,2]的rank值,第二关键字为s[i+2,2]的rank值;

...

首字符为i长度为2k的子串s[i,2k]的第一关键字为s[i,k]的rank值,第二关键字为s[i+k,k]的rank值;

 

由图2,可以看到首字符下标为n-k至n的子串的第二关键字都为0,因此如果按第二关键字排序,必然这些子串都是排在前面的。(第二关键字为0即无法构成以r[i]为首字符的长度为2k的子串,长度都不够,自然字典序会小咯)—— 第1个循环。

我们还可以看到,下面一行的第二关键字的值(非0)都是上一轮的rank值,且上一行中只有首字符下标(sa[i])>=k的子串的rank值才会作为下一行的首字符下标为sa[i]-j的子串的第二关键字,而且显然按sa[i]的顺序rank[sa[i]]是递增的(rank[sa[i]] == i) —— 第2个循环。

图2的y[]值为

k: 1时
y[0] : 8
y[1] : 7
y[2] : 0
y[3] : 2
y[4] : 3
y[5] : 4
y[6] : 5
y[7] : 1
y[8] : 6


k: 2时
y[0] : 7
y[1] : 8
y[2] : 6
y[3] : 1
y[4] : 2
y[5] : 3
y[6] : 4
y[7] : 5
y[8] : 0

y[i] = x的含义是按第二关键字排序后第i小的子串首字符下标为x;

记住y[]里存放的是按第二关键字排序的子串首字符下标!

 1 p=0; for(i=n-k;i<n;i++) y[p++]=i;
 2 for(i=0;i<n;i++) if(sa[i]>=k) y[p++]=sa[i]-k;
for(i=0;i<n;i++) wv[i]=x[y[i]];
/*这里相当于提取出每个子串的第一关键字(前面说过了x[]是保存上一轮的rank值的,也就是子串的第一关键字),放到wv[]里面是方便后面的使用*/

 

以下四行代码是按第一关键字进行的基数排序

wv[]: 存放每个长度为k的子串的第一关键字,wv[i] = x的含义为按第二关键字第i小的子串的第一关键字的值;
ws[]: 存放每个关键字的个数

for(i=0;i<m;i++) ws[i]=0;
for(i=0;i<n;i++) ws[wv[i]]++;
for(i=1;i<m;i++) ws[i]+=ws[i-1];
for(i=n-1;i>=0;i--) sa[--ws[wv[i]]]=y[i];
/*i之所以从n-1开始循环,含义同上,同时注意这里是y[i],因为y[i]里面才存着字符串的下标*/

 最后一行巧妙地将第一关键字与第二关键字结合起来了,注意理解!

 

下面两行就是计算合并之后的rank值了,而合并之后的rank值应该存在x[]里面,但我们计算的时候又必须用到上一层的rank值,也就是现在x[]里面放的东西,如果我既要从x[]里面拿,又要向x[]里面放,怎么办?
当然是先把x[]的东西放到另外一个数组里面,省得乱了。这里就是用交换指针的方式,高效实现了将x[]的东西“复制”到了y[]中。

t=x,x=y,y=t; 
for(p=1,x[sa[0]]=0,i=1;i<n;i++)
    x[sa[i]]=cmp(y,sa[i-1],sa[i],k)?p-1:p++;

/*
  这里就是用x[]存储计算出的各字符串rank的值了,记得我们前面说过,计算sa[]值的时候如果字符串相同是默认前面的更小的,但这里计算rank的时候必须将相同的字符串看作有相同的rank,要不然p==n之后就不会再循环啦
  p的值表示的是此时关键字不同的串的数量
  cmp比较函数,合并的子串相同则返回1,不同返回0;
  注意p和i的初始值需为1,因为循环中存在i-1和p-1,而x[sa[0]]的值也需初始化为0
*/

 

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