graph_algorithms.cpp

// graph_algorithms.cpp
// Compile: g++ -std=c++17 -O2 -Wall -Wextra -o graph_algorithms graph_algorithms.cpp
// Run:     ./graph_algorithms

#include <iostream>
#include <vector>
#include <queue>
#include <stack>
#include <deque>
#include <algorithm>
#include <numeric>
#include <climits>
#include <functional>
#include <iomanip>
#include <utility>

using namespace std;

// ============================================================
// 1. Graph class — adjacency list with weighted edges
// ============================================================
class Graph {
public:
    int n;
    vector<vector<pair<int, int>>> adj;  // adj[u] = {(neighbor, weight), ...}
    bool directed;

    Graph(int n_, bool directed_ = false) : n(n_), adj(n_), directed(directed_) {}

    void addEdge(int u, int v, int w = 1) {
        adj[u].push_back({v, w});
        if (!directed) adj[v].push_back({u, w});
    }

    int size() const { return n; }
};

// ============================================================
// 2. Breadth-First Search — single-source distance in unweighted graph
// ============================================================
vector<int> bfs(const Graph& g, int s) {
    vector<int> dist(g.n, -1);
    queue<int> q;
    q.push(s);
    dist[s] = 0;
    while (!q.empty()) {
        int u = q.front(); q.pop();
        for (const auto& e : g.adj[u]) {
            int v = e.first;
            if (dist[v] == -1) {
                dist[v] = dist[u] + 1;
                q.push(v);
            }
        }
    }
    return dist;
}

// ============================================================
// 3. Depth-First Search — iterative + recursive variants
// ============================================================
void dfsRec(const Graph& g, int u, vector<bool>& visited, vector<int>& order) {
    visited[u] = true;
    order.push_back(u);
    for (const auto& e : g.adj[u])
        if (!visited[e.first]) dfsRec(g, e.first, visited, order);
}

vector<int> dfs(const Graph& g, int s) {
    vector<bool> visited(g.n, false);
    vector<int> order;
    dfsRec(g, s, visited, order);
    return order;
}

// ============================================================
// 4. Topological Sort — Kahn's algorithm (BFS-based)
// ============================================================
vector<int> topoSort(const Graph& g) {
    vector<int> inDeg(g.n, 0);
    for (int u = 0; u < g.n; ++u)
        for (const auto& e : g.adj[u]) ++inDeg[e.first];
    queue<int> q;
    for (int i = 0; i < g.n; ++i) if (inDeg[i] == 0) q.push(i);
    vector<int> order;
    while (!q.empty()) {
        int u = q.front(); q.pop();
        order.push_back(u);
        for (const auto& e : g.adj[u])
            if (--inDeg[e.first] == 0) q.push(e.first);
    }
    return (int)order.size() == g.n ? order : vector<int>();  // empty -> cycle
}

// ============================================================
// 5. Cycle Detection (directed graph) — DFS three colors
//    color[v]: 0 = unvisited, 1 = on stack, 2 = finished
// ============================================================
bool dfsCycle(const Graph& g, int u, vector<int>& color) {
    color[u] = 1;
    for (const auto& e : g.adj[u]) {
        int v = e.first;
        if (color[v] == 1) return true;
        if (color[v] == 0 && dfsCycle(g, v, color)) return true;
    }
    color[u] = 2;
    return false;
}

bool hasCycleDirected(const Graph& g) {
    vector<int> color(g.n, 0);
    for (int i = 0; i < g.n; ++i)
        if (color[i] == 0 && dfsCycle(g, i, color)) return true;
    return false;
}

// ============================================================
// 6. Connected Components (undirected) — DFS labelling
// ============================================================
vector<int> connectedComponents(const Graph& g) {
    vector<int> comp(g.n, -1);
    int cid = 0;
    for (int s = 0; s < g.n; ++s) {
        if (comp[s] != -1) continue;
        stack<int> st;
        st.push(s);
        comp[s] = cid;
        while (!st.empty()) {
            int u = st.top(); st.pop();
            for (const auto& e : g.adj[u])
                if (comp[e.first] == -1) {
                    comp[e.first] = cid;
                    st.push(e.first);
                }
        }
        ++cid;
    }
    return comp;
}

// ============================================================
// 7. Strongly Connected Components — Tarjan's algorithm
// ============================================================
class TarjanSCC {
    const Graph& g;
    int timer = 0, sccCount = 0;
    vector<int> disc, low, sccId;
    vector<bool> onStack;
    stack<int> stk;

