后缀数组--学习笔记（倍增算法）

这几天看了小罗的论文，感觉太神了。

以下先添上一段Staginner大牛的读后笔记：

文章链接：点击打开链接

int wa[maxn],wb[maxn],wv[maxn],ws[maxn];
int cmp(int *r,int a,int b,int l)
{return r[a]==r&&r[a+l]==r[b+l];}  //就像论文所说，由于末尾填了0，所以如果r[a]==r[b]（实际是y[a]==y[b]），说明待合并的两个长为j的字符串，前面那个一定不包含末尾0，因而后面这个的起始位置至多在0的位置，不会再靠后了，因而不会产生数组越界。
//da函数的参数n代表字符串中字符的个数，这里的n里面是包括人为在字符串末尾添加的那个0的，但论文的图示上并没有画出字符串末尾的0。
//da函数的参数m代表字符串中字符的取值范围，是基数排序的一个参数，如果原序列都是字母可以直接取128，如果原序列本身都是整数的话，则m可以取比最大的整数大1的值。
void da(int *r,int *sa,int n,int m)
{
int i,j,p,*x=wa,*y=wb,*t;
//以下四行代码是把各个字符（也即长度为1的字符串）进行基数排序，如果不理解为什么这样可以达到基数排序的效果，不妨自己实际用纸笔模拟一下，我最初也是这样才理解的。
for(i=0;i<m;i++) ws[i]=0;
for(i=0;i<n;i++) ws[x[i]=r[i]]++;  //x[]里面本意是保存各个后缀的rank值的，但是这里并没有去存储rank值，因为后续只是涉及x[]的比较工作，因而这一步可以不用存储真实的rank值，能够反映相对的大小即可。
for(i=1;i<m;i++) ws[i]+=ws[i-1];
for(i=n-1;i>=0;i--) sa[--ws[x[i]]]=i;  //i之所以从n-1开始循环，是为了保证在当字符串中有相等的字符串时，默认靠前的字符串更小一些。
//下面这层循环中p代表rank值不用的字符串的数量，如果p达到n，那么各个字符串的大小关系就已经明了了。
//j代表当前待合并的字符串的长度，每次将两个长度为j的字符串合并成一个长度为2*j的字符串，当然如果包含字符串末尾具体则数值应另当别论，但思想是一样的。
//m同样代表基数排序的元素的取值范围
for(j=1,p=1;p<n;j*=2,m=p)
{
//以下两行代码实现了对第二关键字的排序
for(p=0,i=n-j;i<n;i++) y[p++]=i;  //结合论文的插图，我们可以看到位置在第n-j至n的元素的第二关键字都为0，因此如果按第二关键字排序，必然这些元素都是排在前面的。
for(i=0;i<n;i++) if(sa[i]>=j) y[p++]=sa[i]-j;  //结合论文的插图，我们可以看到，下面一行的第二关键字不为0的部分都是根据上面一行的排序结果得到的，且上一行中只有sa[i]>=j的第sa[i]个字符串（这里以及后面指的“第?个字符串”不是按字典序排名来的，是按照首字符在字符串中的位置来的）的rank才会作为下一行的第sa[i]-j个字符串的第二关键字，而且显然按sa[i]的顺序rank[sa[i]]是递增的，因此完成了对剩余的元素的第二关键字的排序。
//第二关键字基数排序完成后，y[]里存放的是按第二关键字排序的字符串下标
for(i=0;i<n;i++) wv[i]=x[y[i]];  //这里相当于提取出每个字符串的第一关键字（前面说过了x[]是保存rank值的，也就是字符串的第一关键字），放到wv[]里面是方便后面的使用
//以下四行代码是按第一关键字进行的基数排序
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[]中。
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++; //这里就是用x[]存储计算出的各字符串rank的值了，记得我们前面说过，计算sa[]值的时候如果字符串相同是默认前面的更小的，但这里计算rank的时候必须将相同的字符串看作有相同的rank，要不然p==n之后就不会再循环啦。
}
return;
}

