基本上一搜后缀数组网上的模板都是《后缀数组——处理字符串的有力工具》这一篇的注释,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 }
当然初次接触的人看这代码必然不知这是什么鬼,建议可以多找几篇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);
总觉得把整个代码和注释混一起看好烦躁,分解开来一部分一部分看吧~
定义:
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 */