    void dfs(int u) {
        disc[u] = low[u] = ++timer;
        stk.push(u);
        onStack[u] = true;
        for (const auto& e : g.adj[u]) {
            int v = e.first;
            if (!disc[v]) { dfs(v); low[u] = min(low[u], low[v]); }
            else if (onStack[v]) low[u] = min(low[u], disc[v]);
        }
        if (disc[u] == low[u]) {
            while (true) {
                int v = stk.top(); stk.pop();
                onStack[v] = false;
                sccId[v] = sccCount;
                if (v == u) break;
            }
            ++sccCount;
        }
    }
public:
    TarjanSCC(const Graph& g_)
        : g(g_), disc(g_.n, 0), low(g_.n, 0), sccId(g_.n, -1),
          onStack(g_.n, false) {
        for (int i = 0; i < g.n; ++i) if (!disc[i]) dfs(i);
    }
    int count() const { return sccCount; }
    const vector<int>& component() const { return sccId; }
};

// ============================================================
// 8. Articulation Points & Bridges — Tarjan
// ============================================================
class ArticulationBridge {
    const Graph& g;
    int timer = 0;
    vector<int> disc, low;
public:
    vector<bool> isArticulation;
    vector<pair<int, int>> bridges;

    ArticulationBridge(const Graph& g_)
        : g(g_), disc(g_.n, 0), low(g_.n, 0), isArticulation(g_.n, false) {
        for (int i = 0; i < g.n; ++i) if (!disc[i]) dfs(i, -1);
    }
private:
    void dfs(int u, int parent) {
        disc[u] = low[u] = ++timer;
        int children = 0;
        for (const auto& e : g.adj[u]) {
            int v = e.first;
            if (!disc[v]) {
                ++children;
                dfs(v, u);
                low[u] = min(low[u], low[v]);
                if (low[v] > disc[u]) bridges.push_back({u, v});
                if (parent != -1 && low[v] >= disc[u]) isArticulation[u] = true;
            } else if (v != parent) {
                low[u] = min(low[u], disc[v]);
            }
        }
        if (parent == -1 && children > 1) isArticulation[u] = true;
    }
};

// ============================================================
// 9. Disjoint Set Union (Union-Find) with path compression + rank
// ============================================================
class DSU {
    vector<int> parent, rnk;
public:
    DSU(int n) : parent(n), rnk(n, 0) {
        iota(parent.begin(), parent.end(), 0);
    }
    int find(int x) {
        while (parent[x] != x) {
            parent[x] = parent[parent[x]];
            x = parent[x];
        }
        return x;
    }
    bool unite(int x, int y) {
        int rx = find(x), ry = find(y);
        if (rx == ry) return false;
        if (rnk[rx] < rnk[ry]) swap(rx, ry);
        parent[ry] = rx;
        if (rnk[rx] == rnk[ry]) ++rnk[rx];
        return true;
    }
    bool same(int x, int y) { return find(x) == find(y); }
};

// ============================================================
// 10. Dijkstra — Single-source shortest path (non-negative weights)
// ============================================================
vector<long long> dijkstra(const Graph& g, int s) {
    const long long INF = LLONG_MAX / 4;
    vector<long long> dist(g.n, INF);
    using P = pair<long long, int>;
    priority_queue<P, vector<P>, greater<P>> pq;
    dist[s] = 0;
    pq.push({0, s});
    while (!pq.empty()) {
        auto [d, u] = pq.top(); pq.pop();
        if (d > dist[u]) continue;
        for (const auto& e : g.adj[u]) {
            long long nd = d + e.second;
            if (nd < dist[e.first]) {
                dist[e.first] = nd;
                pq.push({nd, e.first});
            }
        }
    }
    return dist;
}

// ============================================================
// 11. 0-1 BFS — shortest path when weights are only 0 or 1
// ============================================================
vector<int> zeroOneBFS(const Graph& g, int s) {
    vector<int> dist(g.n, INT_MAX);
    deque<int> dq;
    dist[s] = 0;
    dq.push_front(s);
    while (!dq.empty()) {
        int u = dq.front(); dq.pop_front();
        for (const auto& e : g.adj[u]) {
            int nd = dist[u] + e.second;
            if (nd < dist[e.first]) {
                dist[e.first] = nd;
                if (e.second == 0) dq.push_front(e.first);
                else dq.push_back(e.first);
            }
        }
    }
    return dist;
}

