倍增算法实现后缀数组详解+实现代码

阅读量：3898 次

发布时间：2019-05-23

本文共 5774 字，大约阅读时间需要 19 分钟。

【前言】

不要被文章的长度吓到，因为罗穗骞的论文要更长更详尽，我只是取了其中的一部分进行学习并做一个学习笔记，便于以后有需要的时候回顾。文章的内容主要是介绍后缀数组的实现，后缀数组的应用部分主要是结合例题来理解。

【后缀数组】

【引入】

AC自动机可以解决多模板匹配问题，但前提是事先知道所有的模板。在实际应用中，很多时候都是无法事先知道查询的，这时需要预处理文本串T，而不是模式串。

假定文本串为BANANA，可以在它的末尾加一个字符$(代表一个没在串中出现过的字符)，然后把它的所有后缀(BANANA$,ANANA$,NANA$,ANA$,NA$,A$)插入到一棵Trie中(前置知识：)，如下图所示：

有了这个后缀Trie，查找一个长度为m的模板时只需要进行一次普通Trie查找，时间复杂度为O(m)。

在实际运用中，会把后缀Trie中没有分支的链合并到一起，得到所谓的后缀树(suffix tree)，但由于后缀树的构造算法比较复杂难懂且容易写错，因此我们就需要用到后缀数组。后缀数组是后缀树的替代品，具有容易编程实现，占用内存空间小，且能够实现后缀树的很多功能而时间复杂度也并不逊色的特点。

在绘制上图时，我们需要把所有后缀按字典序从小到大排序，也就是说，我们需要把每个结点的所有子结点排好序，字典序小的在左边(规定$比所有其他字符都小)，然后每个叶结点里标上该后缀首字符的下标。比如后缀NA开始于BANANA的下标4(注意下标从0开始)，那么NA对应的叶结点标有"4"。为了方便描述，我们把“以下标k开头的后缀”表示为Suffix(k)，也就是说，对于文本串BANANA，Suffix(4)就是NA。

现在只需自左向右把所有叶子的编号排列出来，就可以得到后缀数组(suffix array)。

比如，BANANA的后缀数组SA[]={5,3,1,0,4,2}。

名次数组Rank[i]：保存的是Suffix(i)在所有后缀中从小到大排列的“名次”。

简单的说，后缀数组是“排第几的是谁？”，名次数组是“你排第几？”。容易看出，后缀数组和名次数组为互逆运算。

【后缀数组的实现】

根据定义，后缀数组可以直接通过一次快速排序得到，但是在最坏情况下，直接排序需要的时间是O(n^2logn)(虽然比较次数是O(nlogn)，但是两个字符串的比较是O(n)的)。下面介绍后缀数组的~~两种实现方法：~~倍增算法和DC3算法。

【倍增算法】

主要思想：用倍增的方法对每个字符开始的长度为 $2^k$ 的子字符串进行排序，求出排名，即rank值。k从0开始，每次加1，当 $2^k$ 大于n以后，每个字符开始的长度为 $2^k$ 的子字符串便相当于所有的后缀。并且这些子字符串都一定已经比较出大小，即rank值中没有相同的值，那么此时的rank值就是最后的结果。每一次排序都利用上次长度为 $2^{k-1}$ 的字符串的rank值，那么长度为 $2^{k}$ 的字符串就可以用两个长度为 $2^{k-1}$ 的字符串的排名作为关键字表示，然后进行基数排序，便得到了长度为 $2^{k}$ 的字符串的rank值。以字符串“aabaaaab”为例，整个过程如下图所示。x、y表示长度为 $2^{k}$ 的字符串的两个关键字。

【具体实现】

