Algorithm::Suffix Array 後綴數組

這次去了 IOICamp,從第一次看到到現在過了這麼久,終於學會 Suffix Array 了QAQ。
大概把理解的過程寫下來。

在處理字串問題的時候,我們常常需要處理一些跟子字串啦、匹配字串等等的問題,當我們需要拿一些東西對一個大字串(約10^5)做一些事情的時候,Suffix Array也就是後綴數組也許就可以幫我們做到。

所謂的後綴數組就是,我們把一個字串的所有的後綴拿去排序,然後紀錄一些數值。比如字串”ababa”變成後綴數組就是:

Rank SA Hei Suffix
0 0 0 ababa
1 2 3 aba
2 4 1 a
3 1 0 baba
4 3 2 ba

實際上我們只會做出前面的 Suffix Array(以下簡稱 SA(x) ),而不會真的把子字串做出來拿去排序。
做成這樣的好處就是我們可以利用它進行二分搜等等的操作,因為每個後綴的前綴就是一個子字串,我們就像是可以對所有的子字串進行二分搜,當然還要搭配一些東西來輔助,像是 Hei(x),代表每個子字串和它排名前一名的人,前面有多少字元是相等的,搭配這個資訊就可以作到一些事情,當然也可以真對不同題目去找不同的值。

要如何找到一個 SA(x) 呢?
我們要利用一種特別的找法,也就是倍增法。
當然可以直接用 std::sort 來排序,但是字串的比較是 O(n) 有點糟糕,但我們發現因為其實後綴原本是同一個字串,比較的時候我們可以省下一些時間,利用其他人比較後的結果來排序。

首先我們先定義兩個必要的東西, 一個就是 y = SA(x) 代表第 x 名的後綴實際上是從 y 開始的後綴; 另外一個就是 x = RA(y) 代表 y 開始的後綴是第 x 名。
很容易可以發現他們是反函數關係(x = RA(S(x)) )。
然後我們利用 radix sort 的概念,也就是排序字串的時候先比較第一個字元,照著排好序之後,分成很多組,針對第一個字元一樣的人,再對第二個字元排序…。
這時我們可以發現一個性質就是,一個 i 到 i+2K 的子字串和 j 到 j+2K 的子字串,想知道它的大小關係,可以透過比較到 i+K 和 j+K 的子字串,如果一樣,再去比較 i+K 到 i+2K 和 j+K 到 j+2K 的子字串,反正就是拆成兩塊啦!
我們又知道我們後面的字串已經排好序了,我們就可以利用他們一層一層推出來!

看看這段程式碼:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
void build_suffix_array(int*SA, int*RA, char*S, bool is_dec=true){
static int RA2[N];
int n = strlen(S);
copy_n(S,n,RA);
int n_of_diff_rank = 256; // char has 256 different value.
for(int k=1; k<n; k*=2){
// Comput SA_{2k} with RA_{k}
// sort by rank-k ( ignore the string which len < 1+k )
if(is_dec){ // "aa" < "a"
iota(SA, SA+n, 0); // init SA to [0,1,2...]
radix_sort(SA, SA+n-k, RA+k, n_of_diff_rank);
}else{ // "a" < "aa"
iota(SA,SA+k,n-k);
iota(SA+k,SA+n,0); // init SA to [k,k+1...n-1,0,1,...]
radix_sort(SA+k, SA+n, RA+k, n_of_diff_rank);
}
// sort by rank-0 ( include the string which len < 1+k )
radix_sort(SA, SA+n, RA, n_of_diff_rank);
// Comput RA_{2k} with SA_{2k}
RA2[SA[0]] = 0; n_of_diff_rank = 1;
for(int i=1; i<n; i++){
if(noteq(SA[i-1], SA[i], k, n, RA) )
n_of_diff_rank++;
RA2[SA[i]] = n_of_diff_rank-1;
}
copy_n(RA2, n, RA);
if(n_of_diff_rank == n)break; // all rank are diff => SA is fully sorted
}
}

