Burrows-Wheeler Transform (BWT), Last-to-First (LF) Mapping & FM-Index Substring Search

1. Constructing the BWT from the Suffix Array & Why It Compresses

Given $S[0..N-1]$ with unique smallest sentinel $S[N-1] = \texttt{'\$'}$ and its Suffix Array $\text{SA}[0..N-1]$, we never materialize the $N \times N$ rotation matrix $\mathcal{M}$. Row $r$ of $\mathcal{M}$ starts at $\text{SA}[r]$ and ends at $(\text{SA}[r] - 1 + N) \bmod N$, so the BWT string $L[0..N-1]$ is computed in one linear pass:

Why does $L$ have long runs of identical characters? In English text, suffixes starting with $\texttt{"he "}$ are all grouped adjacently in $\text{SA}$; in almost every one of those suffixes, the character immediately preceding $\texttt{"he "}$ is $\texttt{'t'}$, producing a long contiguous run $\texttt{"tttt...t"}$ inside $L$! Applying Move-to-Front (MTF) coding, Run-Length Encoding (RLE), and Huffman/Arithmetic coding to $L$ is the entire pipeline of bzip2.

2. The Last-to-First (LF) Mapping & $O(N)$ Inversion

Why is the BWT reversible? Consider all rows in $\mathcal{M}$ whose last column $L[r]$ equals a specific character $c$. In those rows, the rest of the row $F[r..] $ spells the suffix immediately following that $c$. Since the rows of $\mathcal{M}$ are sorted by those suffixes $F[r..]$, the relative top-to-bottom order of all $c$'s in $L$ is identical to the relative top-to-bottom order of all rows starting with $c$ in $F$ (where $c$ is followed by the exact same sorted suffixes)!

Let $C[c]$ be the total number of characters in $S$ strictly smaller than $c$ (so the bucket for character $c$ in $F$ starts at row $C[c]$), and let $\text{Occ}(c, r)$ be the number of times character $c$ appears in the prefix $L[0..r-1]$ ($0$-indexed exclusive rank). Then row $r$ in $L$ maps to row $\text{LF}(r)$ in $F$:

Starting from the row $r_0$ where $L[r_0] = \texttt{'\$'}$ (or the row corresponding to $\texttt{'\$'}$ in $F$, which is row $0$), repeatedly stepping $r \leftarrow \text{LF}(r)$ reconstructs $S$ backwards from right to left in $O(N)$ time!

3. FM-Index Backward Search in $O(M)$ Time

Even better, $C[c]$ and $\text{Occ}(c, r)$ allow us to count how many times a pattern $P[0..M-1]$ occurs in $S$ in $O(M)$ steps (faster than binary searching a Suffix Array, which takes $O(M \log N)$)! We process $P$ backwards from $j = M - 1$ down to $0$, maintaining the half-open Suffix Array interval $[sp, ep)$ of suffixes prefixed by $P[j..M-1]$:

After $M$ steps, $ep - sp$ is the exact number of occurrences of $P$ in $S$!

4. Visual Trace: BWT Matrix & LF Mapping on $S = \texttt{"banana\$"}$

5. Production Implementation: BWT Inversion & FM-Index Backward Search

#include <array>
#include <cstdint>
#include <string>
#include <string_view>
#include <utility>
#include <vector>

class FMIndex {
public:
    // Builds BWT and FM-Index from string s (terminated by '' at s.back()) and its SA
    FMIndex(std::string_view s_with_null, const std::vector<int>& sa)
        : n_(static_cast<int>(s_with_null.size())), bwt_(n_, ''), occ_(n_ + 1) {
        c_.fill(0);
        occ_[0].fill(0);
        for (int r = 0; r < n_; ++r) {
            int prev_idx = (sa[r] == 0) ? (n_ - 1) : (sa[r] - 1);
            uint8_t ch = static_cast<uint8_t>(s_with_null[prev_idx]);
            bwt_[r] = static_cast<char>(ch);
            c_[ch]++;
            occ_[r + 1] = occ_[r];
            occ_[r + 1][ch]++;
        }
        // Convert c_ into exclusive prefix sums (bucket start offsets in F)
        int sum = 0;
        for (int ch = 0; ch < 256; ++ch) {
            int cnt = c_[ch];
            c_[ch] = sum;
            sum += cnt;
        }
    }

    // Reconstructs the original string (excluding trailing '') in O(N) via LF-mapping
    std::string invert_bwt() const {
        std::string out(n_ - 1, '');
        int r = 0; // Row 0 in F starts with '', so L[0] is the last char of S
        for (int i = n_ - 2; i >= 0; --i) {
            uint8_t ch = static_cast<uint8_t>(bwt_[r]);
            out[i] = static_cast<char>(ch);
            r = c_[ch] + occ_[r][ch];
        }
        return out;
    }

    // Counts exact occurrences of pat in O(|pat|) time via FM-Index Backward Search
    int count_matches(std::string_view pat) const {
        int sp = 0, ep = n_;
        for (int j = static_cast<int>(pat.size()) - 1; j >= 0 && sp < ep; --j) {
            uint8_t ch = static_cast<uint8_t>(pat[j]);
            sp = c_[ch] + occ_[sp][ch];
            ep = c_[ch] + occ_[ep][ch];
        }
        return std::max(0, ep - sp);
    }

private:
    int n_;
    std::string bwt_;
    std::array<int, 256> c_{};
    std::vector<std::array<int, 256>> occ_;
};