// ============================================================
// 12. Bellman-Ford — handles negative weights, detects negative cycle
// ============================================================
struct BellmanFordResult {
    vector<long long> dist;
    bool hasNegativeCycle;
};

BellmanFordResult bellmanFord(const Graph& g, int s) {
    const long long INF = LLONG_MAX / 4;
    vector<long long> dist(g.n, INF);
    dist[s] = 0;
    for (int i = 0; i < g.n - 1; ++i) {
        bool relaxed = false;
        for (int u = 0; u < g.n; ++u) {
            if (dist[u] == INF) continue;
            for (const auto& e : g.adj[u]) {
                if (dist[u] + e.second < dist[e.first]) {
                    dist[e.first] = dist[u] + e.second;
                    relaxed = true;
                }
            }
        }
        if (!relaxed) break;
    }
    bool neg = false;
    for (int u = 0; u < g.n && !neg; ++u) {
        if (dist[u] == INF) continue;
        for (const auto& e : g.adj[u])
            if (dist[u] + e.second < dist[e.first]) { neg = true; break; }
    }
    return {dist, neg};
}

// ============================================================
// 13. Floyd-Warshall — all-pairs shortest paths O(V^3)
// ============================================================
vector<vector<long long>> floydWarshall(const Graph& g) {
    const long long INF = LLONG_MAX / 4;
    vector<vector<long long>> d(g.n, vector<long long>(g.n, INF));
    for (int i = 0; i < g.n; ++i) d[i][i] = 0;
    for (int u = 0; u < g.n; ++u)
        for (const auto& e : g.adj[u])
            d[u][e.first] = min(d[u][e.first], (long long)e.second);
    for (int k = 0; k < g.n; ++k)
        for (int i = 0; i < g.n; ++i)
            for (int j = 0; j < g.n; ++j)
                if (d[i][k] < INF && d[k][j] < INF
                    && d[i][k] + d[k][j] < d[i][j])
                    d[i][j] = d[i][k] + d[k][j];
    return d;
}

// ============================================================
// 14. Kruskal — Minimum spanning tree via DSU
// ============================================================
struct EdgeWUV { int u, v, w; };

pair<long long, vector<EdgeWUV>> kruskal(int n, vector<EdgeWUV> edges) {
    sort(edges.begin(), edges.end(),
         [](const EdgeWUV& a, const EdgeWUV& b) { return a.w < b.w; });
    DSU dsu(n);
    long long total = 0;
    vector<EdgeWUV> mst;
    for (const auto& e : edges) {
        if (dsu.unite(e.u, e.v)) {
            total += e.w;
            mst.push_back(e);
        }
    }
    return {total, mst};
}

// ============================================================
// 15. Prim — Minimum spanning tree via priority queue
// ============================================================
long long primMST(const Graph& g) {
    long long total = 0;
    vector<bool> inTree(g.n, false);
    using P = pair<int, int>;  // (weight, node)
    priority_queue<P, vector<P>, greater<P>> pq;
    pq.push({0, 0});
    int taken = 0;
    while (!pq.empty() && taken < g.n) {
        auto [w, u] = pq.top(); pq.pop();
        if (inTree[u]) continue;
        inTree[u] = true;
        total += w;
        ++taken;
        for (const auto& e : g.adj[u])
            if (!inTree[e.first]) pq.push({e.second, e.first});
    }
    return total;
}

// ============================================================
// 16. Edmonds-Karp — max flow via BFS-augmenting paths
// ============================================================
class MaxFlow {
public:
    struct Edge { int to, rev, cap; };
    int n;
    vector<vector<Edge>> g;

    MaxFlow(int n_) : n(n_), g(n_) {}

    void addEdge(int u, int v, int cap) {
        g[u].push_back({v, (int)g[v].size(), cap});
        g[v].push_back({u, (int)g[u].size() - 1, 0});
    }

    int bfsAug(int s, int t, vector<int>& parent, vector<int>& parentEdge) {
        parent.assign(n, -1);
        parentEdge.assign(n, -1);
        parent[s] = s;
        queue<pair<int, int>> q;
        q.push({s, INT_MAX});
        while (!q.empty()) {
            auto [u, flow] = q.front(); q.pop();
            for (int i = 0; i < (int)g[u].size(); ++i) {
                const auto& e = g[u][i];
                if (parent[e.to] == -1 && e.cap > 0) {
                    parent[e.to]     = u;
                    parentEdge[e.to] = i;
                    int newFlow = min(flow, e.cap);
                    if (e.to == t) return newFlow;
                    q.push({e.to, newFlow});
                }
            }
        }
        return 0;
    }

