后缀数组

来源：互联网发布：unity3d 上海外包编辑：程序博客网时间：2024/05/18 01:52

后缀数组是用来替代后缀树的一种比较巧妙的数据结构，它比后缀树容易编程实现，能够实现后缀树的很多功能而时间复杂度也不太逊色，并且，它比后缀树所占用的空间小很多。

详细参考百度百科的介绍：http://baike.baidu.com/view/1240197.htm

在概念方面，百科里面介绍的比较清楚，这里就不再重复了。但是里面在具体思路及编程实现方面介绍的有点混乱（个人觉得），因而这里根据自己的理解再做一些补充，以便更好地了解这种数据结构以及用它来解决字符串中重复子串的问题。

字符串：给定一个字符串

bacac

后缀：它的后缀就是从该串的每个位置到结尾的子串。因此对任意字符串s，n=strlen(s)，共有从0到n-1开始的到结尾的不同子串，即共有n个后缀。

开始位置子串

0 bacac

1 acac

2 cac

3 ac

4 c

比较：字符串间的大小比较遵循标准函数strcmp的大小规律按照字典序进行比较。

后缀数组：对n个后缀进行排序，然后将排好序的后缀按顺序记录该后缀的开始位置存入到数组SA[]中就得到了后缀数组

排序序号开始位置子串

0 3 ac

1 1 acac

2 0 bacac

3 4 c

4 2 cac

SA[0]=3;SA[1]=1;SA[2]=0;SA[3]=4;SA[4]=2

名次数组：保存从当前位置开始的后缀（子串）在所有后缀中的排序序号

开始位置排序序号子串

0 2 bacac

1 1 acac

2 4 cac

3 0 ac

4 3 c

rank[0]=2;Srank[1]=1;rank[2]=4;rank[3]=0;rank[4]=3

上面已经用一个简单的例子介绍了后缀数据的基本概念和数据结构。我们可以看出，要得到后缀数组，就需要对后缀进行排序。而对n个后缀进行排序，若用快速排序方法需要O(nlogn)的复杂度，加上串的大小比较的时间复杂度O(n)，因而整个的时间复杂度会达到O(n*n*logn)。这样的时间效率显然不高。

对字符串比较，我们可以先比较串的第一个位置，然后第二个位置，第三个，……这样的思路显然比较吻合基数排序（http://baike.baidu.com/view/1170573.htm）的思想。此外，在这里，不同后缀之间显然还存在一定的联系。因而可以用倍增算法在基于基数排序的思想上来进行排序以使时间复杂度降到O(nlogn)。下面直接以代码来介绍该算法的思路。其中r为输入串数组，sa为后缀数组，m为串中元素的上限值（如ASCII字符的话128），函数da中前四个for循环对首个元素进行了基数排序。x 数组保存的值相当于是 rank 值。下面的操作只是用 x 数组来比较字符的大小。

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){
      int i,j,p,*x=wa,*y=wb,*t;
      for (i=0;i<m;i++) w[i]=0;
      for (i=0;i<n;i++) w[x[i]=r[i]]++;
      for (i=1;i<m;i++) w[i]+=w[i-1];
      for (i=n-1;i>=0;i--) sa[--w[x[i]]]=i;
      for (p=1,j=1;p<n;m=p,j*=2){

          //-------------对第二关键字排序从小到大---------------------
          for (p=0,i=n-j;i<n;i++) y[p++]=i;  //第二关键字为空，前j小即为这些串
          for (i=0;i<n;i++) if (sa[i]>=j) y[p++]=sa[i]-j; //利用前一次的排序结果copy过来，只不过前面已经填了j个元素，则剩下的再依次填

          //-------------根据第二关键字序对第一关键字排序---------------------
          for (i=0;i<m;i++) w[i]=0;
          for (i=0;i<n;i++) w[wv[i]=x[y[i]]]++;
          for (i=1;i<m;i++) w[i]+=w[i-1];
          for (i=n-1;i>=0;i--) sa[--w[wv[i]]]=y[i];

          //----------------计算rank值，p为不同的字符串个数------------------
          for (t=x,x=y,y=t,x[sa[0]]=0,p=1,i=1;i<n;i++)
          x[sa[i]]=cmp(y,sa[i-1],sa[i],j)?p-1:p++;
      }
      return;
}

************************************************************

--------------后缀数组应用：寻找重复系列--------------------

************************************************************

最长公共前缀（LCP）：两个字符串从开始位置出发，持续比较其对应的字符直到有不同或字符串结束，则相同的字符组成的串就是公共前缀，其字符数就是前缀数。例如，LCP("acvec","acefg")=2 ["ac"]。

对于整数i，j，从SA(i)位置开始的串和从SA(j)位置开始的串，它们的公共前缀满足：

设i<j，则LCP(i,j)=min{LCP(k-1,k)|i+1≤k≤j} （LCP Theorem）

这是因为从SA(x)开始的串是严格按照从小到大排序的（http://baike.baidu.com/view/1240197.htm#6）

根据这个性质，我们可以利用RMQ算法得到任意i,j间的LCP值。

1. 这样我们第一步先要求出相邻的串的最长公共前缀，代码如下：

