演算法知識 - Suffix Array 後綴陣列

Suffix Array 介紹

對某一字串的所有後綴進行字典排序，常用於全文索引、數據壓縮算法與生物資訊學。
這裡介紹的演算法寫法時間複雜度為 \(O(n \log n)\)

Suffix Array 原理

Suffix Array 主要是用 \(sa\) and \(rk\) 這兩個陣列組合而成。
且滿足此性質 \(sa[rk[i]] = rk[sa[i]] = i \)

Suffix Array 圖示說明：

QUESTION 1: sa (suffix array) 用途是甚麼呢？

\(sa[i]\) 表示此字串所有後綴排序後第 i 大的 index。

QUESTION 2: rk (rank array) 用途是甚麼呢？

\(rk[i]\) 表示此字串所有後綴排序後，字串中 index 後綴的排名。

QUESTION 3: 甚麼是所有後綴？

從第 i 個字元開始到最後一個字元的字串，i 的範圍是 string[0] ~ string.length()
所有後綴舉例如下，舉例單字為 apple：

• apple
• pple
• ple
• le
• e

Suffix Array 實現與說明

在實現這個做法時需要用到倍增思想。

QUESTION: 倍增是甚麼？

倍增是寫程式中經常用到的一種加速程式執行效率的想法之一，通常是把一個問題拆成兩個子問題且這兩個子問題的問題邏輯相同，於是就使用一次寫法來加速效率。

倍增經典範例：快速冪 OI wiki

回歸焦點

OK，解釋完倍增的想法後就要繼續探討 Suffix Array，我們先對字串中每個字元(長度為\(1\))進行排序，再來我們對字串中以每個字元開始長度為 2 的倍數開始進行排序，這裡我們假設字串除了原先的字元之外，之後的 index 全部都是空字元，方便我們對演算法的編寫。記憶體很大，我們從沒在怕的，科技真讚

我們這裡字串的第一個 index 由 1 開始，方便我們程式編寫，建議使用者也這樣做，只需要在 string A 前面加一空格。

倍增排序圖示說明：

其中如果黑線只有一條而沒有類似於勾勾的黑線(通常都在字串右方)，可以直接把她想像成此字串後面其實還有字元，但都是「空白字元」，如果這樣假設那讀者勢必會更快理解XD。我自己把她這樣解釋後才能理解，我腦袋真笨QQQQ，雜牌軍日常呀。

基數排序

如果不懂基數排序可以先看演算法知識 - Radix Sort 基數排序

由於字串比較的時間複雜度是 \(O(n)\)，跟一般數字比較不一樣，因此如果這裡直接使用 STL 函數庫中的 std::sort 會讓此演算法時間複雜度來到 \(O(n \log^2 n)\)，於是這裡推薦使用非比較型排序，基數排序，透過這種排序方式來躲開字串的時間比較複雜度才可以讓 Suffix Array 的時間複雜度降到 \(O(n \log n)\)

是不是已經被冗長的文字敘述搞到頭暈了？沒關係，這裡準備了程式碼已經逐行註解，希望可以讓各位讀者都能比我更快速的讀懂 Suffix Array

int sa[N] , rk[N<<1] , oldrk[N<<1] , id[N] , cnt[N] ;
// id = 舊的 sa 排名 , cnt 在 radix sort 時幫助
int n , m , maxn , lenA , lenB , flag =0 ;
// n 字串長度 , m 文字的最大長度
string strA ;
// 要進行 suffix array 的字串，請記住這裡的字串第一個 index 為 1