//能够线性计算height[]的值的关键在于h[](height[rank[]])的性质，即h[i]>=h[i-1]-1，下面具体分析一下这个不等式的由来。
//论文里面证明的部分一开始看得我云里雾里，后来画了一下终于搞明白了，我们先把要证什么放在这：对于第i个后缀，设j=sa[rank[i] - 1]，也就是说j是i的按排名来的上一个字符串，按定义来i和j的最长公共前缀就是height[rank[i]]，我们现在就是想知道height[rank[i]]至少是多少，而我们要证明的就是至少是height[rank[i-1]]-1。
//好啦，现在开始证吧。
//首先我们不妨设第i-1个字符串（这里以及后面指的“第?个字符串”不是按字典序排名来的，是按照首字符在字符串中的位置来的）按字典序排名来的前面的那个字符串是第k个字符串，注意k不一定是i-2，因为第k个字符串是按字典序排名来的i-1前面那个，并不是指在原字符串中位置在i-1前面的那个第i-2个字符串。
//这时，依据height[]的定义，第k个字符串和第i-1个字符串的公共前缀自然是height[rank[i-1]]，现在先讨论一下第k+1个字符串和第i个字符串的关系。
//第一种情况，第k个字符串和第i-1个字符串的首字符不同，那么第k+1个字符串的排名既可能在i的前面，也可能在i的后面，但没有关系，因为height[rank[i-1]]就是0了呀，那么无论height[rank[i]]是多少都会有height[rank[i]]>=height[rank[i-1]]-1，也就是h[i]>=h[i-1]-1。
//第二种情况，第k个字符串和第i-1个字符串的首字符相同，那么由于第k+1个字符串就是第k个字符串去掉首字符得到的，第i个字符串也是第i-1个字符串去掉首字符得到的，那么显然第k+1个字符串要排在第i个字符串前面，要么就产生矛盾了。同时，第k个字符串和第i-1个字符串的最长公共前缀是height[rank[i-1]]，那么自然第k+1个字符串和第i个字符串的最长公共前缀就是height[rank[i-1]]-1。
//到此为止，第二种情况的证明还没有完，我们可以试想一下，对于比第i个字符串的字典序排名更靠前的那些字符串，谁和第i个字符串的相似度最高（这里说的相似度是指最长公共前缀的长度）？显然是排名紧邻第i个字符串的那个字符串了呀，即sa[rank[i]-1]。也就是说sa[rank[i]]和sa[rank[i]-1]的最长公共前缀至少是height[rank[i-1]]-1，那么就有height[rank[i]]>=height[rank[i-1]]-1，也即h[i]>=h[i-1]-1。
//证明完这些之后，下面的代码也就比较容易看懂了。
int rank[maxn],height[maxn];
void calheight(int *r,int *sa,int n)
{
int i,j,k=0;
for(i=1;i<=n;i++) rank[sa[i]]=i;  //计算每个字符串的字典序排名
for(i=0;i<n;height[rank[i++]]=k)  //将计算出来的height[rank[i]]的值，也就是k，赋给height[rank[i]]。i是由0循环到n-1，但实际上height[]计算的顺序是由height[rank[0]]计算到height[rank[n-1]]。
for(k?k--:0,j=sa[rank[i]-1];r[i+k]==r[j+k];k++);  //上一次的计算结果是k，首先判断一下如果k是0的话，那么k就不用动了，从首字符开始看第i个字符串和第j个字符串前面有多少是相同的，如果k不为0，按我们前面证明的，最长公共前缀的长度至少是k-1，于是从首字符后面k-1个字符开始检查起即可。
return;
}

//最后再说明一点，就是关于da和calheight的调用问题，实际上在“小罗”写的源程序里面是如下调用的，这样我们也能清晰的看到da和calheight中的int n不是一个概念，同时height数组的值的有效范围是height[1]~height
其中height[1]=0，原因就是sa[0]实际上就是我们补的那个0，所以sa[1]和sa[0]的最长公共前缀自然是0。
da(r,sa,n+1,128);
calheight(r,sa,n);

[b]首先可以很肯定的说，Staginner大牛对小罗的思想已经理解很深透了。因为之前我一直看不懂笔记中的很多地方，但这段时间来，我不断理解了这些话，反复体会大牛的笔记啊。
可以这么说，小罗的这段代码绝对是经过历练的，代码已经优化得相当完美，而这背后全都是由定理保证的。

首先小罗在他的数组后增添了一个元素0，比前面任何一个元素都来小，并且前面不能出现0。而且这个元素也看成原数组的一部分。

我主要来讲讲为什么这么做的原因？其他的只能看上面的代码来体会。

原数组下标是0开始的，起初有n个元素，那么最大下标为n-1，由于加了0，元素个数为n+1，最大下标为n，因此可以看到da数组的第三个参数应该传n+1。也就是0也被作为后缀排序了，那么显然rank
=0，名次从0开始。而且sa[0]=n，这些都是显然的。