void cal(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)
         for (k?k--:0,j=sa[rank[i]-1];r[i+k]==r[j+k];k++);
     return;
}

2. 第二步就是利用RMQ算法来得到从i位置开始的2^j个连续大小的串的LCP

int nmin(int a,int b){
    return a<b?a:b;
}
void rmq(int n){
     int i,j;
     for (i=1;i<=n;i++)
       f[i][0]=height[i];
     for (j=1;j<20;j++)
         for (i=1;i+(1<<j)-1<=n;i++)
         f[i][j]=nmin(f[i][j-1],f[i+(1<<j-1)][j-1]);
     return;
}

3. 计算任意两数间的LCP

int lcp(int a,int b){
    int x=rank[a],y=rank[b];
    if (x>y) {int t=x; x=y; y=t;}
    x++;
    int t=ft[y-x+1];
    return nmin(f[x][t],f[y-(1<<t)+1][t]);
}

最后求最多的重复串的办法为枚举：

for (l=1;l<=n/2;l++)
         for (i=0;i+l<n;i+=l){
             if(s[i] != s[i+l])continue;
             k=lcp(i,i+l);
             r=k/l+1;  //重复的次数
             t=i-(l-k%l);
             for(j = i-1;j > i-l && j >= 0;j --){
                     if(s[j] != s[j+l])break;
                     if(rank[j] < rank[ansp0])
                         ansp0 = j;
                     if(j == t){r ++;ansp0 = t;}  //往回数当到达这个位置时重复次数加1
             }

}

一些例子解释参考：

http://blog.csdn.net/kongming_acm/article/details/6232439

http://www.cnblogs.com/XBWer/archive/2012/05/30/2524987.html

百度之星2012初赛H题

#include <iostream>
#include <cstring>
#include <cmath>
#include <cstdio>
using namespace std;
const int maxn=100010;
const int inf=5000000;
char s[maxn];
int w[maxn],wa[maxn],wb[maxn],wv[maxn];
int sa[maxn],rank[maxn],height[maxn];
int a[maxn],f[maxn][20],n,ft[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<m;i++) w[i]=0;
     for (i=0;i<n;i++) w[x[i]=r[i]]++;
     for (i=1;i<m;i++) w[i]+=w[i-1];
     for (i=n-1;i>=0;i--) sa[--w[x[i]]]=i;
     for (p=1,j=1;p<n;m=p,j*=2){
         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;
         for (i=0;i<m;i++) w[i]=0;
         for (i=0;i<n;i++) w[wv[i]=x[y[i]]]++;
         for (i=1;i<m;i++) w[i]+=w[i-1];
         for (i=n-1;i>=0;i--) sa[--w[wv[i]]]=y[i];
         for (t=x,x=y,y=t,x[sa[0]]=0,p=1,i=1;i<n;i++)
         x[sa[i]]=cmp(y,sa[i-1],sa[i],j)?p-1:p++;
     }
     return;
}
void cal(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)
         for (k?k--:0,j=sa[rank[i]-1];r[i+k]==r[j+k];k++);
     return;
}
int nmin(int a,int b){
    return a<b?a:b;
}
void rmq(int n){
     int i,j;
     for (i=1;i<=n;i++)
       f[i][0]=height[i];
     for (j=1;j<20;j++)
         for (i=1;i+(1<<j)-1<=n;i++)
         f[i][j]=nmin(f[i][j-1],f[i+(1<<j-1)][j-1]);
     return;
}
int lcp(int a,int b){
    int x=rank[a],y=rank[b];
    if (x>y) {int t=x; x=y; y=t;}
    x++;
    int t=ft[y-x+1];
    return nmin(f[x][t],f[y-(1<<t)+1][t]);
}
void solve(int *x){
     int l,i,j,k,max=1,r=0,t;
     int ansp=n-1,ansl=1,ansp0;
     for (l=1;l<=n/2;l++)
         for (i=0;i+l<n;i+=l){
             if(s[i] != s[i+l])continue;
             k=lcp(i,i+l);
             ansp0 = i;
             r=k/l+1;
             t=i-(l-k%l);
             for(j = i-1;j > i-l && j >= 0;j --){
                     if(s[j] != s[j+l])break;
                     if(rank[j] < rank[ansp0])
                         ansp0 = j;
                     if(j == t){r ++;ansp0 = t;}
             }
            if (r>max){
               max=r;
               ansp = ansp0;
               ansl = l;
            }
            else if(r == max && rank[ansp] > rank[ansp0]){
                    ansp = ansp0;
                    ansl = l;
            }
     }
     for(i=ansp;i<ansp+ansl*max;i++)printf("%c",s[i]);printf("\n");
     return;
}
int main(){
    int i,testcase=0;
    for (i=0;i<maxn;i++) ft[i]=int(double(log(i))/log(2.00));
    while (scanf("%s",&s),s[0] != '#'){
          n = strlen(s);
          for (i=0;i<n;i++)
              a[i]=s[i];
          a[n]=0;
          da(a,sa,n+1,128);
          cal(a,sa,n);
          rmq(n);
          printf("Case %d: ", ++testcase);
          solve(a);
    }
    return 0;
}