void build_sa(){ //build suffix array
    int i , m , p , w ; 
    // 在這邊把常在迴圈使用的變數名稱拉出來宣告，以避免不斷宣告造成的效率浪費
    n = strA.length()-1 ; //減去 string[0] 的空格
    m = max(n , 300 ); // 由於 ascii 上限是 255，我們這邊直接開 300
    // 也省去把字元 hash 的麻煩
    memset(cnt,0,sizeof(cnt)); //重設 cnt 陣列以免多次使用時與上次使用的值混亂在其中
    memset(rk,0,sizeof(rk)); //重設 rk 陣列以免多次使用時與上次使用的值混亂在其中

    // **** radix sort 排序開始 ****
    for(i = 1 ; i <= n ; i++) ++cnt[rk[i] = (int)strA[i]] ;
    // rk[i] = (int) strA[i] 每一個字元在字串中排名
    //++cnt[rk[i] = (int)strA[i]] radix sort 分類
    for(i = 1 ; i <= m ; i++) cnt[i] += cnt[i-1] ; 
    //將 cnt 遞增排序方便之後找出數值排在第幾位
    for(i = n ; i >= 1 ; i--) sa[cnt[rk[i]]--] = i ;
    // cnt[rk[i]] 先找出字元在字串中排名再透過 cnt 找出他應該排在第幾位
    // 隨後進行 cnt[rk[i]]--，方便下一個值得排序位置。
    // **** radix sort 排序結束 ****

    for(w = 1 ; w < n ; w <<= 1){ // 倍增思想開始
        memset(cnt,0,sizeof(cnt)); //第一次 radix sort 開始，排序關鍵為後者
        // 最好範例：倍增排序圖中黑色粗線類似於勾勾的斜線，下方解釋 A 圖中
        for(i = 1 ; i <= n ; i++) id[i] = sa[i] ;
        // id 用來記錄現在排序的順序，配合這次的 radix sort 再進行更動
        for(i = 1 ; i <= n ; i++) ++cnt[rk[id[i]+w]] ;
        for(i = 1 ; i <= m ; i++) cnt[i] += cnt[i-1] ;
        for(i = n ; i >= 1 ; i--) sa[cnt[rk[id[i]+w]]--] = id[i] ;
        // radix sort 與一開始相同，只是多增加 w，因為排序關鍵為後者
        // 第一次 radix sort 結束

        // 第二次 radix sort 開始，排序關鍵為前者
        memset(cnt, 0 , sizeof(cnt));
        for(i = 1 ; i <= n ; i++) id[i] = sa[i] ;
        for(i = 1 ; i <= n ; i++) ++cnt[rk[id[i]]] ;
        for(i = 1 ; i <= m ; i++) cnt[i] += cnt[i-1] ;
        for(i = n ; i >= 1 ; i--) sa[cnt[rk[id[i]]]--] = id[i] ;
        // 第二次 radix sort 結束

        // 根據 radix sort 在排列一次 rank 
        memcpy(oldrk , rk , sizeof(rk)); //複製 rank 來幫助 rank 間的交換
        for(p = 0 , i = 1 ; i <= n ; i++){ // p 為 rank 的等級，i 則是資料
            if(oldrk[sa[i]] == oldrk[sa[i-1]] && 
               oldrk[sa[i] + w] == oldrk[sa[i-1] + w])
               //如果與前一個 rank 值是相同則理應現在應該也要相同
               //最好範例：倍增排序圖示說明第一次排序的 rank[4~7], index 從 1 開始，
               //下方解釋 B 圖中
                rk[sa[i]] = p ; //rank 字典排序不變
            else
                rk[sa[i]] = ++p ; //rank 字典排序增加
        }
    }

//debug 輸出測試，以驗證是否正確
//    cout << "Suffix Array is:\n" ;
//    for(int i = 1 ; i <= n ; i++){
//        cout << i << ' ' << strA.substr(sa[i]) << ' ' <<sa[i] << '\n' ;
//    }
}

• 解釋 A 圖
  
• 解釋 B 圖
  

最長共同前綴 Longest Common Prefix Array (LCP)