先把字串的值直接複製到 RA 裡,當作一開始的大小關係的依據,RA 裡面就是紀錄 i 這個位置的後綴是第幾名,一開始有很多重複的值(重複的字母),經過多次的排序,我們的目標就是讓它變成 0 到 N-1 全部都不一樣的排名!利用這些排名,就可以把 SA(x) 給排序好!!

首先先照字串的第二段的值去排序 SA(x) ,如果有長度不夠的字串(也就是位置再 n-k 之後的後綴),就當作它沒有第二段,我們可以先放到旁邊不理它。要照第二段排序只要排序 i 的時候參考 RA(i+K) 就可以照第二段排序了!
因為我這裡的範例用的sort是穩定排序,可以透過把長度不夠的字串放到 SA 前面或是 SA 後面決定不夠長的位置要再 a 前面還是在 z 後面(最上面表格裡的例子就是放在後面的,ababa 然後 aba),不過基本上沒差啦,除非題目有特殊性質要用。

再來第二步就是再照字串第一段的值去排序,要作到這樣就是比較 i 的時候用 RA(i) 去比較就可以了。

到這裡我們已經把依照 K 做出的 SA-2K 做出來了(依照前 2K 個字元排序的 SA ,依此類推也有 RA-2K),再來就是用它來更新 RA 讓它也變成 RA-2K ,作法就是先把「名次」設成 1 ,然後從第0名開始看,如果這一個名次跟上一個名次是一樣的( 用上面那種切兩段比較法來比較 ),那這個人就跟它前一個人並列第 x 名,做完之後我們就得到依照 2K 的 RA-2K 了。

最後,如果名次已經跟 n 一樣多了那代表已經可以排序所有後綴了,不然就把 K 變成兩倍,如果長度超過 n 那代表也不用繼續了,如果都不是那就從第一步重來。

另外比較有趣的是這段 radix sort:

1
2
3
4
5
6
7
8
9
10
11
void radix_sort(int*bg, int*ed, int*rank, int n_of_diff_rank){
static int box[N], tmp[N];
int n = ed-bg;
fill_n(box, n_of_diff_rank, 0);
for(int i=0; i<n; i++)
box[rank[i] ]++;
partial_sum(box, box+n_of_diff_rank, box);
for(int i=n-1; i>=0; i--)
tmp[--box[rank[bg[i]] ]] = bg[i];
copy_n(tmp,n,bg);
}

先把會有的不同的值的範圍先初始化好,然後用 partial_sum 轉成前綴和!!!最後用倒序的順序把東西填回去,這樣就可以作到穩定的排序,蒸的很蚌!

最後就是怎麼做出 hei 數組,單純的做的話就是跟前一個人比較,算出有幾個一樣,當然這樣很慢。
有一個性質是,AAB 和它前面一個位置 BAAB ,你可以看到它只往前多一個 B 而已,所以你前面一名和它前面一名你可以大膽的猜,要不是他的 hei 往前多 1 , 要不然就是它跟它比的人前面多的不一樣變成 0 ,反正它最多 1 所以你的 hei 就是大於等於他的 hei -1 !!!
所以你只要照著位置順序做,每次做下個位置的時候只要從上一個 hei -1 開始找就好,這樣複雜度就是字串最長 n + 每次下一個位置長度最多減一就是好好的 O(2*n)。

1
2
3
4
5
6
7
8
9
10
void build_hei(int*SA, int*RA, char*S, int*Hei){
int k = 0, n = strlen(S);
for(int i=0;i<n;i++){
if(RA[i]){
while(S[i+k]==S[SA[RA[i]-1]+k])k++;
}
Hei[RA[i]] = k;
k = max(0,k-1); // hei(rank(i+1) ) >= hei(rand(i) )-1
}
}