int wa[maxn],wb[maxn],wv[maxn],ws[maxn];int cmp(int *r,int a,int b,int l){return r[a]==r[b]&&r[a+l]==r[b+l];}void da(int *r,int *sa,int n,int m) //n为字符串的长度加1，因为在使用倍增算法前在原字符串后面加一个0{    int i,j,p,*x=wa,*y=wb,*t;    for(i=0;i
   
    =0;i--) sa[--ws[x[i]]]=i;    for(j=1,p=1;p
    
     =j) y[p++]=sa[i]-j;        for(i=0;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

待排序的字符串放在r数组中，从r[0]到r[n-1]，长度为n，且最大值小于m。为了函数操作方便，约定除r[n-1]外所有的r[i]都大于0，r[n-1]=0.函数结束后，结果放在sa数组中，从sa[0]到sa[n-1]。

【步骤解析】

首先，要对长度为1的字符串进行排序。一般来说，在字符串题中r的最大值不会很大，所以这里使用了基数排序。如果r的最大值很大，那么把这段代码改成快速排序。

for(i=0;i
   
    =0;i--) sa[--ws[x[i]]]=i;

以"aabaaaab"为例，执行完第一次排序后的后缀数组为sa[7]=7，sa[5]=6，sa[4]=5，sa[3]=4，sa[2]=3，sa[6]=2，sa[1]=1，sa[0]=0。

a a b a a a a b

0 1 6 2 3 4 5 7

接下来进行若干次基数排序，在实现的时候，这里有一个小优化。基数排序要分两次，第一次是对第二关键字排序，第二次是对第一关键字排序。对第二关键字排序的结果实际上可以利用上一次求得的sa数组直接算出，没有必要再算一次。

for(p=0,i=n-j;i
   
    =j) y[p++]=sa[i]-j;

其中变量j是当前字符串的长度，数组y保存的是对第二关键字排序的结果。

然后对第一关键字进行排序：

for(i=0;i
   
    =0;i--) sa[--ws[wv[i]]]=y[i];

这样便求出了新的 sa 值。在求出 sa 后，下一步是计算 rank 值。这里要注意的是，可能有多个字符串的 rank 值是相同的，所以必须比较两个字符串是否完全相同，y 数组的值已经没有必要保存，为了节省空间，这里用 y 数组保存 rank 值。这里又有一个小优化，将 x 和 y 定义为指针类型，复制整个数组的操作可以用交换指针的值代替，不必将数组中的值一个一个的复制。

for(t=x,x=y,y=t,p=1,x[sa[0]]=0,i=1;i

其中cmp函数的代码为：

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

这里可以看到规定 r[n-1]=0 的好处，如果 r[a]=r[b]，说明以 r[a]或 r[b] 开头的长度为l的字符串肯定不包括字符r[n-1]，所以调用变量 r[a+l]和 r[b+l] 不会导致数组下标越界，这样就不需要做特殊判断。执行完上面的代码后，rank 值保存在 x 数组中，而变量 p 的结果实际上就是不同的字符串的个数。这里可以加一个小优化，如果 p 等于 n，那么函数可以结束。因为在当前长度的字符串中，已经没有相同的字符串，接下来的排序不会改变 rank 值。例如以“aabaaaab”为例时的第四次排序，实际上是没有必要的。对上面的两段代码，循环的初始赋值和终止条件可以这样写：

for(j=1,p=1;p

在第一次排序以后，rank 数组中的最大值小于 p，所以让 m=p。

【复杂度分析】

每次基数排序的时间复杂度为O(n)，排序的次数决定于最长公共子串的长度，最坏情况下，排序次数为logn次，所以总的复杂度为O(nlogn)。

【后缀数组的应用】

【最长公共前缀】

这里先介绍后缀数组的一些性质。

height 数组：定义 height[i]=suffix(sa[i-1]) 和 suffix(sa[i]) 的最长公共前缀，也就是排名相邻的两个后缀的最长公共前缀。那么对于 j 和 k，不妨设 rank[j]<rank[k]，则有以下性质：

suffix(j) 和 suffix(k) 的最长公共前缀为 height[rank[j]+1], height[rank[j]+2], … ,height[rank[k]]中的最小值。

例如，字符串为"aabaaaab"，求后缀"abaaaab"和后缀"aaab"的最长公共前缀，如下图所示：

那么应该如何高效的求出 height 值呢？

如果按 height[2]，height[3]，……，height[n]的顺序计算，最坏情况下时间复杂度为 $O(n^2)$ 。这样做并没有利用字符串的性质。定义 h[i]=height[rank[i]]，也就是 suffix(i)和在它前一名的后缀的最长公共前缀。那么h 数组有性质：h[i]≥h[i-1]-1 。

证明：

设 suffix(k) 是排在 suffix(i-1) 前一名的后缀，则它们的最长公共前缀是 h[i-1]。那么 suffix(k+1)将排在 suffix(i)的前面（这里要求 h[i-1]>1，如果 h[i-1]≤1，原式显然成立）并且 suffix(k+1)和 suffix(i)的最长公共前缀是 h[i-1]-1，所以 suffix(i)和在它前一名的后缀的最长公共前缀至少是 h[i-1]-1。按照 h[1],h[2],……,h[n]的顺序计算,并利用 h 数组的性质，时间复杂度可以降为 O(n)。

【具体实现】

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

【例题】最长公共前缀

例1：给定一个字符串，询问某两个后缀的最长公共前缀。

算法分析：

按照上面所说的做法，求两个后缀的最长公共前缀可以转化为求某个区间上的最小值。对于这个 RMQ 问题可以用 O(nlogn)的时间先预处理，以后每次回答询问的时间为 O(1)。所以对于本问题，预处理时间为 O(nlogn)，每次回答询问的时间为 O(1)。如果 RMQ 问题用 O(n)的时间预处理，那么本问题预处理的时间可以做到 O(n)。

【实现】

#include 
   
    #include 
    
     #include 
     
      using namespace std;const int maxn=1e5+5;char s[maxn];int sa[maxn],wa[maxn],wb[maxn],wv[maxn],c[maxn];int Rank[maxn],height[maxn],dp[maxn][20];int cmp(int *r,int a,int b,int l){return r[a]==r[b]&&r[a+l]==r[b+l];}void da(int n,int m){    int i,j,p,*x=wa,*y=wb,*t;    for(i=0;i
      
       =0;i--) sa[--c[x[i]]]=i;    for(j=1,p=1;p
       
        =j) y[p++]=sa[i]-j;        for(i=0;i
        
         =0;i--) sa[--c[wv[i]]]=y[i]; for(t=x,x=y,y=t,p=1,x[sa[0]]=0,i=1;i
         
          = h[i-1]+1这个性质，先求出前面的后面的就可以由前面推出 j=sa[Rank[i]-1]; while(s[i+k]==s[j+k]) k++; height[Rank[i]]=k; }}void ST_build(int n){ for(int i=0;i

例2：给定一些字符串，每次询问求出两个字符串的最长公共前缀的长度。

例题：UVA 12338 - Anti-Rhyme Pairs

【单个字符串的相关问题】

【可重叠最长重复子串】

【不可重叠最长重复子串】

例题：pku1743

【可重叠的k次最长重复子串】

例题：pku3261

【不相同的子串的个数】

例题：spoj694,spoj705

【最长回文子串】

ural1297

【连续重复子串】

例题：pku2406

【重复次数最多的连续重复子串】

例题：spoj687,pku3693

【两个字符串的相关问题】

【最长公共子串】

例题：pku2774,ural1517

【长度不小于 k 的公共子串的个数】

例题：pku3415

【多个字符串的相关问题】

【不小于k个字符中的最长子串】

例题：pku3294

【每个字符串至少出现两次且不重叠的最长子串】

例题：spoj220

【出现或反转后出现在每个字符串中的最长子串】

例题：PKU3294

【完整代码】

#define maxn 1000001int wa[maxn],wb[maxn],wv[maxn],ws[maxn];int cmp(int *r,int a,int b,int l){return r[a]==r[b]&&r[a+l]==r[b+l];}void da(int *r,int *sa,int n,int m){     int i,j,p,*x=wa,*y=wb,*t;     for(i=0;i
   
    =0;i--) sa[--ws[x[i]]]=i;     for(j=1,p=1;p
    
     =j) y[p++]=sa[i]-j;       for(i=0;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
      
       b) {t=a;a=b;b=t;}    return(height[askRMQ(a+1,b)]);}