最后我们会得到一个sa数组,sa[0]其实是添的那个0，真正有用的是sa[1....n]它分别对应着原数组下标0...n-1的一个全排列，因此求rank数组时，我们只需求rank[0...n-1],分别对应1...n的一个全排列，而rank
=0是显然的。

这些都没什么，小菜一碟吧，真正神奇的是......先看一下cmp函数，里面有r[a+l]==r[b+l],这里a+l和b+l会不会超过n呢？因为我们知道r
=0，而r[n+1...]都是未知数，这个如果越界后果不堪设想。但事实上，这个0添的极其巧妙，它确实导致了a+l和b+l不会超过n。可以简单的想，第一次基数排序后，rank
的值肯定是0.

接下来按照论文的操作，是补第二关键字的0，首先0被补一个0，那么就是00，而之前的第二关键字都是上一次

sa向斜下方挪下来的，第一关键字则直接由上面复制下来，而cmp里只有r[a]==r为真，才会比较r[a+l]==r[b+l]，而r[a]==r[b]即所谓的第一关键字相等，而0是不可能与前面任何值相等的，并且长度为1的后缀合并，能越界的位置只有0所在的位置，所以此时不会越界。并且0被作为前一个元素x的第二关键字，可以设想的是x0绝对是独一无二的，因此它的rank值也是独一无二的，它作为下一次基数排序（长度为2）的第一关键字还有00也算一个，他们也不可能和前面的值相等，而长度为2能够越界的也只有这两个位置，因此也不会越界。

[b]综上，可以由数学归纳法证得，详细的有兴趣的可以自己推导，对理解代码是有很大帮助的。

最后神奇的是calheight函数，起初我认为这个函数是有问题的，因为有一个定理h[i]>=h[i-1]-1，这里h[i]=height[rank[i]]，但其实这个定理是有前提的i>1&&rank[i]>1

而上面那段代码里完全没有体现出rank[i]=1时，对height数组做了什么操作。这个是让人无法理解。

首先它的代码是i从0循环到n-1，也就是rank[i]从1循环到n，height数组的下标是正确的。

因为定理是从下标1开始的，因此i>1,也就是i=1需特别计算。

而小罗代码下标从0开始，因此i=0需要特别计算，可以看到他对k赋了初值0，因此i=0，可以看到他是从后缀i和后缀j的首个字符开始比较，这是非常正确的。

但rank[i]=1的时候呢，实际上height[1]=0是定义的，而且这样定义不会影响这以后的height数组的计算，因为height[1]=0说明此时k=0，那么下一个height计算时，从后缀首字符开始比较，结果必然正确。但上面代码好像没有体现当rank[i]=1时，height[1]=0；但仔细观察代码：会发现当rank[i]=1时，j=sa[rank[i]-1]=sa[0]，这个是什么？别忘了是添加的0对吧。那么此时j=0，然后就像大牛的解释，公共前缀当然为0.

但事实并没有这么简单，为什么这么说呢？当rank[i]=1时，此时的k值为多少？如果k=3，那么两个后缀的比较起始位置都错了，更何况还越界了，这下程序不出错？然而经过大量测试数据表明，rank[i]=1时，k只能为1或0，而经过k?k--:0这句，k=0，那么比较位置从首字符开始，直接得到了height[1]=0；

好神奇，但这是为什么呢？这个需要证明吧，其实我尝试证明了，确实是这样的。

比如rank[3]=1，r[3]=2,首先可以肯定字符里没有小于2的（除了添上去的0），否则2开头的后缀不可能排第一，而k值是由i-1时的height[rank[i-1]]决定的，假设2前面是个x，x>=2，那么

k=height[rank[i-1]]=LCP(rank[i-1]-1,rank[i-1]);如果k>1，也就是比suffix(i-1)小一名的后缀与suffix(i-1)的最长公共前缀>=2，也就是以x2开头，因为名次比suffix(i-1)小1，所以一定会在某个字符出现前者小于后者的情况。

比如x2332x2333，在第5位出现2<3，但这显然不可能因为此时后缀2333不是排第一的，2332x2333才是排第一的，无论怎样都会出现这种情况，因此k<=1。

所以，从上面的代码来看，小罗是下过一番功夫优化的，不然写不出这么美妙的代码。当然这里省略了很多后缀数组和LCP的定理，有兴趣可以自己尝试证明！

感谢小罗这么强大的模板！我也尝试改造了一下，在数组前后均添0，然后在rank[i]=1时特判的代码，这样的好处是从下标1开始，不会被rank和sa的下标不一样绕来绕去，改造代码正在测试中！