[对罗穗骞论文的一点理解]后缀数组——倍增算法

基本定义

首先先来看几个定义：
字符串的大小比较
从第一位开始，按照字典序比较大小，如果相同则比较后一位，否则按照当前这一位的比较结果确定结果，如果比较的时候这一位已经超出其中一个串的长度，则比较它们的长度，长的一个串比较大，如果长度相等，那么它们也相等。
后缀数组
用SA表示，它是一个一维数组，保存的是所有后缀的排序，每一位表示的就是按照大小排序的这个后缀是从哪一位开始的，也可以表示为S[i..N]
名次数组
用Rank表示，表示的是第i个后缀的排名
简单地说，后缀数组存的就是排第i的是谁，名次数组存的就是第i个后缀排第几,借此，容易得出，这两个储存的值是互逆的。
性质：
在求出了Rank之后，就可以用O(1)的时间来比较两个后缀的大小。
同样，如果求出了这两个数组的任意一个，就可以利用O(N)的时间来求出另外一个数组
设一个字符串的长度为N,为了方便比较大小，可以在字符串的结尾处加一个比前面所有字符都小的一个字符。这样，在任意两个后缀中,直接比较这两个后缀的大小，最多也只需要比较N次，也就是在比较第N个字符之后肯定会得到结果

倍增算法的实现

主要思路：
利用倍增的思想，对每一个字符开始的长度为2k的子字符串进行排序，得到Rank数组
k从0开始，每一次加上1，当2k>N的时候，每一个这样的子字符串都相当于是一个后缀，而从开始，已经不断将他们进行了比较，同时会满足Rank中不会有重复的值。
这时，Rank的值就是最终所求的结果
从0开始，每一次k+1，在排序的时候，可以利用k−1的Rank的两个值，通过排序决定当前这一位的排名。
利用这样类似RMQ-ST的做法，就可以求出最终的Rank
以下借用一张罗穗骞神犇论文中的一张图，来表示倍增算法运行的过程。

那么，对于排序来说，应该用什么方法呢？在论文中，提到了一种叫做基数排序的神奇的排序方法。对于这个算法，理论上是可以只跑O(N)的时间的，但事实上对于一些数据来说可能会更慢一些。当然，在大多数的情况下，速度还是很快的。
来看一看代码

void da(int *r, int *sa, int n, int m) {    int i, j, p, *x = wa, *y = wb, *t;    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;    }//基数排序    for(j = 1, p = 1; p < n; j *= 2, m = p) { //倍增        for(p = 0, i = n - j; i < n; i++) {            y[p++] = i;        }        for(i = 0; i < n; i++) {            if(sa[i] >= j) y[p++] = sa[i] - j;        } // 排序第二关键字 借用之前的sa        for(i = 0; i < n; i++) {            wv[i] = x[y[i]];        }        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];        } // 基数排序 排序第一关键字        for(t = x, x = y, y = t, p = 1, x[sa[0]] = 0, i = 1; i < n; i++) {            x[sa[i]] = cmp(y, sa[i - 1], sa[i], j) ? p - 1 : p++;        }    }    return;}

代码有点长，因为是自己重打一遍的原因，而且代码风格比较不一样，有很多空行，实际上真正有意义的也就只有大概30行左右待排序的字符串存在r中，长度为N,而且最大值小于m，当函数结束后，放在sa一开始，先对长度为1的字符串进行排序，同样利用了基数排序。注意，如果r的最大值很大，那么采用快排会更快一些以下为基数排序部分的代码

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;}

这里的x数组保存的就是Rank的值，但是在求的过程中，并不需要直接求出来，所以x只是用来比较字符串的大小的。
接下来要进行若干次基数排序，在这里有一个优化。
基数排序的过程分为两次，第一次是对第一关键字排序，第二次是对第二关键字排序。对于第二次排序，其实可以用上一次的求出的sa数组直接算出，而没有必要去重新算一次。

for(p = 0, i = n - j; i < n; i++) {    y[p++] = i;}for(i = 0; i < n; i++) {    if(sa[i] >= j) y[p++] = sa[i] - j;}

这里的j是字符串的长度，数组y保存的则是对第二关键字排序的结果，接着就要对第一关键字进行排序。

        for(i = 0; i < n; i++) {            wv[i] = x[y[i]];        }        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];    }

求出新的sa的值以后，就要去计算Rank的值，
求Rank的时候，主要要注意，会有多个字符串的Rank是相同的，那么就要直接判断两个字符串是否相等，再来得到结果。
这里有两个在空间和时间上的优化：
因为y数组已经没有了用处，可以用来储存Rank的值，同时，如果交换两个数组，不需要一个元素一个元素去交换，可以直接交换两个数组的指针，就完成了交换。代码只有三行，十分简单

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

其中cmp函数的代码如下

int cmp(int *r, int a, int b, int l) {    return r[a] == r[b] && r[a + l] == r[b + l];}

可以先模拟一遍，这样就更容易理解原文中的一段话。
如果ra=rb，就说明分别以它们开头长度为1的两个字符串中都不包括rn−1，所以调用ra+1、rb+1都不会导致数组的下标越界，最多只会停留在最后一个字符，这样就不再需要另外作特殊的判断。
执行完所有以上代码之后，在这时，Rank的值已经储存在x当中，而p保存了不同的字符串的个数。
这里可以加上一个优化,如果p=n，就证明所有字符串中没有相同的字符串，也就是说接下来的排序已经是没有意义的，不会改变Rank的值。
那么，循环的初始化和结束就写成了一开始的代码那样，这里再看一次。

for(j = 1, p = 1; p < n; j *= 2, m = p) {......}

最后再来看看时间复杂度，倍增部分最多执行O(log2n)，而基数排序需要O(n)的时间，最终时间为O(nlog2n)
原论文中还提到了一种叫做dc3的算法，时间复杂度仅为O(N)，但是太复杂，所以没有写。