suffix array - 使用倍增算法

演算法筆記-String Matching

suffix array ( 後綴陣列 )

使用倍增算法(Prefix Doubling)构造后缀数组

如果采用对每个后缀排序的方法来生成后缀数组，即使采用快速排序，由于每次比较的对象是字符串（设输入字符串的长度为n），因此每一个比较操作的复杂度不再是常数，而是O(n)，快速排序本身的平均情况下时间复杂度为O(nlgn)，因此总的时间复杂度是O(n^2*lgn)，如果考虑到采用快速排序最坏情况下复杂度为O(n^2)，那么最坏时总的复杂度为O(n^3)，显然在输入串长度很大时这种做法不可取。

Prefix Doubling算法（倍增法）是构造后缀数组一个比较实用的算法。其基本思想是先计算出每个后缀的k-前缀的rank值，然后在此基础上计算每个后缀的2k-前缀rank值，k从1开始。直到每个后缀都排出先后顺序为止（任何两个后缀都不会相等，换句话说，每个后缀终究都能排出先后顺序）。在处理2k-前缀时，只需要使用基数排序（radix sort）算法，先后对两位数字排序（可以采用计数排序算法(counting sort)对每一位数字排序）。在最坏情况下，需要做lgn次基数排序，每一次基数排序的操作次数为2*O(n)，因此它的时间复杂度是O(nlgn)。倍增法虽然没有达到像DC3算法的线性复杂度，但是它的优点是实现比较简单，因此在实践中常被采用。

在以下实现中，当k=1时，由于只需要对输入串的每个字符排序，因此在采用基数排序时，只有一位数字需要排序。当k>1时，需要对两位数字排序（为考虑通用性，代码中Tuple.digits数组长度可以取>=1的任何整数，而不限于两位数字的基数排序）。如果rank数组中某个后缀具有最大rank值，而且该rank值等于输入串的长度，这时说明每一个后缀都排出了次序，因而可以提前终止程序。另外，假设字母表中最大的字符为MAX_CHAR。

注：rank数组中的rank值对应SA数组的下标+1。例如T=mississippi#,那么SA=(11,10,7,4,1,0,9,8,6,3,5,2),这里第7个后缀（即ippi$）在SA数组中的下标是2，因此Rank[7]=3。

//reference: ACM/ICPC代碼庫  http://ir.hit.edu.cn/~lfwang/docs/MyLib(For%20ACM).pdf
#include<stdio.h>
#include<string>
#define REP(i, b, n) for (int i = b; i < n; i++)
#define rep(i, n) REP(i, 0, n)
#define DBG 0
#define N 512
using namespace std;
//////////////////////////////////////////////////////
// 后綴數組 O(N * log N), prefix doubling algorithm
// INIT: n = strlen(s) + 1;    (利用\0作為分隔符號)
// CALL: makesa(); lcp();
//////////////////////////////////////////////////////

//////////////////////////////////////////////////////
//sa: suffix array
//height: LCP value
//rank: sa之反函數，i,e., sa[rank[i]]= i
//tmp: 暫存最新之sa
//top: top[i]在迴圈中表示目前計算中排序第i名之字串前方(第1~i-1名)共有幾個字串
//////////////////////////////////////////////////////
int n,sa[N], height[N], rank[N], tmp[N], top[N];
char s[N];      // N > 256
void makesa(){  // O(N * log N)
    int i, j, len, na;
    na = (n < 256 ? 256 : n);
    memset(top, 0, na * sizeof(int));
    for (i = 0; i < n ; i++) top[ rank[i] = s[i]]++;
    for (i = 1; i < na; i++) top[i] += top[i - 1];          //counting sort
    for (i = 0; i < n ; i++) sa[ --top[ rank[i] ] ] = i;
    for (len = 1; len < n; len <<= 1) {             //len每次乘2
        if(DBG){rep(k,n)printf("%d ",sa[k]);printf("\n");}
        for (i = 0; i < n; i++) {                   //計算最新之sa, radix sort
            j = sa[i] - len;
            if(DBG)printf("j= %d - %d = %d(%d)\n",sa[i],len,j,j<0?j+n:j);
            if (j < 0) j += n;
            tmp[ top[ rank[j] ]++ ] = j;
            if(DBG)printf("tmp[%d]=%d\n",top[ rank[j] ]-1,j);
        }
        //此處sa被拿來暫存新的rank
        sa[ tmp[ top[0] = 0 ] ] = j = 0;
        for (i = 1; i < n; i++) {                   //計算top和rank
            if (rank[ tmp[i] ] != rank[ tmp[i-1] ] || rank[ tmp[i]+len ]!=rank[ tmp[i-1]+len ])
                top[++j] = i;
            sa[ tmp[i] ] = j;
        }
        if(DBG)printf("j= %d\n",j);
        memcpy(rank, sa , n * sizeof(int));
        memcpy(sa , tmp, n * sizeof(int));
        if (j >= n - 1) break;
    }
}
//lcp: 計算從每個字串之LCP (from s[0] to s[n-1])
//note: rank[0]>0, 因sa[0]=n (\0)
void lcp(){ // O(4 * N)
    int i, j, k;
    for (j = rank[height[i=k=0]=0]; i < n - 1; i++, k++){
        while (k >= 0 && s[i] != s[ sa[j-1] + k ]){
            height[j] = (k--), j = rank[ sa[j] + 1 ];
        }
    }
}
int main(){
    bool ll=0;
    char arr[5][100]={"ttttttt", "mississippi#", "abcdefghijklmmnopqrstuvwxyz#",
        "yabbadabbado#", "DFDLKJLJldfasdlfjasdfkldjasfldafjdajfdsfjalkdsfaewefsdafdsfa#"
        };
    rep(i,5){
        strcpy(s,arr[i]);
        n=strlen(s)+1;              //+1!
        if(ll)printf("********************************\n");ll=1;
        printf("Text: %s\n",s);
        makesa();
        lcp();
        printf("suffix array: \n");
        rep(i,n)printf(" %d",sa[i]);printf("\n");
        printf("rank array: \n");
        rep(i,n)printf(" %d",rank[i]+1);printf("\n");
        //printf("lcp: ");
        //rep(i,n)printf("%d ",height[i]);printf("\n");
    }
}