Myers O(ND) Shortest Edit Script & Hirschberg's Linear-Space Alignment

1. Myers Edit Graph: Diagonals $k = x - y$ & Greedy Snakes

Given sequences $A[0..N-1]$ and $B[0..M-1]$, consider the grid of vertices $(x, y) \in [0, N] \times [0, M]$. Define diagonal $k$ as the set of grid points satisfying $k = x - y$:

• Deletion from $A$ ($(x, y) \to (x + 1, y)$): Increases $x$ by $1$, moving from diagonal $k - 1$ to diagonal $k$.
• Insertion from $B$ ($(x, y) \to (x, y + 1)$): Increases $y$ by $1$, moving from diagonal $k + 1$ to diagonal $k$.
• Match / Snake ($(x, y) \to (x + 1, y + 1)$ when $A[x] == B[y]$): Increases both $x$ and $y$ by $1$, sliding down the same diagonal $k = x - y$ at $0$ edit cost!

Let $V_d[k]$ be the maximum $x$-coordinate reachable on diagonal $k$ using at most $d$ edits. Because $d$ edits starting from $(0, 0)$ (diagonal $0$) change the diagonal by $\pm 1$ on each edit, after $d$ edits only diagonals $k \in [-d, +d]$ with $k \equiv d \pmod 2$ are reachable! To reach diagonal $k$ at step $d$, we must have taken either a vertical insertion from diagonal $k + 1$ (preserving $x = V_{d-1}[k + 1]$) or a horizontal deletion from diagonal $k - 1$ (advancing $x = V_{d-1}[k - 1] + 1$), followed by a greedy diagonal slide (snake) while $A[x] == B[y]$:

2. Hirschberg's Linear-Space Divide-and-Conquer ($O(N + M)$ Memory)

Whether using classical LCS DP or Myers diff, reconstructing the actual edit script naively requires storing the history of DP rows ($O(N \cdot M)$ or $O((N + M)D)$ space). Dan Hirschberg (1975) showed how to reconstruct the exact optimal alignment in $O(N + M)$ space while only doubling the runtime constant:

1. Split $A$ at its midpoint $x_{\text{mid}} = \lfloor N / 2 \rfloor$.
2. Run 2-row space-optimized LCS/edit-distance forward on $A[0..x_{\text{mid}}-1]$ against $B[0..M-1]$ to obtain row $L_1[0..M]$.
3. Run 2-row space-optimized LCS/edit-distance backwards on the reversed suffix $A[x_{\text{mid}}..N-1]$ against $B[0..M-1]$ to obtain row $L_2[0..M]$.
4. Find $y^* = \arg\max_{0 \le y \le M} (L_1[y] + L_2[y])$ (or $\arg\min$ for edit distance). Then $(x_{\text{mid}}, y^*)$ is guaranteed to lie on an optimal alignment path! Recurse on $(A[0..x_{\text{mid}}-1],\; B[0..y^*-1])$ and $(A[x_{\text{mid}}..N-1],\; B[y^*..M-1])$.

Because the subproblem areas sum to $\frac{N}{2} y^* + \frac{N}{2}(M - y^*) = \frac{1}{2} NM$, the geometric series of work across all recursion levels is $NM(1 + \frac{1}{2} + \frac{1}{4} + \dots) = 2NM = O(NM)$ time in $O(N + M)$ space!

3. Visual Trace: Myers Edit Graph & Hirschberg Midpoint Split

4. Production Implementation: Myers $O(ND)$ Shortest Edit Distance & Hirschberg LCS

#include <algorithm>
#include <string>
#include <string_view>
#include <vector>

// 1. Myers O((N + M) * D) Shortest Edit Script Distance (Ins + Del, where Match = 0)
int myers_edit_distance(std::string_view a, std::string_view b) {
    const int n = static_cast<int>(a.size());
    const int m = static_cast<int>(b.size());
    const int max_d = n + m;
    if (max_d == 0) return 0;

    std::vector<int> v(2 * max_d + 3, 0);
    const int offset = max_d + 1;

    for (int d = 0; d <= max_d; ++d) {
        for (int k = -d; k <= d; k += 2) {
            int x;
            if (k == -d || (k != d && v[offset + k - 1] < v[offset + k + 1])) {
                x = v[offset + k + 1];     // Vertical step (Insert from B)
            } else {
                x = v[offset + k - 1] + 1; // Horizontal step (Delete from A)
            }
            int y = x - k;
            while (x < n && y < m && a[x] == b[y]) {
                ++x;
                ++y;
            }
            v[offset + k] = x;
            if (x >= n && y >= m) return d;
        }
    }
    return max_d;
}

// 2. Hirschberg's O(N * M) Time, O(N + M) Space LCS String Reconstruction
std::vector<int> lcs_row(std::string_view a, std::string_view b, bool rev) {
    const int n = static_cast<int>(a.size());
    const int m = static_cast<int>(b.size());
    std::vector<int> dp(m + 1, 0), nxt(m + 1, 0);
    for (int i = 0; i < n; ++i) {
        char ca = rev ? a[n - 1 - i] : a[i];
        for (int j = 0; j < m; ++j) {
            char cb = rev ? b[m - 1 - j] : b[j];
            nxt[j + 1] = (ca == cb) ? (dp[j] + 1) : std::max(dp[j + 1], nxt[j]);
        }
        dp.swap(nxt);
    }
    return dp;
}

std::string hirschberg_lcs(std::string_view a, std::string_view b) {
    if (a.empty() || b.empty()) return "";
    if (a.size() == 1) {
        return (b.find(a[0]) != std::string_view::npos) ? std::string(a) : "";
    }
    size_t mid = a.size() / 2;
    auto l1 = lcs_row(a.substr(0, mid), b, false);
    auto l2 = lcs_row(a.substr(mid), b, true);

    size_t best_k = 0;
    int best_val = -1;
    for (size_t k = 0; k <= b.size(); ++k) {
        int val = l1[k] + l2[b.size() - k];
        if (val > best_val) {
            best_val = val;
            best_k = k;
        }
    }
    return hirschberg_lcs(a.substr(0, mid), b.substr(0, best_k)) +
           hirschberg_lcs(a.substr(mid), b.substr(best_k));
}