rhrnald's 코딩
Constructing suffix array in linear time 본문
Suffix Array를 $O(n)$시간만에 구성하는 방법입니다. 해당 논문의 내용을 그대로 옮겨왔습니다.
길이 n인 문자열 $T$가 주어져 있을때(1-base index), Suffix Array는 두개의 배열로 구성되어 있습니다.
첫 번째는 $A_{T}$, sort array입니다. $T$의 suffix $S_{i}$를 $T[i]$에서 시작해서 $T[n]$에서 끝나는 부분문자열로 정의할 수 있습니다. 이 때 이 $A_{T}$는 suffix들을 lexicographic 순서로 정렬 하였을때, index들의 배열입니다. $\{S_{1}, \cdots, S_{n}\}$중 lexicographic 순서로 $i$번째로 작은 suffix를 $S_{j}$라 하면, $A_{T}[i]$에는 $j$가 저장됩니다.
두 번째는 $L_{T}$, lcp array입니다. $lcp(a, b)$를 a와 b의 longest common prefix로 정의 하겠습니다. 두 문자열 모두의 prefix가 되는 문자열들 중 가장 긴 것을 의미합니다. 이 때 lcp array는 $L_{T}[i]$에 $|lcp(S_{A_{T}[i]}, S_{A_{T}[i+1]})|$이 저장 되어 있는 배열입니다. $S_{i}$들을 정렬 했을 때 인접한 두 suffix의 lcp의 길이를 의미합니다. $1 \leq i< n$에서 정의 되지만 편의 상 $L_{T}[0]=L_{T}[n]=0$로 정의 하겠습니다.
이제 다음의 class를 짜줄 것입니다. $O(n \log{n})$에 구성하는 방법이 잘 알려져있지만 여기서는 $O(n)$에 구성하는 법을 한번 알아 볼 것입니다.
class SuffixArray {
public:
int len;
vector<int> T, A, L;
void init(string &s) {
len=s.size();
T.resize(len+2);
for(int i=1; i<=len; i++) T[i]=s[i-1];
}
void SA(vector<int> &_A, vector<int> &_L);
};
$SA_{o}$와 $SA_{e}$를 먼저 정의하겠습니다. $SA_{o}=\{A_{o}, L_{o}\}$, $SA_{e}=\{A_{e}, L_{e}\}$로 이루어져 있습니다. 이 때 $A_{o}$와 $A_{e}$는 각각 $S_{odd}$와 $S_{even}$를 나누어서 정렬해준 것입니다. 따라서 $A_{o}$는 홀수 인덱스들만, $A_{e}$는 짝수 인덱스들만 포함하게 됩니다. $L_{o}$와 $L_{e}$는 나누어서 정렬했을 때 마찬가지로 인접한 두 suffix의 lcp길이를 의미합니다.
예시를 한번 보겠습니다. T=aaaabbbbaaabbbaabbb#라 하겠습니다.
홀수번째 인덱스에서 시작하는 suffix들을 lexiographic 순서로 정렬하면 다음과 같습니다.
$S_{1}$=aaaabbbbaaabbbaabbb#
$S_{9}$=aaabbbaabbb#
$S_{15}$=aabbb#
$S_{3}$=aabbbbaaabbbaabbb#
$S_{11}$=abbbaabbb#
$S_{19}$=b#
$S_{7}$=bbaaabbbaabbb#
$S_{13}$=bbaabbb#
$S_{17}$=bbb#
$S_{5}$=bbbbaaabbbaabbb#
이 suffix들을 odd suffix array로 바꾸면 다음과 같이 되는 것입니다. $L_{o}[1]$에는 $S_{1}$과 $S_{9}$의 lcp의 길이가 저장되는 식입니다.
| 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | |
| $A_{o}$ | 1 | 9 | 15 | 3 | 11 | 19 | 7 | 13 | 17 | 5 |
| $L_{o}$ | 3 | 2 | 5 | 1 | 0 | 1 | 4 | 2 | 3 | 0 |
이제 세 단계를 통해 $SA_{T}=\{A_{T}, L_{T}\}$를 구해줄 것입니다.
1. Odd suffix array, $SA_{o}$ 구성
재귀를 이용하여 구해줄 것입니다. 두개씩 묶어서 새로운 문자들에 대응 시킬 것입니다. 최대 길이/2 종류 밖에 안나오므로 radix sort를 이용하여 O(n)시간에 새로운 문자에 대응시켜줄 수 있습니다. 이후에 재귀를 통해 suffix array를 구해주고 원래 문자열에 맞게 변한시켜주면 됩니다.
인덱스는 옮겨주기만 하면 되지만 lcp는 살짝의 처리가 필요한데 마지막 겹치니 않는 쌍에서 한자리 까지는 겹칠 수 있으므로 처리해주어야 합니다.
void ConstructOddSA() {
for(int i=1; i<=olen; i++) v.push_back({{T[i+i-1], T[i+i]}, i});
RadixSort(v);
vector<int> nT; nT.resize(olen+2);
int idx=1; nT[v[0].second]=1;
for(int i=1; i<olen; i++) {
if(v[i].first!=v[i-1].first) idx++;
nT[v[i].second]=idx;
}
SuffixArray nSA;
nSA.init(olen,nT);
nSA.SA(Ao,Lo);
for(int i=1; i<=olen; i++) Ao[i]=Ao[i]+Ao[i]-1;
for(int i=1; i<=olen-1; i++) {
Lo[i]=Lo[i]+Lo[i];
if(T[Ao[i]+Lo[i]]==T[Ao[i+1]+Lo[i]]) Lo[i]++;
}
ORMQ.makeRMQ(Lo);
for(int i=1; i<=olen; i++) oloc[Ao[i]]=i; oloc[len+1]=0;
}
2. Even suffix array, $SA_{e}$ 구성
마찬가지로 재귀를 이용해도 되지 않나 싶지만 안됩니다. 길이 n인 문자열로 suffix array를 만드는데 걸리는시간을 $F(n)$이라 하면 양쪽 모두에서 재귀를 사용할 경우 $F(n)>=F(n/2)+F(n/2)+n\cdot c$와 같이 되어 $O(n \log{n})$의 시간복잡도를 갖게 됩니다. 둘 중 한쪽에서는 재귀를 사용하지 않아야 $F(n)=F(n/2)+c_{1} \cdot n+c_{2} \cdot n$과 같이 되어 $O(n)$의 시간복잡도를 갖는다 말할 수 있습니다.
따라서 이번에는 재귀를 이용하는 대신에 이미 연산을 한 odd suffix array를 이용하여 $SA_{e}$를 구성할 것입니다.
짝수 suffix들은 맨 앞자리를 제거하면 홀수 suffix가 됩니다. 홀수 suffix들의 크기 관계는 이미 구해놨으니 {맨 앞자리, 홀수 suffix} 쌍으로 보아 정렬을 해주면 짝수 suffix들도 정렬 가능해집니다.
lcp 계산은 맨 앞자리가 겹치는지 일단 처리해주고 그 다음 두 홀수 suffix의 lcp를 구해주면 됩니다.
보조정리. $x<y$에 대해 $A_{T}[x]=p$, $A_{T}[y]=q$일 때, $|lcp(S_{p}, S_{q})|=\min_{x \leq i \leq y-1} L_{T}[i]$
식은 복잡해 보이지만 내용은 간단합니다. 위의 예시로 살펴보겠습니다.
$S_{1}$과 $S_{3}$의 lcp의 길이를 구하고 싶습니다. 두 suffix는 $A_{T}$에서 각각 1번째와 4번째에 위치하고 있습니다. 그렇다면 두 suffix의 lcp의 길이는 $L_{T}[1]$~$L_{T}[3]$ 중 최소인 값을 구하면 되는 것입니다. $S_{A_{T}[x]}$와 $S_{A_{T}[y]}$의 prefix가 같다면 lexiographic 순서로 정렬해 놓은 것이므로 사이의 다른 suffix들도 같은 prefix를 갖는다는 사실과 해당 사실의 역을 이용하면 증명이 됩니다.
위의 보조정리와 함께 특정 구간의 suffix들, $lcp(S_{A_T[x]}, ... , S_{A_T[y]})$는 양 끝 두개의 suffix의 lcp인 $lcp(S_{A_T[x]}, S_{A_T[y]})$와 같다는 것을 알 수 있습니다.
보조정리. $MIN(L_{T}, x, y-1)=\min_{x \leq i \leq y-1} L_{T}[i]$는 상수 시간내에 계산 될 수 있다.
RMQ입니다. $O(n)$의 전처리 이후에 특정 구간의 최솟값을 구하는 방법을 앞의 포스팅에서 소개 시켜 드렸습니다. 이 때 인덱스 까지 같이 구해주어야 할 필요가 있습니다.
Range Minimum Query
Range minimum은 배열이 주어져 있을 때, 연속한 구간의 최솟값을 구하는 것을 의미한다. 최솟값을 여러 구간에 대해 반복적으로 구하는것이 Range Minimum Query, 줄여서 RMQ라 부른다. 다양한 방법이 존재하여 P..
rhrnald.tistory.com
따라서 한번의 $O(n)$ 전처리 이후에는 매번 상수 시간에 $L_{e}[i]$를 계산할 수 있게 되고 따라서 총 $O(n)$의 시간에 $A_{e}$와 $L_{e}$를 구할 수 있게 됩니다. 코드는 다음과 같습니다.
class RMQ {
public:
vector<vector<pii>> T;
void makeRMQ(vector<int> &L) {
T.clear();
T.push_back(vector<pii>());
for(int i=0; i<sz(L); i++) T[0].push_back({L[i],i});
for(int i=1,l=1; l<sz(L); i++, l=l+l){
T.push_back(vector<pii>());
for(int j=0; j+l+l<sz(L); j++) {
T[i].push_back(min(T[i-1][j], T[i-1][j+l]));
}
}
}
pii getVal(int l, int r) {
if(l>r) return {INF, -1};
int d=1, i=0;
while(d+d<=r-l+1) d=d+d, i++;
return min(T[i][l], T[i][r-d+1]);
}
};
void ConstructEvenSA() {
v.clear();
for(int i=1; i<=elen; i++) v.push_back({{T[i+i], oloc[i+i+1]}, i+i}); //out of index 예외처리
RadixSort(v);
for(int i=1; i<=elen; i++) Ae[i]=v[i-1].second;
for(int i=1; i<elen;i++) {
if(T[Ae[i]]!=T[Ae[i+1]]) Le[i]=0;
else Le[i]=1+OddRMQ(oloc[Ae[i]+1], oloc[Ae[i+1]+1]).first;
}
ERMQ.makeRMQ(Le);
for(int i=1; i<=elen; i++) eloc[Ae[i]]=i; eloc[len+1]=0;
}
3. $SA_{o}$와 $SA_{e}$ 합치기
마지막으로 지금까지 $O(n)+Recursion(n/2)$시간 걸려서 만든 $SA_{o}$와 $SA_{e}$를 $O(n)$시간에 합쳐주기만 하면 됩니다.
가장 단순히 생각할 수 있는 방법은 merge sort에서 사용하는 방법입니다. 양 쪽의 suffix들은 이미 정렬 되어 있으니 앞에서부터 살펴보면서 더 작은 것을 하나하나 채워나가는 것입니다. 이 때 단순히 두 문자열을 계속 비교를 해주면서 넣으면 비교해주는데 최대 $O(n)$시간이 걸리므로 $O(n^{2})$의 시간을 소모하게 됩니다.
이를 줄이기 위해 살짝 다른 방법을 시도 해볼 것입니다. 일단 odd suffix와 even suffix들을 맨 앞자리가 같은 것들끼리 묶어 줄 수 있습니다. 이 후 묶음 통채로 정렬을 해주는 것입니다. 이후 맨 앞자리가 같았던 것들 중에서 두 번째자리를 비교해주고 다시 정렬해주고 할 수 있습니다. 조금 더 효율적인 거 같지만 시간이 $O(n^{2})$걸리는 것은 마찬가지입니다. 어디서 시간이 많이 투자되었을까?
{pppa, pppb, pppc}와 {ppqx, ppqy, ppqz}를 위의 방법으로 합치는 것을 생각해보겠습니다. 첫번째 자리에서 같고, 두번째 자리에서 같고 세번째 자리 까지 가야지 비교가 됩니다. 결국 최악의 경우, 즉 모든 자리가 같은 경우에는 모든 자리를 비교해야되므로 $O(n^{2})$의 시간이 투자되는 것입니다.
여기서 사용하지 않은 정보가 있습니다. Suffix들의 순서, 즉 $A_{o}$와 $A_{e}$는 자연스럽게 사용하고 있었지만 $L_{o}$와 $L_{e}$에 대한 정보는 전혀 이용하지 않고 있었습니다. 앞의 집합의 lcp의 길이가 3인 것과 뒤의 집합의 lcp의 길이가 3인 것을 알고 있다면 ppp와 ppq만 비교해주면 되는 것입니다. 이 모든 과정을 빠르게 하기 위해 몇가지 새로운 개념 하나를 정의하겠습니다.
Equivalence Class
$E_{l}=\{(x,y) | pref_{l}(S_{x}) =pref_{l}(S_{y})\}$
$1 \leq l \leq n$에 대해 정의 되는 equivalence Class입니다. 두 suffix에 대해 길이 $l$의 prefix, 즉 앞의 $l$개의 문자가 같으면 동치라고 정의하는 것입니다.
1. $A_{T}[p, q]$, ($1 \leq p \leq q \leq n$)이 $E_{l}$에서 equivalence class가 될 조건은 다음과 같습니다.
$L_{T}[p-1] < l$, $L_{T}[q] < l$, and $L_T[i]\geq l$ for all $p \leq i < q$
앞의 두 조건이 없으면 더 큰 equivalence class가 더 많은 인덱스를 포함하게 됩니다.
2. $E_{l}$에서 정의 되는 한 equvialence class에 대해, $E_{l+1}$으로의 equivalence class들로 $O(r)$의 시간에 나누어 줄 수 있습니다. 이 때 r은 나누어지는 class의 개수입니다.
$l+1$번째 자리 기준으로 나누어 주는 것입니다. 이 때 일일히 확인해주면 equivalence class의 원소 개수 만큼 시간이 필요하지만 RMQ를 통해 더 빠른 시간내에 해결하는 것입니다. $A_{T}[p,q]$가 $E_{l}$에서 equivalence class라는 것은 $L_{T}[i]$, $p \leq i < q$의 값들이 전부 $l$ 이상이라는 것을 의미합니다. 이 때 $E_{l+1}$으로 나누어 주기 위해서는 값이 정확히 $l$들인 위치들을 찾아야되는데 이는 RMQ를 통해 찾아줄 수 있습니다. 값이 $l$인 위치의 개수, 즉 $r$번 RMQ를 반복해주면 되므로 $O(r)$시간이 걸리게 됩니다.
3. $A_{T}[p,q]$가 $E_{i}$, $a \leq i \leq b$에만 equivalence class일 필요충분조건은 다음과 같습니다.
$\max\{L_{T}[p-1], L_{T}[q]\}=a-1$, $b=|lcp(S_{A_{T}[p]}, S_{A_{T}[q]})|=\min_{p \leq i \leq q-1} L_{T}[i]$
이 때 a와 b를 equivalence class의 start stage와 end stage라 하겠습니다. RMQ를 통해 빠르게 구해줄 수 있습니다.
여기서 end stage가 결정적인 역할을 합니다. end stage까지 prefix가 전부 같으므로 end stage 또는 end stage+1로 넘어가 비교함으로서 같은 부분을 반복하는 과정을 생략해줄 수 있습니다.
이제 equivalence class 개념을 이용하여 merge하는 기본적인 아이디어를 살펴보겠습니다.
$E_{l}$에서 정의 되며 $pref_{l}$이 같은 두 equivalence class $A_{o}[w,x]$와 $A_{e}[y,z]$를 보겠습니다. 두 equivalence class는 합쳐서 $(x-w+1)+(z-y+1)$개의 suffix를 포함하고 있습니다. Merge를 하게 될 경우 이 suffix들은 $A_{T}$에서 정확히 $[w+y-1,x+z]$ 구간을 차지하게 될 것입니다. 따라서 다른 equivalence class와 비교할 필요 없이 두 class들만 잘 비교하여 합쳐주면 됩니다. 이런식으로 prefix가 같은 equivalence class들 끼리 묶어 가며 merge를 해줄텐데 이런 쌍 들을 couple-pair이라 부르겠습니다.
이제 한 coupled pair $C=<A_{o}[w,x], A_{e}[y,z]>$에 대해 $A_{o}[w,x]$의 end stage를 $k_{o}$, $A_{e}[y,z]$의 end stage를 $k_{e}$라 하겠습니다. 그러면 $A_{o}[w,x]$와 $A_{e}[y,z]$를 합치고나면 $A_{T}[w+x-1,y+z]$가 될텐데 해당 equivalence class의 end stage는 최대 $k=\min\{k_{o}, k_{e}\}$가 됩니다. 이 때 $k$를 limit stage라 정의합니다. 두 equivalence class는 각각 $pref_{k}$가 동일한 suffix들만 모여있을 것입니다. 이 $pref_{k}$를 비교하여 만약 다를 경우 한쪽이 나머지 한쪽을 majorize할 것이므로 반반 나누어 넣어주면 됩니다. 만약 $pref_{k}$가 같을 경우 $E_{k+1}$에서의 equivalence class들로 분할 해준뒤 다시 $pref_{k+1}$이 같은 것들끼리 묶어 준 다음에 위의 과정을 반복해주면 됩니다.
이 때 $pref_{k}$를 비교해주는 과정이 살짝 까다롭습니다. 일단 다음과 같이 coupled pair lcp problem을 다음과 같이 정의하고 넘어가겠습니다.
Coupled Pair Lcp Problem
주어진 coupled pair $C=<A_{o}[w,x], A_{e}[y,z]>$에 대해 limit stage가 $k$로 주어져 있다.
이 때, 상수 시간 내에 $C$의 end stage를 계산할 수 있다. 만약 end stage가 $k$보다 작을 경우 누가 더 작은지도 판단해주어야 한다.
이제 위의 과정을 엄밀하게 진행해보겠습니다. 말보다는 코드로 보여드리겠습니다.
void Merge() {
vector<vector<query>> Q; Q.resize(len);
Q[0].push_back({1, olen, 1, elen});
optr[olen]=elen; eptr[elen]=olen;
for(int k=0; k<len; k++) {
for(query &cur :Q[k]) {
int w=cur.w, x=cur.x, y=cur.y, z=cur.z;
int e,d;
if(k==0) e=0;
else if(k==1) e=1;
else {
if(OddRMQ(w,x).first==k) {
int ww=eloc[Ao[w]+1];
int xx=eloc[Ao[x]+1];
int zz=oloc[Ae[z]+1];
if(efin[ww]) e=min(elst[ww], OddRMQ(eptr[ww], zz).first)+1;
else e=OddRMQ(eptr[EvenRMQ(ww,xx).second], zz).first+1;
} else {
int yy=oloc[Ae[y]+1];
int zz=oloc[Ae[z]+1];
int xx=eloc[Ao[x]+1];
if(ofin[yy]) e=min(olst[yy], EvenRMQ(optr[yy], xx).first)+1;
else e=EvenRMQ(optr[OddRMQ(yy,zz).second], xx).first+1;
}
}
e=min(e,k);
if(k==e) {
int o_idx=w;
int e_idx=y;
while(o_idx<=x || e_idx<=z) {
pii to=OddRMQ(o_idx, x);
pii te=EvenRMQ(e_idx, z);
int o_nxt=(to.first==k)? to.second: x;
int e_nxt=(te.first==k)? te.second: z;
int co=(o_idx<=x)? T[Ao[o_idx]+k] : INF;
int ce=(e_idx<=z)? T[Ae[e_idx]+k] : INF;
if(co==ce) {
int M=min(OddRMQ(o_idx,o_nxt).first, EvenRMQ(e_idx,e_nxt).first);
M=min(M,len-1);
Q[M].push_back({o_idx, o_nxt, e_idx, e_nxt});
if(o_nxt+e_nxt<x+z) L[o_nxt+e_nxt]=k;
optr[o_nxt]=e_nxt;
eptr[e_nxt]=o_nxt;
o_idx=o_nxt+1;
e_idx=e_nxt+1;
}
if(co<ce) {
for(int i=o_idx; i<=o_nxt; i++) A[e_idx+i-1]=Ao[i], ofin[i]=1;
for(int i=o_idx; i<o_nxt; i++) L[e_idx+i-1]=Lo[i];
if(e_idx+o_nxt-1<cur.x+cur.z) L[e_idx+o_nxt-1]=k;
for(int i=o_idx; i<=o_nxt; i++) optr[i]=y, olst[i]=e;
o_idx=o_nxt+1;
}
if(co>ce) {
for(int i=e_idx; i<=e_nxt; i++) A[o_idx+i-1]=Ae[i], efin[i]=1;
for(int i=e_idx; i<e_nxt; i++) L[o_idx+i-1]=Le[i];
if(o_idx+e_nxt-1<cur.x+cur.z) L[o_idx+e_nxt-1]=k;
for(int i=e_idx; i<=e_nxt; i++) eptr[i]=w, elst[i]=e;
e_idx=e_nxt+1;
}
}
} else {
if(T[Ao[w]+e]<T[Ae[y]+e]) {
for(int i=w; i<=x; i++) A[y+i-1]=Ao[i], L[y+i-1]=Lo[i], ofin[i]=1;
for(int i=w; i<=x; i++) optr[i]=y, olst[i]=e;
L[x+y-1]=e;
for(int i=y; i<=z; i++) A[x+i]=Ae[i], efin[i]=1;
for(int i=y; i<z; i++) L[x+i]=Le[i];
for(int i=y; i<=z; i++) eptr[i]=w, elst[i]=e;
} else {
for(int i=y; i<=z; i++) A[w+i-1]=Ae[i], L[w+i-1]=Le[i], efin[i]=1;
for(int i=y; i<=z; i++) eptr[i]=w, elst[i]=e;
L[w+z-1]=e;
for(int i=w; i<=x; i++) A[z+i]=Ao[i], ofin[i]=1;
for(int i=w; i<x; i++) L[z+i]=Lo[i];
for(int i=w; i<=x; i++) optr[i]=y, olst[i]=e;
}
}
}
}
}
해당 방법이 맞는 방법인지 살펴보겠습니다. ptr, lst, fin은 나중에 신경쓰도록 하겠습니다. lcp problem을 해결할 때 필요한 배열들입니다.
짝 지어줄 상대가 없는 equivalence class들이 $A_{T}$에 하나씩 채워지게 됩니다. Suffix의 길이가 유한하므로 유한한 스텝을 거치면 common prefix가 없어질 것이고 결국 모든 equivalence class가 자기 자리에 채워질 것입니다. 이 과정에서 정렬이 잘 되는 것은 쉽게 확인할 수 있습니다.
이제 $L_{T}$가 맞게 계산되는지를 확인해야됩니다. $L_{T}$는 다음 두 경우에 그 값을 계산할 것입니다.
첫 번째 경우는 coupled pair의 마지막 원소입니다. 즉, $C=<A_{o}[w,x], A_{e}[y,z]>$에 대해 이 suffix들은 $A_{T}[w+y-1,x+z]$를 차지할텐데 이 때, $L_{T}[x+z]$를 구해줄 것입니다. $C$가 생성된 과정을 다시 짚어보면, end stage가 $k$인 equivalence class를 $k+1$자리를 기준으로 분할하여 합쳐주었었습니다. 따라서 $L_{T}[x+z]$는 자연스럽게 $k$가 되어야 합니다. 이 때 한 가지 처리를 해주어야 하는데 만약 $C$의 equivalence class들이 이전 coupled pair에서 둘 다 마지막을 차지하고 있었다면 k보다 작은 값을 가지게 될 것입니다. 하지만 이 값은 이전 coupled pair에서 연산이 되었을 것이므로 값을 덮어쓰지 않게 주의만 해주면 됩니다.
두 번째는 짝 지어줄 상대가 없는 equivalence class입니다. 짝 지어줄 상대가 없으면 자기 자리에 위치할 것입니다. 이 때 이 equivalence class에 속하는 모든 원소에 대해 $L_{T}$를 계산해줄 것입니다. 한 equivalence class에 속한다는 뜻은 odd 또는 even suffix array에서 연속한 구간을 차지한다는 뜻이고 인접한 suffix들끼리의 lcp는 이미 $L_{o}$ 혹은 $L_{e}$에 계산되어 있을 것입니다. 마지막 suffix의 경우는 coupled pair에서 마지막을 차지하느냐에 따라 위와 마찬가지로 $k$를 저장해주거나 이미 연산된 값을 그대로 두면 됩니다.
최종적으로 모든 equivalence가 제자리를 찾으므로 모든 $L_{T}$가 연산 되며 그 값이 맞다는 것도 확인할 수 있습니다.
이제 시간복잡도를 살펴보겠습니다. 과정이 복잡해보이는데 결국 (coupled pair의 개수)+(분할된 odd equivalence class의 개수)+(분할된 even equivalence class)의 복잡도를 가지게 됩니다. 전부 최대 $O(n)$의 값을 가지므로 합쳐봤자 $O(n)$의 값을 가지게 되는 것입니다.
이제 마지막으로 coupled pair lcp problem을 해결해보겠습니다.
위의 과정에서는merge에서 query들의 벡터 Q를 굳이 limit stage에 따라 여러개로 만들 필요가 없습니다. 어떤 coupled pair끼리 서로 영향을 미치지 않으며 따라서 어떤 것을 먼저 연산하든 상관없기 때문입니다. 저렇게 순서대로 limit stage에 따라 계산을 해준 것은 coupled pair lcp problem을 해결하기 위해서입니다.
Coupled Pair Lcp Problem
주어진 coupled pair $C=<A_{o}[w,x], A_{e}[y,z]>$에 대해 limit stage가 $k$로 주어져 있다.
이 때, 상수 시간 내에 $C$의 end stage를 계산할 수 있다. 만약 end stage가 $k$보다 작을 경우 누가 더 작은지도 판단해주어야 한다.
일단 limit stage $k$가 0인 경우와 1인 경우는 따로 해결하겠습니다.
$k$가 0인 경우는 초기 coupled pair 뿐입니다. end stage는 0일 수 밖에 없습니다. end stage와 limit stage가 같으므로 비교해줄 필요 없습니다.
$k$가 1인 경우는 $A_{o}[w,x]$와 $A_{e}[y,z]$의 end stage가 둘다 1이상, 적어도 하나는 정확히 1인 경우입니다. 즉, 두 equivalence class는 각각 전부 같은 맨 앞자리를 가지고 있다는 것입니다. Coupled pair 생성하는 과정 상 앞자리가 같은 것들끼리 묶어 주었을 것이기 때문에 두 equivalence의 end stage는 1이상일 것입니다. Limit stage보다 클 수 없으므로 end stage가 1이라는 것을 알 수 있습니다. 마찬가지로 비교해줄 필요는 없습니다.
이제 $k$가 2이상인 경우를 생각해보겠습니다. 다음 식이 성립합니다.
$|lcp(P_{A_{o}}(w,x), P_{A_{e}}(y,z))|=|lcp_{k}(P_{A_{o}}(w,x), P_{A_{e}}(y,z))|$
$=|lcp_{k}(S_{A_{o}[w]}, S_{A_{e}[z]})|$
$=|lcp_{k-1}(S_{A_{o}[w]+1}, S_{A_{e}[z]+1})|+1$
$P$는 $A$의 한 구간의 lcp입니다. 첫번째 등호는 limit stage가 $k$이므로 앞의 $k$개만 비교해주면 되기 때문에 성립합니다. 두번 째 등호는 coupled pair와 end stage정의 상 $A_{o}[w,x]$의 앞 k자리가 전부 같기 때문에 성립합니다. 마지막 등호는 두 suffix의 맨 앞자리가 같기 때문에 앞 한자리를 제거한 뒤 새로운 suffix 두개를 비교하는 것과 같다는 것을 의미합니다.
이제 $w^{\prime}=index_{e}(A_{o}[w]+1)$, $x^{\prime}=index_{e}(A_{o}[x]+1)$, $z^{\prime}=index_{o}(A_{e}[z]+1)$이라 정의하고 $S_{A_{e}}[w^{\prime}]$과 $S_{A_{o}}[z^{\prime}]$를 $k-1$자리 까지 비교해주면 됩니다. Odd suffix에서 맨 앞자리를 제거하면 Even suffix가 되는 것을 이용하여 index를 정의 해준 것입니다.
간단히 말해 두 suffix를 비교하기 위해 앞 한자리 씩 빼고 비교해준다는 것입니다. 이 방법이 가능한 이유는 limit stage가 낮은 coupled pair부터 연산했기 때문인데 앞의 한자리를 뺀 suffix는 이미 처리가 되어 있기 때문이라고 추측해 볼 수 있습니다. 이제 RMQ를 이용하여 두 suffix의 lcp를 구해주고 싶은데 아쉽게도 아직 Odd suffix와 Even suffix끼리는 비교를 하지 못합니다.
일반성을 잃지 않고 $A_{o}[w,x]$의 end stage가 k라 하겠습니다. (두 end stage의 최솟값이 k이므로)
이제 새로운 Even Suffix하나를 찾아줄건데 $A_{o}[w^{\prime}]$와 최대한 많이 겹치면 좋겠습니다. 모든 $Ae[1, len/2]$ 중에서 $A_{o}[w^{\prime}]$와 lcp가 가장 긴 Even suffix $A_{e}[r]$을 찾았다 해보겠습니다. 그러면 다음이 성립하는 것을 알 수 있습니다.
$lcp(S_{A_{o}[w^{\prime}]}, S_{A_{e}[x]}) = lcp(lcp(S_{A_{o}[w^{\prime}]}, S_{A_{e}[r]}),lcp(S_{A_{e}[r]}, S_{A_{e}[x]}))$
따라서
|$lcp(S_{A_{o}[w^{\prime}]}, S_{A_{e}[x]})| = \min(|lcp(S_{A_{o}[w^{\prime}]}, S_{A_{e}[r]})|,|lcp(S_{A_{e}[r]}, S_{A_{e}[x]})|)$
두번째 항은 RMQ를 통해 구할 수 있으므로 첫번째 항만 계산해주면 되는데 이는 $r$을 찾으면서 같이 구하게 될 것입니다. 아직 Suffix Array가 완성이 되지 안흔 상황에서 $r$을 찾는 것은 어려울 수도 있습니다. 하지만 여기서 lcp가 최대인 $r$을 찾아줄 필요가 없습니다. $(k-1)$자리까지만 보았을 때 제일 많이 겹치는 것을 찾아주면 되기 때문에 이는 이전 스테이지에서 처리해준 정보들을 통해 구할 수 있습니다. 두가지 케이스로 나누어서 보겠습니다.
1) Odd suffix $A_{o}[w^{\prime}]$가 자기 위치를 찾은 경우입니다. 이 경우 가장 많이 겹치는 Even suffix는 마지막까지 치열하게 겨루던 equivalence class, 즉 마지막 coupled pair에 포함되어있을 것입니다. 자기 위치를 찾는 경우가 두가지 있습니다. 첫 번째는 end stage와 limit stage가 같았지만 $(k+1)$번째 글자를 비교했지만 상대 equivalence class에는 같은 글자가 없는 경우 입니다. 이 경우 상대 equivalence class의 임의의 suffix를 잡아도 lcp가 k가 되고 이 값이 최대인 것을 알 수 있습니다. 두 번째 경우는 end stage가 limit stage보다 작은 경운데 이 경우 또한 마찬가지로 상대 equivalence class의 임의의 suffix를 잡아도 lcp의 길이가 end stage로 동일하며, 이 값이 최대인 것을 확인할 수 있습니다.
2) 자기 자리를 찾지 못한 경우입니다. 이 경우 아직 어떤 equivalence class속에 포함되어 coupled pair을 이루고 있을 것입니다. 이 경우 limit stage가 $k$이상이어야 한다는 것을 알 수 있습니다. ($k$보다 작을 경우 과정 상 이미 처리되었을 것입니다.) 그런데 현재 end stage가 $k$인 equivalence class의 한 suffix의 앞 한자리를 제거한 suffix $S_{A_{o}}[w^{\prime}]=S_{x}$를 보고 있었습니다. 일단 $S_{x}$를 포함하며 end stage가 $(k-1)$인 equivalence class가 존재한다는 것을 보일 수 있습니다.
현재 equivalence class의 end stage가 $k$라는 뜻은 $k+1$번째 자리가 다른 두 suffix가 포함되어 있다는 뜻입니다. 현재 class의 맨 앞자리를 제거한 suffix들도 같은 $(k-1)$ equivalence class의 일부이며 $k$번째 자리가 다른 두 suffix가 포함되게 됩니다. $S_{x}$도 이 새로운 equivalence class에 포함되어야 하므로 위의 claim이 보여집니다.
$(k-1)$ stage에서 처리가 한번 이루어졌어야 하므로 새로운 coupled pair을 이루고 있다는 것은 두 equivalence class가 최소 k자리는 겹친다는 것입니다. 지금 앞 (k-1)자리까지만 보았을 때 제일 많이 겹치는 것을 찾아주면 되기 때문에 상대 equivalence class의 아무 suffix나 잡으면 최소 k자리 겹치고 따라서 충분합니다.
이제 위의 내용을 $O(1)$의 시간에 처리하고 구해주면 됩니다.
1)의 경우는 전부 배열에 저장해줄 것입니다. 모든 suffix는 한번 자기 위치를 찾게 되므로 상대 equivalence class의 임의의 하나와 lcp값을 저장해주면 됩니다. 총 $O(n)$의 시간이 소모될 것입니다. lst에 lcp의 값을, ptr에 상대 suffix하나를 저장해주면 됩니다.
2)의 경우는 처리 안된 suffix가 포함된 coupled pair의 두 equivalence class를 알면 구할 수 있습니다. 하지만 처리가 아직 안됐는데 이 것들을 매번 저장해줄 경우 시간이 과하게 소비 될 것입니다. 따라서 살짝 테크닉을 사용해줄 것입니다. 일단 처리 안된 모든 coupled pair의 두 equivalence class의 마지막 suffix들끼리 포인터로 연결해줄 것입니다. 최대 $n$개의 coulpled pair이 나오므로 이 과정에는 $O(n)$의 시간이 소모 됩니다. 이제 처리 안된 suffix에 대해서 해당 suffix가 포함된 equivalence class만 구해주면 되는데 이는 RMQ를 통해 구해줄 수 있습니다. $S_{A_{o}}[w^{\prime}]$는 $(k-1)$ stage에서 처리가 이루어졌습니다. 따라서 $k$번째 자리가 같은 suffix들끼리 묶였을 것입니다. 이 때 $S_{A_{o}}[x^{\prime}]$를 볼 수 있는데 이 suffix는 위의 suffix와는 $(k-1)$번째 자리까지는 같지만 $k$번째 자리는 다를 것입니다. (end stage가 $(k+1)$이기 때문에) 따라서 $L_{o}[w^{\prime}, x^{\prime}]$의 값은 전부 $(k-1)$이상이며 값이 정확히 $(k-1)$인 위치들을 기준으로 $E_{k}$로 나뉘는 것을 알 수 있습니다. $
S_{A_{o}}[w^{\prime}]$는 이후 처음으로 $L_{o}$의 값이 정확히 $(k-1)$이 되는 위치까지 같은 equivalence class가 될텐데 이 위치는 RMQ를 통해 구할 수 있는 것입니다. 값이 최솟값과 동일하며 첫 번째 인덱스를 구하는 것인데 이 것이 RMQ의 역할과 정확히 일치합니다. 따라서 $L_{o}[w^{\prime}, x^{\prime}]$에서 RMQ를 사용해주고 포인터를 통해 상대 equivalence class의 아무 suffix를 찾아주면 lcp가 k이상인 $r$을 찾을 수 있게 됩니다. 따라서 coupled pair의 마지막끼리 서로 가르키는 포인터를 ptr에 저장해 주면됩니다. 1)의 ptr과 살짝 다른 역할이지만 동일 배열을 사용해도 충돌이 없으므로 동일배열에 저장하게 되었습니다.
말이 너무 복잡해진 것 같습니다. 그래도 위의 내용을 바탕으로 $O(n)$시간에 Suffix Array를 구성하는 방법이 완성되었습니다.
아래는 BOJ 9248번 Suffix Array에서 AC를 받은 전체 코드입니다.
RMQ만 $O(n)$으로 바꿔주면 되지만 가독성이 떨어지므로 그냥 두었습니다.
#include <stdio.h>
#include <iostream>
#include <vector>
#include <algorithm>
#include <queue>
#include <cmath>
#include <string>
#include <cstring>
#include <map>
#include <set>
#include <complex>
#define sz(v) ((int)((v).size()))
#define all(v) (v).begin(), (v).end()
using namespace std;
typedef long long ll;
const int N = 1000001;
const int INF = 0x3c3c3c3c;
const ll LINF = 1ll*INF*INF*2;
typedef pair<int,int> pii;
typedef pair<pii,int> ppi;
void RadixSort(vector<ppi> &v) {
vector<ppi> nv; nv.resize(sz(v));
int lim=max(sz(v)+sz(v),256)+5;
vector<int> cnt=vector<int>(lim, 0);
for(int i=0;i<sz(v);i++) cnt[v[i].first.second]++;
for(int i=1;i<lim;i++) cnt[i]+=cnt[i-1];
for(int i=sz(v)-1; i>=0; i--) nv[--cnt[v[i].first.second]]=v[i];
cnt.clear(); cnt.resize(lim, 0);
for(int i=0;i<sz(v);i++) cnt[nv[i].first.first]++;
for(int i=1;i<lim;i++) cnt[i]+=cnt[i-1];
for(int i=sz(v)-1; i>=0; i--) v[--cnt[nv[i].first.first]]=nv[i];
}
class RMQ {
public:
vector<vector<pii>> T;
void makeRMQ(vector<int> &L) {
T.clear();
T.push_back(vector<pii>());
for(int i=0; i<sz(L); i++) T[0].push_back({L[i],i});
for(int i=1,l=1; l<sz(L); i++, l=l+l){
T.push_back(vector<pii>());
for(int j=0; j+l+l<sz(L); j++) {
T[i].push_back(min(T[i-1][j], T[i-1][j+l]));
}
}
}
pii getVal(int l, int r) {
if(l>r) return {INF, -1};
int d=1, i=0;
while(d+d<=r-l+1) d=d+d, i++;
return min(T[i][l], T[i][r-d+1]);
}
};
class SuffixArray {
public:
int len,olen,elen;
vector<int> T, A, L;
vector<int> Ao,Ae,Lo,Le,oloc,eloc,optr,eptr,ofin,efin,olst,elst;
RMQ ORMQ, ERMQ;
vector<ppi> v;
void init(string &s) {
len=s.size();
T.resize(len+2);
for(int i=1; i<=len; i++) T[i]=s[i-1];
}
void init(int _len, vector<int> &_T) {
len=_len;
T=_T;
}
pii OddRMQ(int l, int r) {
if(l==r) return {len-Ao[l]+1,r};
if(l>r) swap(l,r);r--;
return ORMQ.getVal(l,r);
}
pii EvenRMQ(int l, int r) {
if(l==r) return {len-Ae[l]+1,r};
if(l>r) swap(l,r);r--;
return ERMQ.getVal(l,r);
}
struct query{int w,x,y,z;};
void ConstructOddSA() {
for(int i=1; i<=olen; i++) v.push_back({{T[i+i-1], T[i+i]}, i});
RadixSort(v);
vector<int> nT; nT.resize(olen+2);
int idx=1; nT[v[0].second]=1;
for(int i=1; i<olen; i++) {
if(v[i].first!=v[i-1].first) idx++;
nT[v[i].second]=idx;
}
SuffixArray nSA;
nSA.init(olen,nT);
nSA.SA(Ao,Lo);
for(int i=1; i<=olen; i++) Ao[i]=Ao[i]+Ao[i]-1;
for(int i=1; i<=olen-1; i++) {
Lo[i]=Lo[i]+Lo[i];
if(T[Ao[i]+Lo[i]]==T[Ao[i+1]+Lo[i]]) Lo[i]++;
}
ORMQ.makeRMQ(Lo);
for(int i=1; i<=olen; i++) oloc[Ao[i]]=i; oloc[len+1]=0;
}
void ConstructEvenSA() {
v.clear();
for(int i=1; i<=elen; i++) v.push_back({{T[i+i], oloc[i+i+1]}, i+i});
RadixSort(v);
for(int i=1; i<=elen; i++) Ae[i]=v[i-1].second;
for(int i=1; i<elen;i++) {
if(T[Ae[i]]!=T[Ae[i+1]]) Le[i]=0;
else Le[i]=1+OddRMQ(oloc[Ae[i]+1], oloc[Ae[i+1]+1]).first;
}
ERMQ.makeRMQ(Le);
for(int i=1; i<=elen; i++) eloc[Ae[i]]=i; eloc[len+1]=0;
}
void Merge() {
vector<vector<query>> Q; Q.resize(len);
Q[0].push_back({1, olen, 1, elen});
optr[olen]=elen; eptr[elen]=olen;
for(int k=0; k<len; k++) {
for(query &cur :Q[k]) {
int w=cur.w, x=cur.x, y=cur.y, z=cur.z;
int e,d;
if(k==0) e=0;
else if(k==1) e=1;
else {
if(OddRMQ(w,x).first==k) {
int ww=eloc[Ao[w]+1];
int xx=eloc[Ao[x]+1];
int zz=oloc[Ae[z]+1];
if(efin[ww]) e=min(elst[ww], OddRMQ(eptr[ww], zz).first)+1;
else e=OddRMQ(eptr[EvenRMQ(ww,xx).second], zz).first+1;
} else {
int yy=oloc[Ae[y]+1];
int zz=oloc[Ae[z]+1];
int xx=eloc[Ao[x]+1];
if(ofin[yy]) e=min(olst[yy], EvenRMQ(optr[yy], xx).first)+1;
else e=EvenRMQ(optr[OddRMQ(yy,zz).second], xx).first+1;
}
}
e=min(e,k);
if(k==e) {
int o_idx=w;
int e_idx=y;
while(o_idx<=x || e_idx<=z) {
pii to=OddRMQ(o_idx, x);
pii te=EvenRMQ(e_idx, z);
int o_nxt=(to.first==k)? to.second: x;
int e_nxt=(te.first==k)? te.second: z;
int co=(o_idx<=x)? T[Ao[o_idx]+k] : INF;
int ce=(e_idx<=z)? T[Ae[e_idx]+k] : INF;
if(co==ce) {
int M=min(OddRMQ(o_idx,o_nxt).first, EvenRMQ(e_idx,e_nxt).first);
M=min(M,len-1);
Q[M].push_back({o_idx, o_nxt, e_idx, e_nxt});
if(o_nxt+e_nxt<x+z) L[o_nxt+e_nxt]=k;
optr[o_nxt]=e_nxt;
eptr[e_nxt]=o_nxt;
o_idx=o_nxt+1;
e_idx=e_nxt+1;
}
if(co<ce) {
for(int i=o_idx; i<=o_nxt; i++) A[e_idx+i-1]=Ao[i], ofin[i]=1;
for(int i=o_idx; i<o_nxt; i++) L[e_idx+i-1]=Lo[i];
if(e_idx+o_nxt-1<cur.x+cur.z) L[e_idx+o_nxt-1]=k;
for(int i=o_idx; i<=o_nxt; i++) optr[i]=y, olst[i]=e;
o_idx=o_nxt+1;
}
if(co>ce) {
for(int i=e_idx; i<=e_nxt; i++) A[o_idx+i-1]=Ae[i], efin[i]=1;
for(int i=e_idx; i<e_nxt; i++) L[o_idx+i-1]=Le[i];
if(o_idx+e_nxt-1<cur.x+cur.z) L[o_idx+e_nxt-1]=k;
for(int i=e_idx; i<=e_nxt; i++) eptr[i]=w, elst[i]=e;
e_idx=e_nxt+1;
}
}
} else {
if(T[Ao[w]+e]<T[Ae[y]+e]) {
for(int i=w; i<=x; i++) A[y+i-1]=Ao[i], L[y+i-1]=Lo[i], ofin[i]=1;
for(int i=w; i<=x; i++) optr[i]=y, olst[i]=e;
L[x+y-1]=e;
for(int i=y; i<=z; i++) A[x+i]=Ae[i], efin[i]=1;
for(int i=y; i<z; i++) L[x+i]=Le[i];
for(int i=y; i<=z; i++) eptr[i]=w, elst[i]=e;
} else {
for(int i=y; i<=z; i++) A[w+i-1]=Ae[i], L[w+i-1]=Le[i], efin[i]=1;
for(int i=y; i<=z; i++) eptr[i]=w, elst[i]=e;
L[w+z-1]=e;
for(int i=w; i<=x; i++) A[z+i]=Ao[i], ofin[i]=1;
for(int i=w; i<x; i++) L[z+i]=Lo[i];
for(int i=w; i<=x; i++) optr[i]=y, olst[i]=e;
}
}
}
}
}
void SA(vector<int> &_A, vector<int> &_L) {
if(len==1) {
A.clear(); L.clear();
A.resize(3,0); L.resize(3,0);
A[1]=1;
_A=A;
_L=L;
return;
}
olen=(len+1)/2, elen=len/2;
A.resize(len+1); L.resize(len+1);
Ae.resize(elen+2,0);Le.resize(elen+2,0);
optr.resize(olen+2,0);eptr.resize(elen+2,0);
ofin.resize(olen+2,0);efin.resize(elen+2,0);
olst.resize(olen+2,0);elst.resize(elen+2,0);
oloc.resize(len+2);eloc.resize(len+2);
//Construct odd SA
ConstructOddSA();
ConstructEvenSA();
Merge();
_A=A;
_L=L;
}
};
int main(void){
string s;
cin >> s;
int len=s.size();
vector<int> T, A, L;
SuffixArray sa;
sa.init(s);
sa.SA(A,L);
for(int i=1; i<=len; i++) printf("%d ", A[i]);
printf("\nx ");
for(int i=1; i<len; i++) printf("%d ", L[i]);
// printf("\n");for(int i=1; i<=len; i++) cout<<A[i] <<':'<<s.substr(A[i]-1) <<endl;
}
추가 내용.
RadixSort
void RadixSort(vector<ppi> &v) {
vector<ppi> nv; nv.resize(sz(v));
int lim=max(sz(v),256)+5;
vector<int> cnt=vector<int>(lim, 0);
for(int i=0;i<sz(v);i++) cnt[v[i].first.second]++;
for(int i=1;i<lim;i++) cnt[i]+=cnt[i-1];
for(int i=sz(v)-1; i>=0; i--) nv[--cnt[v[i].first.second]]=v[i];
cnt.clear(); cnt.resize(lim, 0);
for(int i=0;i<sz(v);i++) cnt[nv[i].first.first]++;
for(int i=1;i<lim;i++) cnt[i]+=cnt[i-1];
for(int i=sz(v)-1; i>=0; i--) v[--cnt[nv[i].first.first]]=nv[i];
}배열을 정렬 해줄 때 배열에 포함된 값들의 범위가 상수이거나 $O(n)$ 스케일이면 $O(n)$의 시간에 정렬할 수 있도록 해줍니다. Tuple을 정렬하는 경우 Tuple에 포함된 수의 개수가 $k$개라면 $O(nk)$의 시간이 걸리게 됩니다.
비교 기반 정렬이 아니라 가능한 방법입니다. 각 값들별로 개수를 세서 차례대로 줄 세워주면 됩니다.
출처. Dong Kyue Kim, Jeong Seop Sim, Heejin Park, Kunsoo Park, Constructing suffix arrays in linear time, Journal of Discrete Algorithms, Volume 3, Issues 2–4, 2005, Pages 126-142
'코딩 > 알고리즘&자료구조' 카테고리의 다른 글
| 스도쿠 잡기술 (0) | 2019.12.04 |
|---|---|
| Knuth Optimization (0) | 2019.11.16 |
| 고인물 복잡도 (0) | 2019.11.06 |
| Range Minimum Query (1) | 2019.11.04 |