    int run(int s, int t) {
        int totalFlow = 0;
        vector<int> parent, parentEdge;
        int pushed;
        while ((pushed = bfsAug(s, t, parent, parentEdge)) > 0) {
            totalFlow += pushed;
            int cur = t;
            while (cur != s) {
                int prev = parent[cur];
                int ei = parentEdge[cur];
                g[prev][ei].cap -= pushed;
                g[cur][g[prev][ei].rev].cap += pushed;
                cur = prev;
            }
        }
        return totalFlow;
    }
};

// ============================================================
// 17. Bipartite Matching — Hungarian-style augmenting paths
// ============================================================
class BipartiteMatching {
    int nL, nR;
    vector<vector<int>> adj;  // adj[L-side u] = list of R-side neighbors
    vector<int> matchR;       // matchR[R] = L matched to R
    vector<bool> usedL;
    bool tryMatch(int u) {
        for (int v : adj[u]) {
            if (usedL[v]) continue;
            usedL[v] = true;
            if (matchR[v] == -1 || tryMatch(matchR[v])) {
                matchR[v] = u;
                return true;
            }
        }
        return false;
    }
public:
    BipartiteMatching(int nL_, int nR_) : nL(nL_), nR(nR_), adj(nL_), matchR(nR_, -1) {}
    void addEdge(int u, int v) { adj[u].push_back(v); }
    int run() {
        int total = 0;
        for (int u = 0; u < nL; ++u) {
            usedL.assign(nR, false);
            if (tryMatch(u)) ++total;
        }
        return total;
    }
    const vector<int>& matchingR() const { return matchR; }
};

// ============================================================
// 18. LCA — Lowest Common Ancestor via binary lifting
// ============================================================
class LCA {
    int n, LOG;
    vector<vector<int>> up;   // up[k][v] = 2^k-th ancestor of v
    vector<int> depth;
public:
    LCA(int n_, const vector<vector<int>>& tree, int root)
        : n(n_), LOG(1), depth(n_, 0) {
        while ((1 << LOG) < n) ++LOG;
        up.assign(LOG + 1, vector<int>(n, root));
        function<void(int, int)> dfs = [&](int u, int par) {
            up[0][u] = par;
            for (int v : tree[u]) {
                if (v == par) continue;
                depth[v] = depth[u] + 1;
                dfs(v, u);
            }
        };
        dfs(root, root);
        for (int k = 1; k <= LOG; ++k)
            for (int v = 0; v < n; ++v)
                up[k][v] = up[k - 1][up[k - 1][v]];
    }
    int query(int u, int v) {
        if (depth[u] < depth[v]) swap(u, v);
        int diff = depth[u] - depth[v];
        for (int k = 0; k <= LOG; ++k)
            if ((diff >> k) & 1) u = up[k][u];
        if (u == v) return u;
        for (int k = LOG; k >= 0; --k)
            if (up[k][u] != up[k][v]) {
                u = up[k][u];
                v = up[k][v];
            }
        return up[0][u];
    }
    int dist(int u, int v) { return depth[u] + depth[v] - 2 * depth[query(u, v)]; }
};

// ============================================================
// Demonstrations
// ============================================================
static void hr() { cout << string(60, '-') << "\n"; }