定義 height = Longest Common Prefix Array
道理其實相當簡單，我們的 Suffix Array 是字典排序，於是我們可以推出一公式 \(height[i] = lcp(sa[i],sa[i-1]\)，也就是讓第 i 名的後綴去跟前一名後綴算出最長共同前綴。

\(height[1] = 0\)，由於我們的 string index 是從 1 開始，所以 1 只能夠跟 0 比無意義因而設成 0。

因為 height 他每個 index 都是獨立並沒有相關性，我們的比較方式是根據字串中的每個後綴由字串 index 順序，去找出他的 sa -1 位置去找出 lcp 長度，下一次則是找 index +1 的後綴與他的 sa-1 位置去找出 lcp 長度，由於比較的字串只是刪除上一個前綴的第一個字元，於是 lcp 長度最差則是 lcp -1(前提為 lcp > 0)，所以可以直接從 index + lcp 的長度直接比較，可以減少重複比較，以達到降低時間複雜度。

時間複雜度是 \(O(n)\)。

如果還是聽不太懂，感覺有點紙上談兵，就讓我用程式碼來解釋吧！

void build_lcp(){
    int lcp[N] = {} ;
    int max_lcp = 0 ; // max_lcp 最大長度
    //k 為現在 i 名的後綴與前一名後綴長度算出的最長共同前綴
    for(int i = 1 , k = 0 ; i <= n ; i++){
        if(k) k-- ; 
        // 由於我們下次是把字串 index 在往後一個 index，所以 lcp 最差的情況則會是 x-1
        // 因為其實要比較的字串只是刪除上一個後綴的第一個字元
        
        while(strA[i+k] == strA[sa[rk[i]-1]+k]) ++k ;
        //比較順序較為特殊，因為他每個陣列都可以分開進行討論，
        //於是我們比較字串的 index 開始的後綴與他的 sa 前一項找出 lcp
        // strA[i] = 字串的 index 開始的後綴
        // strA[sa[rk[i]-1]] = 字串的 index 開始的後綴的 sa 前一項
        // while 裡面的 +k 則是減少重複比較，如果上次的 lcp 已經找出長度為 x 的 lcp

        lcp[rk[i]] = k ;
    }

    for(int i = 1 ; i <= n ; i++){
        if((sa[i] < lenA && sa[i-1] < n+1 && sa[i-1] > lenA ) ||
            (sa[i] > lenA && sa[i-1] < n+1 && sa[i-1] < lenA))
            max_lcp = max(max_lcp , lcp[i]);
    }
}

以上，Suffix Array 與 LCP 到此介紹結束，謝謝各位不嫌棄我的文筆閱讀到此。

Suffix Array 應用

我在學習此演算法的過程中如果有題目讓我應用到此演算法我則收錄在此，連結則是我的詳解如果想要去看題目的就搜尋一下吧！

最長共同前綴 Longest Common Prefix Array

• UVa 760

參考連結：

快速冪 OI wiki
UVa 760 – DNA Sequencing
常用ASCII码详细对照表 （0—255）

心得

學習演算法真的是 CP 值不高又很浪費時間的東西呢(無誤XD)，除了 ACM ICPC 之外，你很難知道甚麼時候可以用到，又或者當你想用時卻已經忘記怎麼用了，也有可能到時已經有現成套件可以用，那你為甚麼還要學習演算法呢？

因為有趣，可以增加思維，我自認我從演算法學習到了許多不同新事物，透過演算法來認識這個世界，這些演算法都是優秀的電腦科學家研究出來，透過他們的思維來引道我來看待事物，我相信，會讓我變得更加優秀。

也希望我在學習演算法的路上可以更加優秀，不被打倒，學習速度更快更容易讓我獲得新知識擴展自己腦袋的世界觀，也同時保持著舊有知識不遺忘。

最後要感謝 OI WIKI，他讓對於英文不好的我可以有更多機會去接觸演算法，將各個演算法收錄在此且無私的開放給大家學習，讓我有機會可以學習到此演算法。

也謝謝我自己的努力付出，即使我沒有機會成為一個優秀的人但或許可以成為他人的基礎吧XD。

題目程式碼

其實就是沒有附上說明的程式碼，我也花了三個小時撰寫文章，看得懂跟說得出來真的是兩回事，也許還有人覺得我說的很差，不過我已經盡力解釋。希望可以幫助到各位。

也希望大家想要學習的知識都可以快速學習而成，不要因為自身的環境不足而學習不到。

於是這裡就補述一些變數宣告與資料範圍限制與標頭檔了。

#include <iostream>
#include <bits/stdc++.h>
#define LOCAL
#define N 2000

int sa[N] , rk[N<<1] , oldrk[N<<1] , id[N] , cnt[N] ;
int n , m  ;
string strA ;

void build_sa(){
    int i , m , p , w ;
    n = strA.length()-1 ;
    m = max(n , 300 );
    memset(cnt,0,sizeof(cnt));
    memset(rk,0,sizeof(rk));
    for(i = 1 ; i <= n ; i++) ++cnt[rk[i] = (int)strA[i]] ;
    for(i = 1 ; i <= m ; i++) cnt[i] += cnt[i-1] ;
    for(i = n ; i >= 1 ; i--) sa[cnt[rk[i]]--] = i ;

    for(w = 1 ; w < n ; w <<= 1){
        memset(cnt,0,sizeof(cnt));
        for(i = 1 ; i <= n ; i++) id[i] = sa[i] ;
        for(i = 1 ; i <= n ; i++) ++cnt[rk[id[i]+w]] ;
        for(i = 1 ; i <= m ; i++) cnt[i] += cnt[i-1] ;
        for(i = n ; i >= 1 ; i--) sa[cnt[rk[id[i]+w]]--] = id[i] ;

        memset(cnt, 0 , sizeof(cnt));
        for(i = 1 ; i <= n ; i++) id[i] = sa[i] ;
        for(i = 1 ; i <= n ; i++) ++cnt[rk[id[i]]] ;
        for(i = 1 ; i <= m ; i++) cnt[i] += cnt[i-1] ;
        for(i = n ; i >= 1 ; i--) sa[cnt[rk[id[i]]]--] = id[i] ;

        memcpy(oldrk , rk , sizeof(rk));
        for(p = 0 , i = 1 ; i <= n ; i++){
            if(oldrk[sa[i]] == oldrk[sa[i-1]] &&
               oldrk[sa[i] + w] == oldrk[sa[i-1] + w])
                rk[sa[i]] = p ;
            else
                rk[sa[i]] = ++p ;
        }
    }

    //debug
//    cout << "Suffix Array is:\n" ;
//    for(int i = 1 ; i <= n ; i++){
//        cout << i << ' ' << strA.substr(sa[i]) << ' ' <<sa[i] << '\n' ;
//    }
}

void build_lcp(){
    int lcp[N] = {} ;
    int max_lcp = 0 ;
    for(int i = 1 , k = 0 ; i <= n ; i++){
        if(k) k-- ;
        while(strA[i+k] == strA[sa[rk[i]-1]+k]) ++k ;
        lcp[rk[i]] = k ;
    }
}

用紙筆去模擬 Suffix Array and LCP

因為只是自己在思考過程中隨手做的筆記，覺得不是對讀者這麼有幫助，因此放在文章最下面。

有時候用想或用看的可能不太能夠懂這演算法，那就用寫的吧！這是大衛我自己學習的手稿，紀錄者我的學習過程。