【后缀数组真难懂啊啊……就20+行的代码搞了几天才理解……不知是不是我太沙茶了】
【1】一些定义:
字符串:广义的字符串是指“元素类型有序,且元素值有一定范围的序列”,其元素不一定非要是字符,可以是数字等,因此整数、二进制数等也是字符串;
字符集:字符串的元素值的范围称为字符集,其大小记为SZ。
字符串的长度:字符串中元素的个数,一般记为N,长度为N的字符串A第一次提到时一般用A[0..N-1]来表示;
前缀:字符串A[0..N-1]的从A[0]开始的若干个连续的字符组成的字符串称为A的前缀,以下“前缀i”或者“编号为i的前缀”指的都是A[0..i];
后缀:字符串A[0..N-1]的到A[N-1]终止的若干个连续的字符组成的字符串称为A的后缀,以下“后缀i”或者“编号为i的后缀”指的都是A[i..N-1];
对于一个长度为N的字符串,将其N个后缀按字典序大小进行排序,得到两个数组sa[i]和rank[i],sa[i]为排在第i位的后缀的编号(也就是一般说的ord[i]),rank[i]为排在后缀i排在的位置(称为后缀i的名次)。sa、rank值的范围均为[0..N-1]。sa和rank互逆,即sa[i]=j等价于rank[j]=i,或者说成sa[rank[i]]=rank[sa[i]]=i。这里,sa称为
后缀数组,rank称为
名次数组。
【2】用倍增算法求后缀数组:
在论文里,后缀数组有两种求法:倍增算法和DC3算法,前者的时间复杂度为O(NlogN),但常数较小,后者的时间复杂度为O(N),但常数较大,在实际应用中,两者的总时间相差不大,且后者比前者难理解得多(本沙茶理解前者都用了几天时间……后者就木敢看了)。这里就总结一下倍增算法吧囧……
首先,贴一下本沙茶的用倍增算法求后缀数组的模板:
void suffix_array()
{
int p, v0, v1, v00, v01;
re(i, SZ) S[i] = 0;
re(i, n) rank[i] = A[i];
re(i, n) S[A[i]]++;
re2(i, 1, SZ) S[i] += S[i - 1];
rre(i, n) sa[--S[A[i]]] = i;
for (int j=1; j<n; j<<=1) {
p = 0; re2(i, n-j, n) tmp[p++] = i;
re(i, n) if (sa[i] >= j) tmp[p++] = sa[i] - j;
re(i, SZ) S[i] = 0;
re(i, n) S[rank[i]]++;
re2(i, 1, SZ) S[i] += S[i - 1];
rre(i, n) sa[--S[rank[tmp[i]]]] = tmp[i];
tmp[sa[0]] = p = 0;
re2(i, 1, n) {
v0 = sa[i - 1]; v1 = sa[i];
if (v0 + j < n) v00 = rank[v0 + j]; else v00 = -1;
if (v1 + j < n) v01 = rank[v1 + j]; else v01 = -1;
if (rank[v0] == rank[v1] && v00 == v01) tmp[sa[i]] = p; else tmp[sa[i]] = ++p;
}
re(i, n) rank[i] = tmp[i];
SZ = ++p;
}
}
这里A是待求sa和rank的字符串。
<1>倍增算法的思想:
记R[i][j]为A[i..i+2
j-1](如果越界,则后面用@填充)在A的所有长度为2
j的子串(越界则后面用@填充)中的名次(rank)值。倍增算法就是按阶段求出所有R[i][j]的值,直到2
j>N为止。首先,R[i][0]的就是字符A[i]在A[0..N-1]中的名次,是可以直接用计数排序来实现的。然后,若R[0..N-1][j-1]已知,则可以按照以下方法求出R[0..N-1][j]的值:对每个i(0<=i<N),构造一个二元组<X
i, Y
i>,其中X
i=R[i][j-1],Y
i=R[i+2
j][j-1](若i+2
j>=N,则Y
i=-∞),然后对这N个二元组按照第一关键字为X,第二关键字为Y(若两者都相等则判定为相等)进行排序(可以用基数排序来实现),排序后,<X
i, Y
i>的名次就是的R[i][j]的值。
<2>一开始,对A中的各个字符进行计数排序:
re(i, SZ) S[i] = 0;
re(i, n) rank[i] = A[i];
re(i, n) S[A[i]]++;
re2(i, 1, SZ) S[i] += S[i - 1];
rre(i, n) sa[--S[A[i]]] = i;
这个木有神马好说的,在搞懂了基数排序之后可以秒掉。唯一不同的是这里加了一句:rank[i]=A[i],这里的rank[i]是初始的i的名次,MS不符合rank[i]的定义和sa与rank间的互逆性。这里就要解释一下了囧。因为在求sa的过程中,rank值可能不符合定义,因为长度为2
j的子串可能会有相等的,此时它们的rank值也要相等,而sa值由于有下标的限制所以不可能有相等的。因此,在过程中,rank其实是用来代替A的子串的,这样rank值只需要表示一个“相对顺序”就行了,也就是:rank[i0]>(=, <)rank[i1],当且仅当A[i0..i0+2
j-1]>(=, <)A[i1..i1+2
j-1]。这样,可以直接将A[i]值作为初始的rank[i]值。
<3>j(代替2
j)的值从1开始不断倍增,对二元组进行基数排序求出新阶段的sa值:
for (int j=1; j<n; j<<=1) {
p = 0; re2(i, n-j, n) tmp[p++] = i;
re(i, n) if (sa[i] >= j) tmp[p++] = sa[i] - j;
re(i, SZ) S[i] = 0;
re(i, n) S[rank[i]]++;
re2(i, 1, SZ) S[i] += S[i - 1];
rre(i, n) sa[--S[rank[tmp[i]]]] = tmp[i];
注意这个基数排序的过程是很特别的。首先,它并不是对A在进行排序,而是对上一阶段求出的rank在进行排序。因为前面已经说过,在求sa的过程中,rank就是用来代替A的对应长度的子串的,由于不能直接对子串进行排序(那样的话时间开销很恐怖的),所以只能对rank进行排序。另外,这里在对二元组<x, y>的第二关键字(y)进行排序的过程中加了优化:这些y其实就是把上一阶段的sa整体左移了j,右边空出的部分全部用@(空串)填充得到的,由于空串的字典序肯定最小,因此将右边的空串按照下标顺序先写入临时sa(代码中用tmp表示的就是临时sa,也就是对第二关键字y排序后的ord结果),然后,上一阶段的sa如果左移后还木有消失的(也就是sa值大于等于j的),再按顺序写入临时sa,就得到了排序结果。剩下的对x的排序结果就是上一阶段的sa,唯一不同的是对于x相同的,按照临时名次递增的顺序。
<4>求出新阶段的rank值:
tmp[sa[0]] = p = 0;
re2(i, 1, n) {
v0 = sa[i - 1]; v1 = sa[i];
if (v0 + j < n) v00 = rank[v0 + j]; else v00 = -1;
if (v1 + j < n) v01 = rank[v1 + j]; else v01 = -1;
if (rank[v0] == rank[v1] && v00 == v01) tmp[sa[i]] = p; else tmp[sa[i]] = ++p;
}
re(i, n) rank[i] = tmp[i];
SZ = ++p;
由于下一阶段需要使用本阶段的rank值,因此在求出了本阶段的sa值以后,需要求rank值。(代码中的tmp起了临时rank的作用,目的是节省空间)
因为sa值已经求出,因此只要依次扫描sa就可以得到rank值,唯一要做的工作就是找到哪些子串是相等的,它们的rank值应该相等,除此之外,rank值只要依次加1即可。判定相等的方法:只需判定rank[i]和rank[i+j]是否都对应相等即可。若rank[i+j]越界,用-∞(当然任何一个负数都行,代码中用了-1)来表示。
最后还有一个优化:由于本阶段的名次的范围只有[0..p]这么多,下一阶段的“字符集”(其实就是rank集)的大小SZ可以设为p+1,这样可以省一些时间。
这样后缀数组sa和名次数组rank就全部求完了。
以后还有一些更重要的东东就是AC自动机、后缀数组等的应用问题,算了,以后再搞吧囧。