最後說個簡單的例題,比如要你求一個字串的不相等的子字串有幾個?
作法就是把 SA 求出來之後做 Hei ,然後所有子字串的可能個數扣掉 Hei 的總和(剛好就是重複的子字串數)。

下面就是完整的程式碼。
實做上可以再壓一點變數之類的,不過為了寫清楚就這樣寫了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
#include <bits/stdc++.h>
using namespace std;
const int N = 1e5+1;
void debug_suffix(int *SA, int *RA, int *Hei, char *S, const char*s="suffix array:"){
int n = strlen(S);
cout<<s<<endl;
cout<<"SA Hei:"<<endl;
for(int i=0;i<n;i++)cout<<SA[i]<<' '<<Hei[i]<<' '<<setw(n)<<(S+SA[i])<<endl;
cout<<"RA:"<<endl;
for(int i=0;i<n;i++)cout<<RA[i]<<' '<<setw(n)<<(S+i)<<endl;
cout<<"-------"<<endl;
}
bool noteq(int a, int b, int k, int n, int*RA){
return RA[a] != RA[b] || a+k>=n || b+k>=n || RA[a+k] != RA[b+k];
}
void radix_sort(int*bg, int*ed, int*rank, int n_of_diff_rank){
static int box[N], tmp[N];
int n = ed-bg;
fill_n(box, n_of_diff_rank, 0);
for(int i=0; i<n; i++)
box[rank[i] ]++;
partial_sum(box, box+n_of_diff_rank, box);
for(int i=n-1; i>=0; i--)
tmp[--box[rank[bg[i]] ]] = bg[i];
copy_n(tmp,n,bg);
}
void build_suffix_array(int*SA, int*RA, char*S, bool is_dec=true){
static int RA2[N];
int n = strlen(S);
copy_n(S,n,RA);
int n_of_diff_rank = 256; // char has 256 different value.
for(int k=1; k<n; k*=2){
// Comput SA_{2k} with RA_{k}
// sort by rank-k ( ignore the string which len < 1+k )
if(is_dec){ // "aa" < "a"
iota(SA, SA+n, 0); // init SA to [0,1,2...]
radix_sort(SA, SA+n-k, RA+k, n_of_diff_rank);
}else{ // "a" < "aa"
iota(SA,SA+k,n-k);
iota(SA+k,SA+n,0); // init SA to [k,k+1...n-1,0,1,...]
radix_sort(SA+k, SA+n, RA+k, n_of_diff_rank);
}
// sort by rank-0 ( include the string which len < 1+k )
radix_sort(SA, SA+n, RA, n_of_diff_rank);
// Comput RA_{2k} with SA_{2k}
RA2[SA[0]] = 0; n_of_diff_rank = 1;
for(int i=1; i<n; i++){
if(noteq(SA[i-1], SA[i], k, n, RA) )
n_of_diff_rank++;
RA2[SA[i]] = n_of_diff_rank-1;
}
copy_n(RA2, n, RA);
if(n_of_diff_rank == n)break; // all rank are diff => SA is fully sorted
}
}
void build_hei(int*SA, int*RA, char*S, int*Hei){
int k = 0, n = strlen(S);
for(int i=0;i<n;i++){
if(RA[i]){
while(S[i+k]==S[SA[RA[i]-1]+k])k++;
}
Hei[RA[i]] = k;
k = max(0,k-1); // hei(rank(i+1) ) >= hei(rand(i) )-1
}
}
char S[N];
int SA[N], RA[N], Hei[N];
int main(){
int t;
cin>>t;
while(t--){
cin>>S;
build_suffix_array(SA,RA,S,true);
build_hei(SA,RA,S,Hei);
debug_suffix(SA,RA,Hei,S,"Suffix Array(DES)");
build_suffix_array(SA,RA,S,false);
build_hei(SA,RA,S,Hei);
debug_suffix(SA,RA,Hei,S,"Suffix Array(AES)");
}
}