-
Notifications
You must be signed in to change notification settings - Fork 0
/
FM_index.cpp
54 lines (52 loc) · 1.36 KB
/
FM_index.cpp
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
#pragma once
#include "../structure/WaveletMatrix.cpp"
#include "../template/template.cpp"
#include "BWT.cpp"
template <class T, class C>
class FMIndex {
ll N, base;
T bwt;
vector<ll> c;
WaveletMatrix<T, C> WM;
SuffixArray<T> SA;
public:
T ST;
PL occ(T &S) {
for (auto i : S)
if ((ll)i < base || (ll)i - base >= len(c)) return PL(0, 0);
ll sp = 0, ep = N;
rev(i, len(S)) {
sp = c[(ll)S[i] - base] + WM.rank(S[i], sp);
ep = c[(ll)S[i] - base] + WM.rank(S[i], ep);
if (sp >= ep) return PL(0, 0);
}
return PL(sp, ep);
}
vector<ll> locate(T &S) {
vector<bool> v(len(ST) + 1);
PL range = occ(S);
for (ll i = range.first; i < range.second; i++) v[SA[i]] = true;
vector<ll> res;
rep(i, len(ST) + 1) if (v[i]) res.emplace_back(i);
return res;
}
FMIndex(T S) : N(len(S) + 1), ST(S + '$'), WM("", 0), SA(S) {
bwt = BWT(S, SA);
WM = WaveletMatrix<T, C>(bwt, 8);
ll mn = inf, mx = -inf;
for (C i : ST) {
chmin(mn, (ll)i);
chmax(mx, (ll)i);
}
c.resize(mx - mn + 2);
for (C i : ST) {
c[(ll)i - mn + 1]++;
}
rep(i, mx - mn + 1) c[i + 1] += c[i];
base = mn;
}
};
/*
@brief FM Index
@docs docs/FM_index.md
*/