int main() {
    cout << fixed << setprecision(2);

    // ---- BFS demo ----
    hr(); cout << "1. BFS shortest distance (unweighted)\n";
    {
        Graph g(6);
        g.addEdge(0, 1); g.addEdge(0, 2);
        g.addEdge(1, 3); g.addEdge(2, 3);
        g.addEdge(3, 4); g.addEdge(4, 5);
        auto d = bfs(g, 0);
        cout << "  dist from 0 = ";
        for (int x : d) cout << x << " ";
        cout << "\n";
    }

    // ---- DFS demo ----
    hr(); cout << "2. DFS pre-order from 0\n";
    {
        Graph g(6);
        g.addEdge(0, 1); g.addEdge(0, 2);
        g.addEdge(1, 3); g.addEdge(2, 4); g.addEdge(4, 5);
        auto o = dfs(g, 0);
        cout << "  order: ";
        for (int x : o) cout << x << " ";
        cout << "\n";
    }

    // ---- Topological sort ----
    hr(); cout << "3. Topological sort (Kahn): 0->1, 0->2, 1->3, 2->3, 3->4, 2->5\n";
    {
        Graph g(6, true);
        g.addEdge(0, 1); g.addEdge(0, 2);
        g.addEdge(1, 3); g.addEdge(2, 3);
        g.addEdge(3, 4); g.addEdge(2, 5);
        auto t = topoSort(g);
        cout << "  order: ";
        for (int x : t) cout << x << " ";
        cout << "\n";
    }

    // ---- Cycle detection ----
    hr(); cout << "4. Cycle detection (directed)\n";
    {
        Graph g(4, true);
        g.addEdge(0, 1); g.addEdge(1, 2); g.addEdge(2, 0); g.addEdge(2, 3);
        cout << "  has cycle = " << boolalpha << hasCycleDirected(g) << "\n";
    }
    {
        Graph g(4, true);
        g.addEdge(0, 1); g.addEdge(1, 2); g.addEdge(2, 3);
        cout << "  has cycle (DAG) = " << boolalpha << hasCycleDirected(g) << "\n";
    }

    // ---- SCC (Tarjan) ----
    hr(); cout << "5. SCC (Tarjan): 0->1->2->0, 1->3, 3->4, 4->5, 5->3\n";
    {
        Graph g(6, true);
        g.addEdge(0, 1); g.addEdge(1, 2); g.addEdge(2, 0);
        g.addEdge(1, 3);
        g.addEdge(3, 4); g.addEdge(4, 5); g.addEdge(5, 3);
        TarjanSCC scc(g);
        cout << "  components count = " << scc.count() << "\n";
        cout << "  comp[v] = ";
        for (int x : scc.component()) cout << x << " ";
        cout << "\n";
    }

    // ---- Articulation & Bridges ----
    hr(); cout << "6. Articulation points & bridges\n";
    {
        Graph g(7);
        g.addEdge(0, 1); g.addEdge(1, 2); g.addEdge(2, 0);  // triangle
        g.addEdge(2, 3); g.addEdge(3, 4); g.addEdge(4, 5);  // chain
        g.addEdge(5, 6); g.addEdge(6, 4);                   // small cycle
        ArticulationBridge ab(g);
        cout << "  articulation = ";
        for (int i = 0; i < g.n; ++i)
            if (ab.isArticulation[i]) cout << i << " ";
        cout << "\n  bridges = ";
        for (auto& p : ab.bridges) cout << "(" << p.first << "," << p.second << ") ";
        cout << "\n";
    }

    // ---- Connected components ----
    hr(); cout << "7. Connected components\n";
    {
        Graph g(6);
        g.addEdge(0, 1); g.addEdge(2, 3); g.addEdge(3, 4);
        auto cc = connectedComponents(g);
        cout << "  comp = ";
        for (int x : cc) cout << x << " ";
        cout << "\n";
    }

    // ---- Dijkstra ----
    hr(); cout << "8. Dijkstra (single-source shortest path)\n";
    {
        Graph g(5, true);
        g.addEdge(0, 1, 4); g.addEdge(0, 2, 1);
        g.addEdge(2, 1, 2); g.addEdge(1, 3, 1);
        g.addEdge(2, 3, 5); g.addEdge(3, 4, 3);
        auto d = dijkstra(g, 0);
        cout << "  dist from 0 = ";
        for (auto x : d) cout << x << " ";
        cout << "\n";
    }

    // ---- 0-1 BFS ----
    hr(); cout << "9. 0-1 BFS\n";
    {
        Graph g(5, true);
        g.addEdge(0, 1, 0); g.addEdge(0, 2, 1);
        g.addEdge(1, 2, 0); g.addEdge(2, 3, 1); g.addEdge(3, 4, 0);
        auto d = zeroOneBFS(g, 0);
        cout << "  dist from 0 = ";
        for (auto x : d) cout << x << " ";
        cout << "\n";
    }

    // ---- Bellman-Ford (negative edge demo) ----
    hr(); cout << "10. Bellman-Ford (with negative edge)\n";
    {
        Graph g(5, true);
        g.addEdge(0, 1, 6);  g.addEdge(0, 2, 7);
        g.addEdge(1, 2, 8);  g.addEdge(1, 3, 5);  g.addEdge(1, 4, -4);
        g.addEdge(2, 3, -3); g.addEdge(2, 4, 9);
        g.addEdge(3, 1, -2);
        g.addEdge(4, 0, 2);  g.addEdge(4, 3, 7);
        auto r = bellmanFord(g, 0);
        cout << "  dist from 0 = ";
        for (auto x : r.dist) cout << x << " ";
        cout << "\n  has negative cycle = " << boolalpha << r.hasNegativeCycle << "\n";
    }

    // ---- Floyd-Warshall ----
    hr(); cout << "11. Floyd-Warshall (all pairs)\n";
    {
        Graph g(4, true);
        g.addEdge(0, 1, 5);  g.addEdge(0, 3, 10);
        g.addEdge(1, 2, 3);  g.addEdge(2, 3, 1);
        const long long INF = LLONG_MAX / 4;
        auto d = floydWarshall(g);
        for (int i = 0; i < g.n; ++i) {
            cout << "  ";
            for (int j = 0; j < g.n; ++j) {
                if (d[i][j] >= INF) cout << setw(5) << "INF";
                else                cout << setw(5) << d[i][j];
            }
            cout << "\n";
        }
    }

    // ---- Kruskal MST ----
    hr(); cout << "12. Kruskal MST\n";
    {
        vector<EdgeWUV> edges = {
            {0, 1, 4}, {0, 2, 3}, {1, 2, 1}, {1, 3, 2},
            {2, 3, 4}, {3, 4, 2}, {4, 5, 6}
        };
        auto [w, mst] = kruskal(6, edges);
        cout << "  total weight = " << w << "\n  MST edges:";
        for (auto& e : mst) cout << " (" << e.u << "," << e.v << "," << e.w << ")";
        cout << "\n";
    }

    // ---- Prim MST ----
    hr(); cout << "13. Prim MST\n";
    {
        Graph g(6);
        g.addEdge(0, 1, 4); g.addEdge(0, 2, 3);
        g.addEdge(1, 2, 1); g.addEdge(1, 3, 2);
        g.addEdge(2, 3, 4); g.addEdge(3, 4, 2);
        g.addEdge(4, 5, 6);
        cout << "  total weight = " << primMST(g) << "\n";
    }

    // ---- DSU ----
    hr(); cout << "14. Disjoint Set Union\n";
    {
        DSU d(6);
        d.unite(0, 1); d.unite(2, 3); d.unite(1, 2);
        cout << "  same(0,3)? " << boolalpha << d.same(0, 3) << "\n";
        cout << "  same(0,4)? " << boolalpha << d.same(0, 4) << "\n";
    }

    // ---- Edmonds-Karp Max Flow ----
    hr(); cout << "15. Edmonds-Karp max flow (CLRS classic)\n";
    {
        MaxFlow mf(6);
        // 0=s, 5=t
        mf.addEdge(0, 1, 16); mf.addEdge(0, 2, 13);
        mf.addEdge(1, 2, 10); mf.addEdge(2, 1, 4);
        mf.addEdge(1, 3, 12); mf.addEdge(3, 2, 9);
        mf.addEdge(2, 4, 14); mf.addEdge(4, 3, 7);
        mf.addEdge(3, 5, 20); mf.addEdge(4, 5, 4);
        cout << "  max flow s->t = " << mf.run(0, 5) << "\n";
    }

    // ---- Bipartite Matching ----
    hr(); cout << "16. Bipartite matching (4 jobs, 4 workers)\n";
    {
        BipartiteMatching m(4, 4);
        m.addEdge(0, 0); m.addEdge(0, 1);
        m.addEdge(1, 1); m.addEdge(1, 2);
        m.addEdge(2, 0); m.addEdge(2, 3);
        m.addEdge(3, 2); m.addEdge(3, 3);
        cout << "  max matching = " << m.run() << "\n";
        cout << "  R-side match: ";
        for (int v : m.matchingR()) cout << v << " ";
        cout << "\n";
    }

    // ---- LCA ----
    hr(); cout << "17. LCA (binary lifting)\n";
    {
        // Tree:
        //       0
        //      / \
        //     1   2
        //    / \   \
        //   3   4   5
        //  /
        // 6
        vector<vector<int>> tree(7);
        auto addT = [&](int u, int v){ tree[u].push_back(v); tree[v].push_back(u); };
        addT(0, 1); addT(0, 2); addT(1, 3); addT(1, 4);
        addT(2, 5); addT(3, 6);
        LCA lca(7, tree, 0);
        cout << "  LCA(6,4) = " << lca.query(6, 4)
             << " dist = " << lca.dist(6, 4) << "\n";
        cout << "  LCA(6,5) = " << lca.query(6, 5)
             << " dist = " << lca.dist(6, 5) << "\n";
        cout << "  LCA(3,4) = " << lca.query(3, 4)
             << " dist = " << lca.dist(3, 4) << "\n";
    }

    return 0;
}