后缀数组

【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自动机、后缀数组等的应用问题，算了，以后再搞吧囧。

 

原文转自：http://www.cppblog.com/MatoNo1/archive/2011/10/23/158926.html