Undirected Graph Algorithms
(graph/undirected.hpp)

• View this file on GitHub
• Last update: 2026-08-13 01:41:40+09:00
• Include: #include "graph/undirected.hpp"

Overview

graph/undirected.hpp includes algorithms whose main interpretation is undirected, plus the shortest-path bundle. In that bundle, use the direction-respecting algorithms on graphs built with add_edge; the DAG shortest-path helper is directed-only.

Use this header when edges represent two-way movement or an endpoint constraint where direction should not matter.

Included Headers

Header	Graph orientation	Contents
graph/shortest_path.hpp	Mixed shortest-path bundle	Use BFS, 0-1 BFS, Dijkstra, Bellman-Ford, and Warshall-Floyd on undirected graphs built with add_edge; DAG shortest path is directed-only.
graph/dfs.hpp	Direction-respecting	Iterative DFS forests with parent paths, timestamps, and traversal orders.
graph/lowlink.hpp	Undirected only	Articulation points and bridges.
graph/biconnected_components.hpp	Undirected only	Vertex-biconnected blocks, articulation points, and block incidence.
graph/block_cut_tree.hpp	Undirected only	Block-cut forest and original-vertex-to-node mappings.
graph/two_edge_connected_components.hpp	Undirected only	Two-edge-connected components, bridges, and the contracted bridge forest.
graph/three_edge_connected_components.hpp	Undirected only	Linear-time three-edge-connected vertex decomposition.
graph/st_numbering.hpp	Undirected only	Bipolar numbering between specified source and sink vertices.
graph/kruskal.hpp	Undirected only	Minimum spanning forest.
graph/matrix_tree_theorem.hpp	Undirected or directed	Counts weighted spanning trees with a Laplacian cofactor.
graph/bipartite.hpp	Direction ignored / explicit bipartite sides	Two-colorability, maximum matching, minimum vertex cover, maximum independent set, and minimum edge cover.
graph/general_matching.hpp	Undirected only	Maximum-cardinality matching and minimum edge cover in general undirected graphs.
graph/general_weighted_matching.hpp	Undirected only	Maximum-total-weight matching in general undirected graphs.
graph/maximum_clique.hpp	Direction ignored	Exact maximum clique, maximum independent set, and minimum vertex cover with bitset branch-and-bound.
graph/chromatic_number.hpp	Direction ignored	Exact chromatic number for graphs with at most 20 vertices.
graph/minimum_steiner_tree.hpp	Undirected only	Exact edge- and vertex-weighted minimum Steiner-tree costs and reconstruction for a small terminal set.
graph/replacement_paths.hpp	Undirected positive-weight graphs	Edge- and vertex-failure replacement distances along one fixed shortest path.
graph/namori.hpp	Undirected Namori graph	Ordered cycles and the trees attached to them.
graph/connected_components.hpp	Direction ignored	Weak/ordinary connected components.
graph/count_four_cycles.hpp	Direction ignored	Counts four-cycles globally and for every edge, including parallel-edge choices.
graph/cycle_detection.hpp	Directed and undirected variants	Use find_undirected_cycle(g) for undirected graphs.
graph/enumerate_cliques.hpp	Direction ignored	Enumerates every nonempty clique through a callback.
graph/enumerate_triangles.hpp	Direction ignored	Enumerates every triangle through a callback.
graph/eulerian_trail.hpp	Directed and undirected variants	Use undirected_eulerian_trail(g) for undirected graphs.
graph/grid.hpp	Undirected graph builder	Builds 4/8-neighbor grid graphs.

Complexity

This header is an include bundle and provides no runtime operation by itself. See the included algorithm pages for public interfaces and complexities.

Depends on

• DSU (Disjoint Set Union) (ds/dsu/dsu.hpp)
• Bellman-Ford (graph/bellman_ford.hpp)
• BFS (graph/bfs.hpp)
• Biconnected Components (graph/biconnected_components.hpp)
• Bipartite Graph (graph/bipartite.hpp)
• Block-Cut Tree (graph/block_cut_tree.hpp)
• Chordal Graph Recognition (graph/chordal_graph_recognition.hpp)
• Chromatic Number (graph/chromatic_number.hpp)
• Complement-Graph Connected Components (graph/complement_connected_components.hpp)
• Connected Components (graph/connected_components.hpp)
• Count Four Cycles (graph/count_four_cycles.hpp)
• Cow Game (Difference Constraints) (graph/cow_game.hpp)
• Cycle Detection (graph/cycle_detection.hpp)
• DAG Shortest Path (graph/dag_shortest_path.hpp)
• DFS (graph/dfs.hpp)
• Dijkstra (graph/dijkstra.hpp)
• Enumerate Cliques (graph/enumerate_cliques.hpp)
• Enumerate Triangles (graph/enumerate_triangles.hpp)
• Eulerian Trail (graph/eulerian_trail.hpp)
• General Matching (graph/general_matching.hpp)
• General Weighted Matching (graph/general_weighted_matching.hpp)
• Graph (graph/graph.hpp)
• Grid (graph/grid.hpp)
• K-Shortest Walk (graph/k_shortest_walk.hpp)
• Kruskal (graph/kruskal.hpp)
• LowLink (graph/lowlink.hpp)
• Matrix-Tree Theorem (graph/matrix_tree_theorem.hpp)
• Maximum Clique, Independent Set, and Vertex Cover (graph/maximum_clique.hpp)
• Minimum Steiner Tree (graph/minimum_steiner_tree.hpp)
• Namori Graph Decomposition (graph/namori.hpp)
• Replacement Paths (graph/replacement_paths.hpp)
• Shortest Path (graph/shortest_path.hpp)
• st-Numbering (graph/st_numbering.hpp)
• Three-Edge-Connected Components (graph/three_edge_connected_components.hpp)
• Topological Sort (graph/topological_sort.hpp)
• Two-Edge-Connected Components (graph/two_edge_connected_components.hpp)
• Warshall-Floyd (graph/warshall_floyd.hpp)
• 0-1 BFS (graph/zero_one_bfs.hpp)
• Matrix Linear Algebra (math/matrix/linear_algebra.hpp)
• Dense Matrix (math/matrix/matrix.hpp)

Required by

• Graph All (graph/all.hpp)

Verified with

• verify/graph/cow_game.test.cpp
• verify/graph/graph_algorithms.test.cpp
• verify/graph/range_edge_graph.test.cpp

Code

#ifndef M1UNE_GRAPH_UNDIRECTED_HPP
#define M1UNE_GRAPH_UNDIRECTED_HPP 1

#include "bipartite.hpp"
#include "biconnected_components.hpp"
#include "block_cut_tree.hpp"
#include "chordal_graph_recognition.hpp"
#include "chromatic_number.hpp"
#include "complement_connected_components.hpp"
#include "connected_components.hpp"
#include "count_four_cycles.hpp"
#include "cycle_detection.hpp"
#include "dfs.hpp"
#include "enumerate_cliques.hpp"
#include "enumerate_triangles.hpp"
#include "eulerian_trail.hpp"
#include "general_matching.hpp"
#include "general_weighted_matching.hpp"
#include "graph.hpp"
#include "grid.hpp"
#include "kruskal.hpp"
#include "lowlink.hpp"
#include "maximum_clique.hpp"
#include "matrix_tree_theorem.hpp"
#include "minimum_steiner_tree.hpp"
#include "namori.hpp"
#include "replacement_paths.hpp"
#include "shortest_path.hpp"
#include "st_numbering.hpp"
#include "three_edge_connected_components.hpp"
#include "two_edge_connected_components.hpp"

#endif  // M1UNE_GRAPH_UNDIRECTED_HPP

#line 1 "graph/undirected.hpp"



#line 1 "graph/bipartite.hpp"



#include <algorithm>
#include <cassert>
#include <cstddef>
#include <cstdint>
#include <limits>
#include <optional>
#include <queue>
#include <utility>
#include <vector>

#line 1 "graph/graph.hpp"



#include <array>
#line 8 "graph/graph.hpp"

namespace m1une {
namespace graph {

template <class T = int>
struct Edge {
    using cost_type = T;

    int from;
    int to;
    T cost;
    int id;
    bool alive;

    Edge() : from(-1), to(-1), cost(T()), id(-1), alive(true) {}
    Edge(int from_, int to_, T cost_ = T(1), int id_ = -1, bool alive_ = true)
        : from(from_), to(to_), cost(cost_), id(id_), alive(alive_) {}

    int other(int v) const {
        assert(v == from || v == to);
        return from ^ to ^ v;
    }
};

template <class T = int>
struct Graph {
    using edge_type = Edge<T>;
    using cost_type = T;

   private:
    struct EdgePositions {
        std::array<std::pair<int, int>, 2> value{};
        int size = 0;

        void push_back(std::pair<int, int> position) {
            assert(size < 2);
            value[size++] = position;
        }
    };

    int _n;
    int _edge_count;
    std::vector<std::vector<edge_type>> _g;
    std::vector<EdgePositions> _edge_positions;

   public:
    Graph() : _n(0), _edge_count(0) {}
    explicit Graph(int n) : _n(n), _edge_count(0), _g(n) {
        assert(0 <= n);
    }

    int size() const {
        return _n;
    }

    bool empty() const {
        return _n == 0;
    }

    int edge_count() const {
        return _edge_count;
    }

    int add_vertex() {
        _g.emplace_back();
        return _n++;
    }

    int add_directed_edge(int from, int to, T cost = T(1)) {
        assert(0 <= from && from < _n);
        assert(0 <= to && to < _n);
        int id = _edge_count++;
        int idx = int(_g[from].size());
        _g[from].push_back(edge_type(from, to, cost, id));
        _edge_positions.emplace_back();
        _edge_positions.back().push_back({from, idx});
        return id;
    }

    int add_edge(int u, int v, T cost = T(1)) {
        assert(0 <= u && u < _n);
        assert(0 <= v && v < _n);
        int id = _edge_count++;
        int u_idx = int(_g[u].size());
        _g[u].push_back(edge_type(u, v, cost, id));
        int v_idx = int(_g[v].size());
        _g[v].push_back(edge_type(v, u, cost, id));
        _edge_positions.emplace_back();
        _edge_positions.back().push_back({u, u_idx});
        _edge_positions.back().push_back({v, v_idx});
        return id;
    }

    void set_edge_alive(int id, bool alive) {
        assert(0 <= id && id < _edge_count);
        for (int i = 0; i < _edge_positions[id].size; ++i) {
            auto [v, idx] = _edge_positions[id].value[i];
            _g[v][idx].alive = alive;
        }
    }

    void erase_edge(int id) {
        set_edge_alive(id, false);
    }

    void revive_edge(int id) {
        set_edge_alive(id, true);
    }

    bool is_edge_alive(int id) const {
        assert(0 <= id && id < _edge_count);
        assert(_edge_positions[id].size != 0);
        auto [v, idx] = _edge_positions[id].value[0];
        return _g[v][idx].alive;
    }

    const std::vector<edge_type>& operator[](int v) const {
        assert(0 <= v && v < _n);
        return _g[v];
    }

    std::vector<edge_type>& operator[](int v) {
        assert(0 <= v && v < _n);
        return _g[v];
    }

    const std::vector<std::vector<edge_type>>& adjacency() const {
        return _g;
    }

    std::vector<std::vector<edge_type>>& adjacency() {
        return _g;
    }

    std::vector<edge_type> edges(bool include_inactive = false) const {
        std::vector<edge_type> result;
        result.reserve(_edge_count);
        std::vector<char> used(_edge_count, false);
        for (int v = 0; v < _n; v++) {
            for (const auto& e : _g[v]) {
                if (!include_inactive && !e.alive) continue;
                if (0 <= e.id && e.id < _edge_count) {
                    if (used[e.id]) continue;
                    used[e.id] = true;
                }
                result.push_back(e);
            }
        }
        return result;
    }

    Graph reversed() const {
        Graph result(_n);
        result._edge_count = _edge_count;
        result._edge_positions.assign(_edge_count, {});
        for (int v = 0; v < _n; v++) {
            for (const auto& e : _g[v]) {
                int idx = int(result._g[e.to].size());
                result._g[e.to].push_back(edge_type(e.to, e.from, e.cost, e.id, e.alive));
                if (0 <= e.id && e.id < _edge_count) result._edge_positions[e.id].push_back({e.to, idx});
            }
        }
        return result;
    }
};

}  // namespace graph
}  // namespace m1une


#line 15 "graph/bipartite.hpp"

namespace m1une {
namespace graph {

struct BipartiteResult {
    bool is_bipartite;
    std::vector<int> color;
    std::vector<int> left_vertices;
    std::vector<int> right_vertices;
    std::vector<int> left_id;
    std::vector<int> right_id;
};

template <class T>
BipartiteResult bipartite(const Graph<T>& g) {
    int n = g.size();
    BipartiteResult result;
    result.is_bipartite = true;
    result.color.assign(n, -1);
    result.left_id.assign(n, -1);
    result.right_id.assign(n, -1);

    std::vector<std::vector<int>> adjacency(n);
    for (const auto& e : g.edges()) {
        adjacency[e.from].push_back(e.to);
        adjacency[e.to].push_back(e.from);
    }

    std::queue<int> que;
    for (int s = 0; s < n; s++) {
        if (result.color[s] != -1) continue;
        result.color[s] = 0;
        que.push(s);
        while (!que.empty()) {
            int v = que.front();
            que.pop();
            for (int to : adjacency[v]) {
                if (result.color[to] == -1) {
                    result.color[to] = result.color[v] ^ 1;
                    que.push(to);
                } else if (result.color[to] == result.color[v]) {
                    result.is_bipartite = false;
                    return result;
                }
            }
        }
    }

    for (int v = 0; v < n; v++) {
        if (result.color[v] == 0) {
            result.left_id[v] = int(result.left_vertices.size());
            result.left_vertices.push_back(v);
        } else {
            result.right_id[v] = int(result.right_vertices.size());
            result.right_vertices.push_back(v);
        }
    }

    return result;
}

template <class T>
bool is_bipartite(const Graph<T>& g) {
    return bipartite(g).is_bipartite;
}

struct BipartiteVertexSet {
    std::vector<int> left;
    std::vector<int> right;

    int size() const {
        return int(left.size() + right.size());
    }
};

struct BipartiteMatching {
    struct Edge {
        int left;
        int right;
        int id;
        bool alive;
    };

    struct Pair {
        int left;
        int right;
        int edge_id;
    };

   private:
    int _left_size;
    int _right_size;
    std::vector<Edge> _edges;
    std::vector<std::vector<int>> _adj;
    std::vector<std::vector<int>> _radj;
    std::vector<int> _left_match;
    std::vector<int> _right_match;
    std::vector<int> _left_match_edge;
    std::vector<int> _right_match_edge;
    bool _calculated;

    void invalidate() {
        _calculated = false;
    }

    void ensure_matching() {
        if (!_calculated) max_matching();
    }

   public:
    BipartiteMatching() : BipartiteMatching(0, 0) {}

    BipartiteMatching(int left_size, int right_size)
        : _left_size(left_size),
          _right_size(right_size),
          _adj(left_size),
          _radj(right_size),
          _left_match(left_size, -1),
          _right_match(right_size, -1),
          _left_match_edge(left_size, -1),
          _right_match_edge(right_size, -1),
          _calculated(false) {
        assert(0 <= left_size);
        assert(0 <= right_size);
    }

    int left_size() const {
        return _left_size;
    }

    int right_size() const {
        return _right_size;
    }

    int edge_count() const {
        return int(_edges.size());
    }

    int add_edge(int left, int right) {
        assert(0 <= left && left < _left_size);
        assert(0 <= right && right < _right_size);
        int id = int(_edges.size());
        _edges.push_back(Edge{left, right, id, true});
        _adj[left].push_back(id);
        _radj[right].push_back(id);
        invalidate();
        return id;
    }

    Edge get_edge(int i) const {
        assert(0 <= i && i < int(_edges.size()));
        return _edges[i];
    }

    std::vector<Edge> edges(bool include_inactive = false) const {
        std::vector<Edge> result;
        result.reserve(_edges.size());
        for (const auto& e : _edges) {
            if (include_inactive || e.alive) result.push_back(e);
        }
        return result;
    }

    void set_edge_alive(int id, bool alive) {
        assert(0 <= id && id < int(_edges.size()));
        _edges[id].alive = alive;
        invalidate();
    }

    void erase_edge(int id) {
        set_edge_alive(id, false);
    }

    void revive_edge(int id) {
        set_edge_alive(id, true);
    }

    bool is_edge_alive(int id) const {
        assert(0 <= id && id < int(_edges.size()));
        return _edges[id].alive;
    }

    int max_matching() {
        _left_match.assign(_left_size, -1);
        _right_match.assign(_right_size, -1);
        _left_match_edge.assign(_left_size, -1);
        _right_match_edge.assign(_right_size, -1);

        std::vector<int> dist(_left_size);
        auto bfs = [&]() -> bool {
            std::queue<int> que;
            bool found = false;
            for (int l = 0; l < _left_size; l++) {
                if (_left_match[l] == -1) {
                    dist[l] = 0;
                    que.push(l);
                } else {
                    dist[l] = -1;
                }
            }

            while (!que.empty()) {
                int l = que.front();
                que.pop();
                for (int id : _adj[l]) {
                    const auto& e = _edges[id];
                    if (!e.alive) continue;
                    int next_left = _right_match[e.right];
                    if (next_left == -1) {
                        found = true;
                    } else if (dist[next_left] == -1) {
                        dist[next_left] = dist[l] + 1;
                        que.push(next_left);
                    }
                }
            }
            return found;
        };

        auto dfs = [&](auto self, int l) -> bool {
            for (int id : _adj[l]) {
                const auto& e = _edges[id];
                if (!e.alive) continue;
                int next_left = _right_match[e.right];
                if (next_left != -1 && (dist[next_left] != dist[l] + 1 || !self(self, next_left))) {
                    continue;
                }
                _left_match[l] = e.right;
                _right_match[e.right] = l;
                _left_match_edge[l] = id;
                _right_match_edge[e.right] = id;
                return true;
            }
            dist[l] = -1;
            return false;
        };

        int result = 0;
        while (bfs()) {
            for (int l = 0; l < _left_size; l++) {
                if (_left_match[l] == -1 && dfs(dfs, l)) result++;
            }
        }

        _calculated = true;
        return result;
    }

    int matching_size() {
        ensure_matching();
        int result = 0;
        for (int right : _left_match) {
            if (right != -1) result++;
        }
        return result;
    }

    std::vector<int> left_match() {
        ensure_matching();
        return _left_match;
    }

    std::vector<int> right_match() {
        ensure_matching();
        return _right_match;
    }

    std::vector<Pair> matching() {
        ensure_matching();
        std::vector<Pair> result;
        for (int l = 0; l < _left_size; l++) {
            if (_left_match[l] != -1) result.push_back(Pair{l, _left_match[l], _left_match_edge[l]});
        }
        return result;
    }

    BipartiteVertexSet minimum_vertex_cover() {
        ensure_matching();

        std::vector<char> visited_left(_left_size, false), visited_right(_right_size, false);
        std::queue<int> que;
        for (int l = 0; l < _left_size; l++) {
            if (_left_match[l] == -1) {
                visited_left[l] = true;
                que.push(l);
            }
        }

        while (!que.empty()) {
            int l = que.front();
            que.pop();
            for (int id : _adj[l]) {
                const auto& e = _edges[id];
                if (!e.alive || _left_match_edge[l] == id || visited_right[e.right]) continue;
                visited_right[e.right] = true;
                int next_left = _right_match[e.right];
                if (next_left != -1 && !visited_left[next_left]) {
                    visited_left[next_left] = true;
                    que.push(next_left);
                }
            }
        }

        BipartiteVertexSet result;
        for (int l = 0; l < _left_size; l++) {
            if (!visited_left[l]) result.left.push_back(l);
        }
        for (int r = 0; r < _right_size; r++) {
            if (visited_right[r]) result.right.push_back(r);
        }
        return result;
    }

    BipartiteVertexSet maximum_independent_set() {
        auto cover = minimum_vertex_cover();
        std::vector<char> in_left_cover(_left_size, false), in_right_cover(_right_size, false);
        for (int l : cover.left) in_left_cover[l] = true;
        for (int r : cover.right) in_right_cover[r] = true;

        BipartiteVertexSet result;
        for (int l = 0; l < _left_size; l++) {
            if (!in_left_cover[l]) result.left.push_back(l);
        }
        for (int r = 0; r < _right_size; r++) {
            if (!in_right_cover[r]) result.right.push_back(r);
        }
        return result;
    }

    std::optional<std::vector<int>> minimum_edge_cover() {
        ensure_matching();

        std::vector<int> result;
        std::vector<char> covered_left(_left_size, false), covered_right(_right_size, false);
        std::vector<char> used_edge(_edges.size(), false);

        auto use_edge = [&](int id) {
            if (used_edge[id]) return;
            used_edge[id] = true;
            result.push_back(id);
            covered_left[_edges[id].left] = true;
            covered_right[_edges[id].right] = true;
        };

        for (int l = 0; l < _left_size; l++) {
            if (_left_match_edge[l] != -1) use_edge(_left_match_edge[l]);
        }

        for (int l = 0; l < _left_size; l++) {
            if (covered_left[l]) continue;
            int id = -1;
            for (int edge_id : _adj[l]) {
                if (_edges[edge_id].alive) {
                    id = edge_id;
                    break;
                }
            }
            if (id == -1) return std::nullopt;
            use_edge(id);
        }

        for (int r = 0; r < _right_size; r++) {
            if (covered_right[r]) continue;
            int id = -1;
            for (int edge_id : _radj[r]) {
                if (_edges[edge_id].alive) {
                    id = edge_id;
                    break;
                }
            }
            if (id == -1) return std::nullopt;
            use_edge(id);
        }

        return result;
    }
};

struct BipartiteMatchingGraph {
    BipartiteResult parts;
    BipartiteMatching matching;
    std::vector<int> original_edge_id;

    int left_vertex(int left) const {
        assert(0 <= left && left < int(parts.left_vertices.size()));
        return parts.left_vertices[left];
    }

    int right_vertex(int right) const {
        assert(0 <= right && right < int(parts.right_vertices.size()));
        return parts.right_vertices[right];
    }

    int original_edge(int edge_id) const {
        assert(0 <= edge_id && edge_id < int(original_edge_id.size()));
        return original_edge_id[edge_id];
    }
};

template <class T>
std::optional<BipartiteMatchingGraph> make_bipartite_matching(const Graph<T>& g) {
    auto parts = bipartite(g);
    if (!parts.is_bipartite) return std::nullopt;

    BipartiteMatchingGraph result;
    result.parts = parts;
    result.matching = BipartiteMatching(int(parts.left_vertices.size()), int(parts.right_vertices.size()));

    for (const auto& e : g.edges()) {
        int left, right;
        if (parts.color[e.from] == 0) {
            left = parts.left_id[e.from];
            right = parts.right_id[e.to];
        } else {
            left = parts.left_id[e.to];
            right = parts.right_id[e.from];
        }
        int id = result.matching.add_edge(left, right);
        if (int(result.original_edge_id.size()) <= id) result.original_edge_id.resize(id + 1);
        result.original_edge_id[id] = e.id;
    }

    return result;
}

struct BipartiteEdgeColoringResult {
    int color_count;
    std::vector<int> color;
};

namespace detail {

struct BipartiteEdgeColoringGroups {
    int count;
    std::vector<int> group;
};

inline BipartiteEdgeColoringGroups group_vertices(
    const std::vector<int>& degree,
    int maximum_degree
) {
    BipartiteEdgeColoringGroups result;
    result.count = 0;
    result.group.assign(degree.size(), -1);
    int current_degree = 0;
    for (int vertex = 0; vertex < int(degree.size()); vertex++) {
        if (degree[vertex] == 0) continue;
        if (result.count == 0 || current_degree + degree[vertex] > maximum_degree) {
            result.count++;
            current_degree = 0;
        }
        result.group[vertex] = result.count - 1;
        current_degree += degree[vertex];
    }
    return result;
}

class BipartiteEdgeColoringSolver {
   private:
    int _side_size;
    int _original_edge_count;
    std::vector<int> _left;
    std::vector<int> _right;
    std::vector<int> _color;
    std::vector<int> _used_stamp;
    int _stamp;

    int other_endpoint(int vertex, int edge) const {
        if (vertex < _side_size) return _side_size + _right[edge];
        return _left[edge];
    }

    std::vector<int> perfect_matching(const std::vector<int>& edge_ids) const {
        std::vector<std::vector<int>> adjacency(_side_size);
        for (int edge : edge_ids) adjacency[_left[edge]].push_back(edge);

        std::vector<int> right_match(_side_size, -1);
        std::vector<int> left_match_edge(_side_size, -1);
        int matching_size = 0;
        for (int left = 0; left < _side_size; left++) {
            for (int edge : adjacency[left]) {
                int right = _right[edge];
                if (right_match[right] != -1) continue;
                right_match[right] = left;
                left_match_edge[left] = edge;
                matching_size++;
                break;
            }
        }

        std::vector<int> distance(_side_size);
        std::vector<int> next_edge(_side_size);
        std::vector<int> left_stack;
        std::vector<int> path_edges;
        left_stack.reserve(_side_size);
        path_edges.reserve(_side_size);

        while (matching_size < _side_size) {
            std::queue<int> queue;
            std::fill(distance.begin(), distance.end(), -1);
            for (int left = 0; left < _side_size; left++) {
                if (left_match_edge[left] != -1) continue;
                distance[left] = 0;
                queue.push(left);
            }

            bool reachable_free_right = false;
            while (!queue.empty()) {
                int left = queue.front();
                queue.pop();
                for (int edge : adjacency[left]) {
                    int next_left = right_match[_right[edge]];
                    if (next_left == -1) {
                        reachable_free_right = true;
                    } else if (distance[next_left] == -1) {
                        distance[next_left] = distance[left] + 1;
                        queue.push(next_left);
                    }
                }
            }
            assert(reachable_free_right);

            std::fill(next_edge.begin(), next_edge.end(), 0);
            int augmented = 0;
            for (int root = 0; root < _side_size; root++) {
                if (left_match_edge[root] != -1 || distance[root] == -1) continue;
                left_stack.clear();
                path_edges.clear();
                left_stack.push_back(root);
                bool found = false;

                while (!left_stack.empty() && !found) {
                    int left = left_stack.back();
                    bool advanced = false;
                    while (next_edge[left] < int(adjacency[left].size())) {
                        int edge = adjacency[left][next_edge[left]++];
                        int right = _right[edge];
                        int next_left = right_match[right];
                        if (next_left == -1) {
                            left_match_edge[left] = edge;
                            right_match[right] = left;
                            for (int index = int(path_edges.size()) - 1; index >= 0; index--) {
                                int path_edge = path_edges[index];
                                int path_left = left_stack[index];
                                left_match_edge[path_left] = path_edge;
                                right_match[_right[path_edge]] = path_left;
                            }
                            found = true;
                            break;
                        }
                        if (distance[next_left] != distance[left] + 1) continue;
                        path_edges.push_back(edge);
                        left_stack.push_back(next_left);
                        advanced = true;
                        break;
                    }
                    if (found || advanced) continue;
                    distance[left] = -1;
                    left_stack.pop_back();
                    if (path_edges.size() == left_stack.size() && !path_edges.empty()) {
                        path_edges.pop_back();
                    }
                }
                if (found) augmented++;
            }
            assert(augmented > 0);
            matching_size += augmented;
        }

        return left_match_edge;
    }

    std::pair<std::vector<int>, std::vector<int>> split_even(
        const std::vector<int>& edge_ids
    ) {
        std::vector<std::vector<int>> incidence(std::size_t(2) * _side_size);
        for (int edge : edge_ids) {
            incidence[_left[edge]].push_back(edge);
            incidence[_side_size + _right[edge]].push_back(edge);
        }

        _stamp++;
        assert(_stamp > 0);
        std::vector<int> next_edge(std::size_t(2) * _side_size, 0);
        std::vector<int> first;
        std::vector<int> second;
        first.reserve(edge_ids.size() / 2);
        second.reserve(edge_ids.size() / 2);

        for (int start = 0; start < 2 * _side_size; start++) {
            while (true) {
                while (next_edge[start] < int(incidence[start].size()) &&
                       _used_stamp[incidence[start][next_edge[start]]] == _stamp) {
                    next_edge[start]++;
                }
                if (next_edge[start] == int(incidence[start].size())) break;

                int vertex = start;
                bool parity = false;
                do {
                    while (next_edge[vertex] < int(incidence[vertex].size()) &&
                           _used_stamp[incidence[vertex][next_edge[vertex]]] == _stamp) {
                        next_edge[vertex]++;
                    }
                    assert(next_edge[vertex] < int(incidence[vertex].size()));
                    int edge = incidence[vertex][next_edge[vertex]++];
                    _used_stamp[edge] = _stamp;
                    if (!parity) {
                        first.push_back(edge);
                    } else {
                        second.push_back(edge);
                    }
                    parity = !parity;
                    vertex = other_endpoint(vertex, edge);
                } while (vertex != start);
                assert(!parity);
            }
        }
        assert(first.size() == second.size());
        return {std::move(first), std::move(second)};
    }

    void color_regular(const std::vector<int>& edge_ids, int degree, int offset) {
        assert(std::size_t(_side_size) * std::size_t(degree) == edge_ids.size());
        if (degree == 0) return;
        if (degree == 1) {
            for (int edge : edge_ids) {
                if (edge < _original_edge_count) _color[edge] = offset;
            }
            return;
        }

        if (degree % 2 == 1) {
            std::vector<int> matching = perfect_matching(edge_ids);
            _stamp++;
            assert(_stamp > 0);
            for (int edge : matching) {
                _used_stamp[edge] = _stamp;
                if (edge < _original_edge_count) _color[edge] = offset;
            }
            std::vector<int> remaining;
            remaining.reserve(edge_ids.size() - matching.size());
            for (int edge : edge_ids) {
                if (_used_stamp[edge] != _stamp) remaining.push_back(edge);
            }
            color_regular(remaining, degree - 1, offset + 1);
            return;
        }

        auto [first, second] = split_even(edge_ids);
        color_regular(first, degree / 2, offset);
        color_regular(second, degree / 2, offset + degree / 2);
    }

   public:
    BipartiteEdgeColoringSolver(
        int side_size,
        int original_edge_count,
        std::vector<int> left,
        std::vector<int> right
    )
        : _side_size(side_size),
          _original_edge_count(original_edge_count),
          _left(std::move(left)),
          _right(std::move(right)),
          _color(original_edge_count, -1),
          _used_stamp(_left.size(), 0),
          _stamp(0) {}

    std::vector<int> solve(int degree) {
        std::vector<int> edge_ids(_left.size());
        for (int edge = 0; edge < int(edge_ids.size()); edge++) edge_ids[edge] = edge;
        color_regular(edge_ids, degree, 0);
        for (int color : _color) assert(0 <= color && color < degree);
        return _color;
    }
};

}  // namespace detail

// Returns an optimal edge coloring of a bipartite multigraph.
inline BipartiteEdgeColoringResult bipartite_edge_coloring(
    int left_size,
    int right_size,
    const std::vector<std::pair<int, int>>& edges
) {
    assert(left_size >= 0);
    assert(right_size >= 0);
    assert(edges.size() <= std::size_t(std::numeric_limits<int>::max()));

    std::vector<int> left_degree(left_size, 0);
    std::vector<int> right_degree(right_size, 0);
    int maximum_degree = 0;
    for (auto [left, right] : edges) {
        assert(0 <= left && left < left_size);
        assert(0 <= right && right < right_size);
        left_degree[left]++;
        right_degree[right]++;
        maximum_degree = std::max(maximum_degree, left_degree[left]);
        maximum_degree = std::max(maximum_degree, right_degree[right]);
    }

    BipartiteEdgeColoringResult result;
    result.color_count = maximum_degree;
    if (edges.empty()) return result;

    detail::BipartiteEdgeColoringGroups left_groups =
        detail::group_vertices(left_degree, maximum_degree);
    detail::BipartiteEdgeColoringGroups right_groups =
        detail::group_vertices(right_degree, maximum_degree);
    int side_size = std::max(left_groups.count, right_groups.count);

    std::vector<int> contracted_left;
    std::vector<int> contracted_right;
    contracted_left.reserve(std::size_t(3) * edges.size());
    contracted_right.reserve(std::size_t(3) * edges.size());
    std::vector<int> contracted_left_degree(side_size, 0);
    std::vector<int> contracted_right_degree(side_size, 0);
    for (auto [left, right] : edges) {
        int contracted_left_vertex = left_groups.group[left];
        int contracted_right_vertex = right_groups.group[right];
        contracted_left.push_back(contracted_left_vertex);
        contracted_right.push_back(contracted_right_vertex);
        contracted_left_degree[contracted_left_vertex]++;
        contracted_right_degree[contracted_right_vertex]++;
    }

    int left = 0;
    int right = 0;
    while (true) {
        while (left < side_size && contracted_left_degree[left] == maximum_degree) left++;
        while (right < side_size && contracted_right_degree[right] == maximum_degree) right++;
        if (left == side_size || right == side_size) break;
        contracted_left.push_back(left);
        contracted_right.push_back(right);
        contracted_left_degree[left]++;
        contracted_right_degree[right]++;
    }
    assert(left == side_size && right == side_size);
    assert(contracted_left.size() == std::size_t(side_size) * std::size_t(maximum_degree));

    detail::BipartiteEdgeColoringSolver solver(
        side_size,
        int(edges.size()),
        std::move(contracted_left),
        std::move(contracted_right)
    );
    result.color = solver.solve(maximum_degree);
    return result;
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/biconnected_components.hpp"



#line 6 "graph/biconnected_components.hpp"

#line 8 "graph/biconnected_components.hpp"

namespace m1une {
namespace graph {

struct BiconnectedComponentsResult {
    std::vector<std::vector<int>> components;
    std::vector<std::vector<int>> edge_components;
    std::vector<int> component_of_edge;
    std::vector<std::vector<int>> vertex_components;
    std::vector<int> articulation;
    std::vector<int> ord;
    std::vector<int> low;

    int component_count() const {
        return int(components.size());
    }

    bool is_articulation(int vertex) const {
        assert(0 <= vertex && vertex < int(vertex_components.size()));
        return vertex_components[vertex].size() >= 2;
    }
};

// Decomposes an undirected graph into maximal vertex-biconnected blocks.
// Every active edge belongs to exactly one block. Isolated vertices form
// singleton blocks, and articulation vertices occur in multiple blocks.
template <class T>
BiconnectedComponentsResult biconnected_components(const Graph<T>& graph) {
    const int n = graph.size();
    const int edge_count = graph.edge_count();

    BiconnectedComponentsResult result;
    result.component_of_edge.assign(edge_count, -1);
    result.vertex_components.assign(n, {});
    result.ord.assign(n, -1);
    result.low.assign(n, -1);

    std::vector<int> edge_from(edge_count, -1);
    std::vector<int> edge_to(edge_count, -1);
    std::vector<int> incidence_count(edge_count, 0);
    std::vector<int> alive_degree(n, 0);
    for (int vertex = 0; vertex < n; vertex++) {
        for (const Edge<T>& edge : graph[vertex]) {
            if (!edge.alive) continue;
            assert(0 <= edge.id && edge.id < edge_count);
            alive_degree[vertex]++;
            if (incidence_count[edge.id] == 0) {
                edge_from[edge.id] = edge.from;
                edge_to[edge.id] = edge.to;
            }
            incidence_count[edge.id]++;
        }
    }
#ifndef NDEBUG
    for (int edge_id = 0; edge_id < edge_count; edge_id++) {
        if (incidence_count[edge_id] == 0) continue;
        assert(incidence_count[edge_id] == 2);
        assert(edge_from[edge_id] != edge_to[edge_id]);
    }
#endif

    std::vector<int> parent(n, -1);
    std::vector<int> parent_edge(n, -1);
    std::vector<int> next_edge(n, 0);
    std::vector<int> dfs_stack;
    std::vector<int> edge_stack;
    std::vector<int> vertex_mark(n, -1);
    int timer = 0;

    auto add_singleton = [&](int vertex) {
        const int component = result.component_count();
        result.components.push_back(std::vector<int>(1, vertex));
        result.edge_components.emplace_back();
        result.vertex_components[vertex].push_back(component);
    };

    auto extract_component = [&](int stopping_edge) {
        const int component = result.component_count();
        result.components.emplace_back();
        result.edge_components.emplace_back();
        std::vector<int>& vertices = result.components.back();
        std::vector<int>& edges = result.edge_components.back();

        while (true) {
            assert(!edge_stack.empty());
            const int edge_id = edge_stack.back();
            edge_stack.pop_back();
            edges.push_back(edge_id);
            result.component_of_edge[edge_id] = component;

            const int endpoints[2] = {edge_from[edge_id], edge_to[edge_id]};
            for (int vertex : endpoints) {
                if (vertex_mark[vertex] == component) continue;
                vertex_mark[vertex] = component;
                vertices.push_back(vertex);
            }
            if (edge_id == stopping_edge) break;
        }
        for (int vertex : vertices) {
            result.vertex_components[vertex].push_back(component);
        }
    };

    for (int root = 0; root < n; root++) {
        if (result.ord[root] != -1) continue;
        if (alive_degree[root] == 0) {
            result.ord[root] = result.low[root] = timer++;
            add_singleton(root);
            continue;
        }

        result.ord[root] = result.low[root] = timer++;
        dfs_stack.push_back(root);
        while (!dfs_stack.empty()) {
            const int vertex = dfs_stack.back();
            if (next_edge[vertex] < int(graph[vertex].size())) {
                const Edge<T>& edge = graph[vertex][next_edge[vertex]++];
                if (!edge.alive || edge.id == parent_edge[vertex]) continue;
                const int to = edge.to;
                if (result.ord[to] == -1) {
                    parent[to] = vertex;
                    parent_edge[to] = edge.id;
                    edge_stack.push_back(edge.id);
                    result.ord[to] = result.low[to] = timer++;
                    dfs_stack.push_back(to);
                } else if (result.ord[to] < result.ord[vertex]) {
                    edge_stack.push_back(edge.id);
                    if (result.ord[to] < result.low[vertex]) {
                        result.low[vertex] = result.ord[to];
                    }
                }
                continue;
            }

            dfs_stack.pop_back();
            const int parent_vertex = parent[vertex];
            if (parent_vertex == -1) {
                assert(edge_stack.empty());
                continue;
            }
            if (result.low[vertex] < result.low[parent_vertex]) {
                result.low[parent_vertex] = result.low[vertex];
            }
            if (result.ord[parent_vertex] <= result.low[vertex]) {
                extract_component(parent_edge[vertex]);
            }
        }
    }

    for (int vertex = 0; vertex < n; vertex++) {
        if (result.is_articulation(vertex)) result.articulation.push_back(vertex);
    }
    return result;
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/block_cut_tree.hpp"



#line 6 "graph/block_cut_tree.hpp"

#line 8 "graph/block_cut_tree.hpp"

namespace m1une {
namespace graph {

struct BlockCutTreeResult {
    std::vector<std::vector<int>> forest;
    std::vector<int> node_of_block;
    std::vector<int> node_of_articulation;
    std::vector<int> node_of_vertex;
    std::vector<int> block_of_node;
    std::vector<int> articulation_of_node;

    int node_count() const {
        return int(forest.size());
    }

    int block_count() const {
        return int(node_of_block.size());
    }

    bool is_block_node(int node) const {
        assert(0 <= node && node < node_count());
        return block_of_node[node] != -1;
    }

    bool is_articulation_node(int node) const {
        assert(0 <= node && node < node_count());
        return articulation_of_node[node] != -1;
    }
};

// Builds the block-cut forest of a biconnected-components decomposition.
// Block nodes have IDs [0, block_count); articulation nodes follow them.
inline BlockCutTreeResult block_cut_tree(
    const BiconnectedComponentsResult& biconnected
) {
    const int vertex_count = int(biconnected.vertex_components.size());
    const int block_count = biconnected.component_count();

    BlockCutTreeResult result;
    result.node_of_block.resize(block_count);
    result.node_of_articulation.assign(vertex_count, -1);
    result.node_of_vertex.assign(vertex_count, -1);
    result.forest.resize(block_count);
    result.block_of_node.resize(block_count);
    result.articulation_of_node.assign(block_count, -1);
    for (int block = 0; block < block_count; block++) {
        result.node_of_block[block] = block;
        result.block_of_node[block] = block;
    }

    for (int vertex = 0; vertex < vertex_count; vertex++) {
        const std::vector<int>& blocks = biconnected.vertex_components[vertex];
        assert(!blocks.empty());
        if (blocks.size() == 1) {
            assert(0 <= blocks[0] && blocks[0] < block_count);
            result.node_of_vertex[vertex] = result.node_of_block[blocks[0]];
            continue;
        }

        const int node = result.node_count();
        result.node_of_articulation[vertex] = node;
        result.node_of_vertex[vertex] = node;
        result.forest.emplace_back();
        result.block_of_node.push_back(-1);
        result.articulation_of_node.push_back(vertex);
        for (int block : blocks) {
            assert(0 <= block && block < block_count);
            const int block_node = result.node_of_block[block];
            result.forest[node].push_back(block_node);
            result.forest[block_node].push_back(node);
        }
    }
    return result;
}

template <class T>
BlockCutTreeResult block_cut_tree(const Graph<T>& graph) {
    return block_cut_tree(biconnected_components(graph));
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/chordal_graph_recognition.hpp"



#line 9 "graph/chordal_graph_recognition.hpp"

#line 11 "graph/chordal_graph_recognition.hpp"

namespace m1une {
namespace graph {

struct ChordalGraphResult {
    bool is_chordal;
    std::vector<int> perfect_elimination_order;
    std::vector<int> induced_cycle;
};

namespace internal {

class MaximumCardinalitySearch {
    std::vector<int> _head;
    std::vector<int> _next;
    std::vector<int> _previous;
    std::vector<int> _weight;

    void erase(int vertex) {
        const int weight = _weight[vertex];
        if (_previous[vertex] == -1) {
            _head[weight] = _next[vertex];
        } else {
            _next[_previous[vertex]] = _next[vertex];
        }
        if (_next[vertex] != -1) _previous[_next[vertex]] = _previous[vertex];
    }

    void insert(int vertex) {
        const int weight = _weight[vertex];
        _previous[vertex] = -1;
        _next[vertex] = _head[weight];
        if (_head[weight] != -1) _previous[_head[weight]] = vertex;
        _head[weight] = vertex;
    }

   public:
    explicit MaximumCardinalitySearch(int size)
        : _head(size + 1, -1),
          _next(size, -1),
          _previous(size, -1),
          _weight(size, 0) {
        for (int vertex = 0; vertex < size; vertex++) insert(vertex);
    }

    std::vector<int> run(const std::vector<std::vector<int>>& adjacency) {
        const int size = int(adjacency.size());
        std::vector<int> order;
        order.reserve(size);
        std::vector<char> selected(size, false);
        std::vector<int> seen_neighbor(size, -1);
        int maximum_weight = 0;

        while (int(order.size()) < size) {
            while (_head[maximum_weight] == -1) maximum_weight--;
            const int vertex = _head[maximum_weight];
            erase(vertex);
            selected[vertex] = true;
            order.push_back(vertex);

            for (int to : adjacency[vertex]) {
                if (to == vertex || selected[to] || seen_neighbor[to] == vertex) continue;
                seen_neighbor[to] = vertex;
                erase(to);
                _weight[to]++;
                insert(to);
                maximum_weight = std::max(maximum_weight, _weight[to]);
            }
        }
        return order;
    }
};

inline std::vector<int> chordless_cycle(
    const std::vector<std::vector<int>>& adjacency, int vertex, int first,
    int second
) {
    const int size = int(adjacency.size());
    std::vector<char> forbidden(size, false);
    for (int to : adjacency[vertex]) forbidden[to] = true;
    forbidden[vertex] = true;
    forbidden[first] = false;
    forbidden[second] = false;

    std::vector<int> parent(size, -1);
    std::queue<int> queue;
    parent[first] = first;
    queue.push(first);
    while (!queue.empty() && parent[second] == -1) {
        const int current = queue.front();
        queue.pop();
        for (int to : adjacency[current]) {
            if (forbidden[to] || parent[to] != -1) continue;
            parent[to] = current;
            queue.push(to);
        }
    }
    assert(parent[second] != -1);

    std::vector<int> path;
    for (int current = second; current != first; current = parent[current]) {
        path.push_back(current);
    }
    path.push_back(first);
    std::reverse(path.begin(), path.end());

    std::vector<int> cycle;
    cycle.reserve(path.size() + 1);
    cycle.push_back(vertex);
    cycle.insert(cycle.end(), path.begin(), path.end());
    return cycle;
}

}  // namespace internal

// Recognizes a chordal graph. On success, returns a perfect elimination
// ordering; on failure, returns an induced cycle of length at least four.
template <class T>
ChordalGraphResult chordal_graph_recognition(const Graph<T>& graph) {
    const int size = graph.size();
    std::vector<std::vector<int>> adjacency(size);
    for (const Edge<T>& edge : graph.edges()) {
        if (edge.from == edge.to) continue;
        adjacency[edge.from].push_back(edge.to);
        adjacency[edge.to].push_back(edge.from);
    }

    std::vector<int> order = internal::MaximumCardinalitySearch(size).run(adjacency);
    std::vector<int> position(size);
    for (int index = 0; index < size; index++) position[order[index]] = index;

    std::vector<int> parent(size, -1);
    std::vector<std::vector<int>> children(size);
    for (int vertex = 0; vertex < size; vertex++) {
        for (int to : adjacency[vertex]) {
            if (position[to] < position[vertex] &&
                (parent[vertex] == -1 || position[parent[vertex]] < position[to])) {
                parent[vertex] = to;
            }
        }
        if (parent[vertex] != -1) children[parent[vertex]].push_back(vertex);
    }

    std::vector<int> adjacent_stamp(size, -1);
    for (int center = 0; center < size; center++) {
        for (int to : adjacency[center]) adjacent_stamp[to] = center;
        for (int vertex : children[center]) {
            for (int to : adjacency[vertex]) {
                if (position[to] >= position[center] || adjacent_stamp[to] == center) continue;
                return ChordalGraphResult{
                    false,
                    {},
                    internal::chordless_cycle(adjacency, vertex, to, center),
                };
            }
        }
    }

    std::reverse(order.begin(), order.end());
    return ChordalGraphResult{true, std::move(order), {}};
}

template <class T>
bool is_chordal(const Graph<T>& graph) {
    return chordal_graph_recognition(graph).is_chordal;
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/chromatic_number.hpp"



#line 5 "graph/chromatic_number.hpp"
#include <bit>
#line 9 "graph/chromatic_number.hpp"

#line 11 "graph/chromatic_number.hpp"

namespace m1une {
namespace graph {

namespace detail {

struct ChromaticResidues {
    static constexpr std::array<std::uint32_t, 14> mod = {
        1000000007, 1000000009, 998244353, 985661441, 943718401, 935329793, 918552577,
        897581057,  880803841,  754974721, 645922817, 595591169, 469762049, 167772161,
    };

    std::array<std::uint32_t, 14> value;

    explicit ChromaticResidues(std::uint32_t x = 0) {
        value.fill(x);
    }

    void multiply(std::uint32_t x) {
        for (int i = 0; i < int(mod.size()); i++) {
            value[i] = std::uint32_t(std::uint64_t(value[i]) * x % mod[i]);
        }
    }
};

}  // namespace detail

template <class T>
int chromatic_number(const Graph<T>& g) {
    int n = g.size();
    assert(n <= 20);
    if (n == 0) return 0;

    std::vector<std::uint32_t> adjacent(n, 0);
    for (const auto& e : g.edges()) {
        if (e.from == e.to) continue;
        adjacent[e.from] |= std::uint32_t(1) << e.to;
        adjacent[e.to] |= std::uint32_t(1) << e.from;
    }

    std::uint32_t subset_count = std::uint32_t(1) << n;
    std::vector<std::uint32_t> independent_count(subset_count, 0);
    independent_count[0] = 1;
    for (std::uint32_t mask = 1; mask < subset_count; mask++) {
        int v = std::countr_zero(mask);
        std::uint32_t rest = mask ^ (std::uint32_t(1) << v);
        independent_count[mask] =
            independent_count[rest] + independent_count[rest & ~adjacent[v]];
    }

    std::vector<detail::ChromaticResidues> power(subset_count, detail::ChromaticResidues(1));
    for (int colors = 1; colors <= n; colors++) {
        std::array<std::uint32_t, 14> sum = {};
        for (std::uint32_t mask = 0; mask < subset_count; mask++) {
            power[mask].multiply(independent_count[mask]);
            bool positive = ((n - std::popcount(mask)) & 1) == 0;
            for (int i = 0; i < int(sum.size()); i++) {
                std::uint32_t x = power[mask].value[i];
                if (positive) {
                    sum[i] += x;
                    if (sum[i] >= detail::ChromaticResidues::mod[i]) {
                        sum[i] -= detail::ChromaticResidues::mod[i];
                    }
                } else {
                    sum[i] = (sum[i] >= x ? sum[i] - x
                                          : sum[i] + detail::ChromaticResidues::mod[i] - x);
                }
            }
        }
        for (std::uint32_t x : sum) {
            if (x != 0) return colors;
        }
    }
    return n;
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/complement_connected_components.hpp"



#line 6 "graph/complement_connected_components.hpp"

#line 1 "graph/connected_components.hpp"



#line 6 "graph/connected_components.hpp"

#line 1 "ds/dsu/dsu.hpp"



#line 5 "ds/dsu/dsu.hpp"
#include <numeric>
#line 8 "ds/dsu/dsu.hpp"

namespace m1une {
namespace ds {

struct Dsu {
   private:
    int _n;
    // parent_or_size[i] is the parent of i if it's >= 0.
    // If it's < 0, then i is a root and -parent_or_size[i] is the size of the group.
    std::vector<int> parent_or_size;

    // Returns {new leader, absorbed leader}. The absorbed leader is -1 when
    // both vertices already belong to the same component.
    std::pair<int, int> merge_leaders(int a, int b) {
        int x = leader(a), y = leader(b);
        if (x == y) return {x, -1};
        if (-parent_or_size[x] < -parent_or_size[y]) std::swap(x, y);
        parent_or_size[x] += parent_or_size[y];
        parent_or_size[y] = x;
        return {x, y};
    }

   public:
    Dsu() : _n(0) {}
    explicit Dsu(int n) : _n(n), parent_or_size(n, -1) {}

    // Merges the group containing 'a' with the group containing 'b'.
    // Returns the leader of the merged group.
    int merge(int a, int b) {
        return merge_leaders(a, b).first;
    }

    // Invokes callback(new_leader, absorbed_leader) after an actual merge.
    // Returns the leader of the merged group.
    template <class Callback>
    int merge(int a, int b, Callback&& callback) {
        std::pair<int, int> merged = merge_leaders(a, b);
        if (merged.second != -1) callback(merged.first, merged.second);
        return merged.first;
    }

    // Returns true if 'a' and 'b' belong to the same group.
    bool same(int a, int b) {
        return leader(a) == leader(b);
    }

    // Returns the leader (representative) of the group containing 'a'.
    int leader(int a) {
        if (parent_or_size[a] < 0) return a;
        // Path compression
        return parent_or_size[a] = leader(parent_or_size[a]);
    }

    // Returns the size of the group containing 'a'.
    int size(int a) {
        return -parent_or_size[leader(a)];
    }

    // Returns a list of all groups, where each group is a vector of its elements.
    std::vector<std::vector<int>> groups() {
        std::vector<int> leader_buf(_n), group_size(_n);
        for (int i = 0; i < _n; i++) {
            leader_buf[i] = leader(i);
            group_size[leader_buf[i]]++;
        }
        std::vector<std::vector<int>> result(_n);
        for (int i = 0; i < _n; i++) {
            result[i].reserve(group_size[i]);
        }
        for (int i = 0; i < _n; i++) {
            result[leader_buf[i]].push_back(i);
        }
        result.erase(std::remove_if(result.begin(), result.end(), [&](const std::vector<int>& v) { return v.empty(); }),
                     result.end());
        return result;
    }
};

}  // namespace ds
}  // namespace m1une


#line 9 "graph/connected_components.hpp"

namespace m1une {
namespace graph {

struct ConnectedComponents {
    int count;
    std::vector<int> comp;
    std::vector<std::vector<int>> groups;

    bool same(int u, int v) const {
        assert(0 <= u && u < int(comp.size()));
        assert(0 <= v && v < int(comp.size()));
        return comp[u] == comp[v];
    }
};

template <class T>
ConnectedComponents connected_components(const Graph<T>& g) {
    int n = g.size();
    m1une::ds::Dsu dsu(n);
    for (const auto& e : g.edges()) dsu.merge(e.from, e.to);

    ConnectedComponents result;
    result.comp.assign(n, 0);
    std::vector<int> leader_to_comp(n, -1);
    for (int v = 0; v < n; v++) {
        int leader = dsu.leader(v);
        if (leader_to_comp[leader] == -1) {
            leader_to_comp[leader] = int(result.groups.size());
            result.groups.push_back({});
        }
        int c = leader_to_comp[leader];
        result.comp[v] = c;
        result.groups[c].push_back(v);
    }
    result.count = int(result.groups.size());

    return result;
}

}  // namespace graph
}  // namespace m1une


#line 8 "graph/complement_connected_components.hpp"

namespace m1une {
namespace graph {

// Computes connected components after complementing the underlying simple
// undirected graph, without constructing the complement graph.
template <class T>
ConnectedComponents complement_connected_components(const Graph<T>& graph) {
    const int size = graph.size();
    std::vector<std::vector<int>> adjacency(size);
    for (const Edge<T>& edge : graph.edges()) {
        if (edge.from == edge.to) continue;
        adjacency[edge.from].push_back(edge.to);
        adjacency[edge.to].push_back(edge.from);
    }

    const int sentinel = size;
    std::vector<int> next(size + 1);
    std::vector<int> previous(size + 1);
    if (size == 0) {
        next[sentinel] = previous[sentinel] = sentinel;
    } else {
        next[sentinel] = 0;
        previous[sentinel] = size - 1;
        for (int vertex = 0; vertex < size; vertex++) {
            next[vertex] = (vertex + 1 == size ? sentinel : vertex + 1);
            previous[vertex] = (vertex == 0 ? sentinel : vertex - 1);
        }
    }

    auto erase = [&](int vertex) {
        next[previous[vertex]] = next[vertex];
        previous[next[vertex]] = previous[vertex];
    };

    ConnectedComponents result;
    result.comp.assign(size, -1);
    std::vector<int> neighbor_stamp(size, -1);
    std::queue<int> queue;

    while (next[sentinel] != sentinel) {
        const int root = next[sentinel];
        erase(root);
        const int component = int(result.groups.size());
        result.groups.emplace_back();
        result.groups.back().push_back(root);
        result.comp[root] = component;
        queue.push(root);

        while (!queue.empty()) {
            const int vertex = queue.front();
            queue.pop();
            for (int to : adjacency[vertex]) neighbor_stamp[to] = vertex;

            int candidate = next[sentinel];
            while (candidate != sentinel) {
                const int following = next[candidate];
                if (neighbor_stamp[candidate] != vertex) {
                    erase(candidate);
                    result.comp[candidate] = component;
                    result.groups.back().push_back(candidate);
                    queue.push(candidate);
                }
                candidate = following;
            }
        }
    }
    result.count = int(result.groups.size());
    return result;
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/count_four_cycles.hpp"



#line 6 "graph/count_four_cycles.hpp"
#include <tuple>
#line 9 "graph/count_four_cycles.hpp"

#line 11 "graph/count_four_cycles.hpp"

namespace m1une {
namespace graph {

namespace four_cycle_detail {

// Counts C4s containing one particular copy of each edge in a simple graph
// whose edge weights represent parallel-edge multiplicities.
inline std::vector<long long> count_simple_per_edge(
    int vertex_count,
    std::vector<int> first,
    std::vector<int> second,
    const std::vector<long long>& multiplicity
) {
    const int edge_count = int(first.size());
    assert(second.size() == first.size());
    assert(multiplicity.size() == first.size());

    std::vector<int> degree(vertex_count, 0);
    for (int edge = 0; edge < edge_count; edge++) {
        degree[first[edge]]++;
        degree[second[edge]]++;
    }

    int maximum_degree = 0;
    for (int value : degree) maximum_degree = std::max(maximum_degree, value);
    std::vector<int> degree_start(maximum_degree + 2, 0);
    for (int value : degree) degree_start[value + 1]++;
    for (int value = 0; value <= maximum_degree; value++) {
        degree_start[value + 1] += degree_start[value];
    }
    std::vector<int> cursor = degree_start;
    std::vector<int> order(vertex_count);
    for (int vertex = 0; vertex < vertex_count; vertex++) {
        order[cursor[degree[vertex]]++] = vertex;
    }
    std::vector<int> rank(vertex_count);
    for (int i = 0; i < vertex_count; i++) rank[order[i]] = i;
    for (int edge = 0; edge < edge_count; edge++) {
        first[edge] = rank[first[edge]];
        second[edge] = rank[second[edge]];
        if (first[edge] < second[edge]) {
            std::swap(first[edge], second[edge]);
        }
    }

    std::vector<int> start(vertex_count + 1, 0);
    for (int vertex = 0; vertex < vertex_count; vertex++) {
        start[vertex + 1] = start[vertex] + degree[order[vertex]];
    }
    std::vector<int> end = start;
    std::vector<int> edge_at(2 * edge_count);
    std::vector<int> to(2 * edge_count);
    for (int edge = 0; edge < edge_count; edge++) {
        int position = end[first[edge]]++;
        edge_at[position] = edge;
        to[position] = second[edge];
    }

    std::vector<int> downward_end = end;
    for (int vertex = 0; vertex < vertex_count; vertex++) {
        for (int i = start[vertex]; i < downward_end[vertex]; i++) {
            int edge = edge_at[i];
            int neighbor = to[i];
            int position = end[neighbor]++;
            edge_at[position] = edge;
            to[position] = vertex;
        }
    }

    std::vector<long long> path_count(vertex_count, 0);
    std::vector<long long> result(edge_count, 0);
    for (int vertex = vertex_count - 1; vertex >= 0; vertex--) {
        for (int i = start[vertex]; i < end[vertex]; i++) {
            int first_edge = edge_at[i];
            int middle = to[i];
            end[middle]--;
            for (int j = start[middle]; j < end[middle]; j++) {
                int second_edge = edge_at[j];
                int opposite = to[j];
                path_count[opposite] +=
                    multiplicity[first_edge] * multiplicity[second_edge];
            }
        }

        for (int i = start[vertex]; i < end[vertex]; i++) {
            int first_edge = edge_at[i];
            int middle = to[i];
            for (int j = start[middle]; j < end[middle]; j++) {
                int second_edge = edge_at[j];
                int opposite = to[j];
                long long other_paths =
                    path_count[opposite] -
                    multiplicity[first_edge] * multiplicity[second_edge];
                result[first_edge] +=
                    other_paths * multiplicity[second_edge];
                result[second_edge] +=
                    other_paths * multiplicity[first_edge];
            }
        }

        for (int i = start[vertex]; i < end[vertex]; i++) {
            int middle = to[i];
            for (int j = start[middle]; j < end[middle]; j++) {
                path_count[to[j]] = 0;
            }
        }
    }
    return result;
}

}  // namespace four_cycle_detail

// Returns, for every graph edge id, the number of C4 subgraphs containing it.
// Parallel active edges are distinct choices; inactive edges receive zero.
template <class T>
std::vector<long long> count_four_cycles_per_edge(const Graph<T>& graph) {
    struct ActiveEdge {
        int first;
        int second;
        int id;
    };

    std::vector<ActiveEdge> active_edges;
    active_edges.reserve(graph.edge_count());
    for (const Edge<T>& edge : graph.edges()) {
        assert(edge.from != edge.to);
        assert(0 <= edge.id && edge.id < graph.edge_count());
        if (edge.from == edge.to) continue;
        active_edges.push_back(ActiveEdge{
            std::min(edge.from, edge.to),
            std::max(edge.from, edge.to),
            edge.id
        });
    }
    std::sort(
        active_edges.begin(),
        active_edges.end(),
        [](const ActiveEdge& left, const ActiveEdge& right) {
            return std::tie(left.first, left.second) <
                   std::tie(right.first, right.second);
        }
    );

    std::vector<int> first;
    std::vector<int> second;
    std::vector<long long> multiplicity;
    std::vector<int> group_of_edge(graph.edge_count(), -1);
    first.reserve(active_edges.size());
    second.reserve(active_edges.size());
    multiplicity.reserve(active_edges.size());
    for (const ActiveEdge& edge : active_edges) {
        if (first.empty() || first.back() != edge.first ||
            second.back() != edge.second) {
            first.push_back(edge.first);
            second.push_back(edge.second);
            multiplicity.push_back(0);
        }
        multiplicity.back()++;
        group_of_edge[edge.id] = int(first.size()) - 1;
    }

    std::vector<long long> simple_result =
        four_cycle_detail::count_simple_per_edge(
            graph.size(),
            std::move(first),
            std::move(second),
            multiplicity
        );
    std::vector<long long> result(graph.edge_count(), 0);
    for (const ActiveEdge& edge : active_edges) {
        result[edge.id] = simple_result[group_of_edge[edge.id]];
    }
    return result;
}

template <class T>
long long count_four_cycles(const Graph<T>& graph) {
    std::vector<long long> per_edge = count_four_cycles_per_edge(graph);
    long long incidence_count = 0;
    for (long long count : per_edge) incidence_count += count;
    assert(incidence_count % 4 == 0);
    return incidence_count / 4;
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/cycle_detection.hpp"



#line 7 "graph/cycle_detection.hpp"

#line 9 "graph/cycle_detection.hpp"

namespace m1une {
namespace graph {

struct Cycle {
    std::vector<int> vertices;
    std::vector<int> edge_ids;

    bool empty() const {
        return vertices.empty();
    }
};

inline Cycle restore_cycle(int from, int to, int closing_edge, const std::vector<int>& parent,
                           const std::vector<int>& parent_edge) {
    Cycle result;
    result.vertices.push_back(to);

    std::vector<int> middle_vertices;
    std::vector<int> middle_edges;
    for (int v = from; v != to; v = parent[v]) {
        middle_vertices.push_back(v);
        middle_edges.push_back(parent_edge[v]);
    }
    std::reverse(middle_vertices.begin(), middle_vertices.end());
    std::reverse(middle_edges.begin(), middle_edges.end());

    result.vertices.insert(result.vertices.end(), middle_vertices.begin(), middle_vertices.end());
    result.vertices.push_back(to);
    result.edge_ids.insert(result.edge_ids.end(), middle_edges.begin(), middle_edges.end());
    result.edge_ids.push_back(closing_edge);
    return result;
}

template <class T>
Cycle find_directed_cycle(const Graph<T>& g) {
    int n = g.size();
    std::vector<int> color(n, 0), parent(n, -1), parent_edge(n, -1);
    struct Frame {
        int vertex;
        std::size_t next_edge;
    };

    std::vector<Frame> stack;
    stack.reserve(n);
    for (int start = 0; start < n; start++) {
        if (color[start] != 0) continue;
        color[start] = 1;
        stack.push_back(Frame{start, 0});
        while (!stack.empty()) {
            Frame& frame = stack.back();
            const int vertex = frame.vertex;
            const auto& adjacency = g[vertex];
            while (
                frame.next_edge < adjacency.size() &&
                !adjacency[frame.next_edge].alive
            ) {
                frame.next_edge++;
            }
            if (frame.next_edge == adjacency.size()) {
                color[vertex] = 2;
                stack.pop_back();
                continue;
            }

            const auto& edge = adjacency[frame.next_edge++];
            const int to = edge.to;
            const int edge_id = edge.id;
            if (color[to] == 0) {
                parent[to] = vertex;
                parent_edge[to] = edge_id;
                color[to] = 1;
                stack.push_back(Frame{to, 0});
            } else if (color[to] == 1) {
                return restore_cycle(vertex, to, edge_id, parent, parent_edge);
            }
        }
    }
    return Cycle();
}

template <class T>
Cycle find_undirected_cycle(const Graph<T>& g) {
    int n = g.size();
    std::vector<int> color(n, 0), parent(n, -1), parent_edge(n, -1);
    struct Frame {
        int vertex;
        std::size_t next_edge;
    };

    std::vector<Frame> stack;
    stack.reserve(n);
    for (int start = 0; start < n; start++) {
        if (color[start] != 0) continue;
        color[start] = 1;
        stack.push_back(Frame{start, 0});
        while (!stack.empty()) {
            Frame& frame = stack.back();
            const int vertex = frame.vertex;
            const auto& adjacency = g[vertex];
            while (
                frame.next_edge < adjacency.size() &&
                (
                    !adjacency[frame.next_edge].alive ||
                    adjacency[frame.next_edge].id == parent_edge[vertex]
                )
            ) {
                frame.next_edge++;
            }
            if (frame.next_edge == adjacency.size()) {
                color[vertex] = 2;
                stack.pop_back();
                continue;
            }

            const auto& edge = adjacency[frame.next_edge++];
            const int to = edge.to;
            const int edge_id = edge.id;
            if (color[to] == 0) {
                parent[to] = vertex;
                parent_edge[to] = edge_id;
                color[to] = 1;
                stack.push_back(Frame{to, 0});
            } else if (color[to] == 1) {
                return restore_cycle(vertex, to, edge_id, parent, parent_edge);
            }
        }
    }
    return Cycle();
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/dfs.hpp"



#line 6 "graph/dfs.hpp"
#include <concepts>
#include <functional>
#line 10 "graph/dfs.hpp"

#line 12 "graph/dfs.hpp"

namespace m1une {
namespace graph {

struct DfsResult {
    std::vector<int> depth;
    std::vector<int> parent;
    std::vector<int> parent_edge;
    std::vector<int> root;
    std::vector<int> tin;
    std::vector<int> tout;
    std::vector<int> preorder;
    std::vector<int> postorder;
    std::vector<int> roots;

    bool reachable(int vertex) const {
        assert(0 <= vertex && vertex < int(depth.size()));
        return depth[vertex] != -1;
    }

    int component_count() const {
        return int(roots.size());
    }

    std::vector<int> path(int target) const {
        assert(reachable(target));
        std::vector<int> result;
        for (int vertex = target; vertex != -1; vertex = parent[vertex]) {
            result.push_back(vertex);
        }
        std::reverse(result.begin(), result.end());
        return result;
    }

    bool is_ancestor(int ancestor, int vertex) const {
        assert(0 <= ancestor && ancestor < int(depth.size()));
        assert(0 <= vertex && vertex < int(depth.size()));
        if (!reachable(ancestor) || !reachable(vertex)) return false;
        return tin[ancestor] <= tin[vertex] && tout[vertex] <= tout[ancestor];
    }
};

namespace dfs_detail {

template <class Callback>
concept DfsCallback =
    std::invocable<Callback&, int, int> ||
    std::invocable<Callback&, int>;

template <DfsCallback Callback>
void invoke_callback(Callback& callback, int vertex, int parent) {
    if constexpr (std::invocable<Callback&, int, int>) {
        std::invoke(callback, vertex, parent);
    } else {
        std::invoke(callback, vertex);
    }
}

template <class T, class Callback>
DfsResult run_dfs(
    const Graph<T>& graph,
    const std::vector<int>& sources,
    bool complete_forest,
    Callback& callback
) {
    const int n = graph.size();
    DfsResult result;
    result.depth.assign(n, -1);
    result.parent.assign(n, -1);
    result.parent_edge.assign(n, -1);
    result.root.assign(n, -1);
    result.tin.assign(n, -1);
    result.tout.assign(n, -1);
    result.preorder.reserve(n);
    result.postorder.reserve(n);
    result.roots.reserve(n);

    struct Frame {
        int vertex;
        int next_edge;
    };
    std::vector<Frame> stack;
    stack.reserve(n);
    int timer = 0;

    auto traverse = [&](int source) {
        assert(0 <= source && source < n);
        if (result.reachable(source)) return;

        result.depth[source] = 0;
        result.root[source] = source;
        result.tin[source] = ++timer;
        result.preorder.push_back(source);
        result.roots.push_back(source);
        invoke_callback(callback, source, -1);
        stack.push_back(Frame{source, 0});

        while (!stack.empty()) {
            Frame& frame = stack.back();
            int vertex = frame.vertex;
            if (frame.next_edge == int(graph[vertex].size())) {
                result.tout[vertex] = ++timer;
                result.postorder.push_back(vertex);
                stack.pop_back();
                continue;
            }

            const Edge<T>& edge = graph[vertex][frame.next_edge++];
            if (!edge.alive || result.reachable(edge.to)) continue;
            result.depth[edge.to] = result.depth[vertex] + 1;
            result.parent[edge.to] = vertex;
            result.parent_edge[edge.to] = edge.id;
            result.root[edge.to] = result.root[vertex];
            result.tin[edge.to] = ++timer;
            result.preorder.push_back(edge.to);
            invoke_callback(callback, edge.to, vertex);
            stack.push_back(Frame{edge.to, 0});
        }
    };

    for (int source : sources) traverse(source);
    if (complete_forest) {
        for (int vertex = 0; vertex < n; vertex++) traverse(vertex);
    }
    return result;
}

}  // namespace dfs_detail

template <class T>
DfsResult dfs(const Graph<T>& graph, const std::vector<int>& sources) {
    auto callback = [](int) {};
    return dfs_detail::run_dfs(graph, sources, false, callback);
}

template <class T>
DfsResult dfs(const Graph<T>& graph, int source) {
    return dfs(graph, std::vector<int>{source});
}

template <class T>
DfsResult dfs(const Graph<T>& graph) {
    auto callback = [](int) {};
    return dfs_detail::run_dfs(
        graph,
        std::vector<int>(),
        true,
        callback
    );
}

template <class T, class Callback>
requires dfs_detail::DfsCallback<Callback>
DfsResult dfs(
    const Graph<T>& graph,
    const std::vector<int>& sources,
    Callback&& callback
) {
    return dfs_detail::run_dfs(graph, sources, false, callback);
}

template <class T, class Callback>
requires dfs_detail::DfsCallback<Callback>
DfsResult dfs(const Graph<T>& graph, int source, Callback&& callback) {
    return dfs(
        graph,
        std::vector<int>{source},
        std::forward<Callback>(callback)
    );
}

template <class T, class Callback>
requires dfs_detail::DfsCallback<Callback>
DfsResult dfs(const Graph<T>& graph, Callback&& callback) {
    return dfs_detail::run_dfs(
        graph,
        std::vector<int>(),
        true,
        callback
    );
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/enumerate_cliques.hpp"



#line 8 "graph/enumerate_cliques.hpp"

#line 10 "graph/enumerate_cliques.hpp"

namespace m1une {
namespace graph {

// Invokes callback once for every nonempty clique. The callback receives a
// const reference to a temporary vector that is reused after it returns.
template <class T, class Callback>
void enumerate_cliques(const Graph<T>& graph, Callback&& callback) {
    const int n = graph.size();
    std::vector<std::vector<int>> adjacency(n);
    for (const Edge<T>& edge : graph.edges()) {
        assert(edge.from != edge.to);
        if (edge.from == edge.to) continue;
        adjacency[edge.from].push_back(edge.to);
        adjacency[edge.to].push_back(edge.from);
    }

    for (std::vector<int>& neighbors : adjacency) {
        std::sort(neighbors.begin(), neighbors.end());
#ifndef NDEBUG
        for (int i = 1; i < int(neighbors.size()); i++) {
            assert(neighbors[i - 1] != neighbors[i]);
        }
#endif
        neighbors.erase(
            std::unique(neighbors.begin(), neighbors.end()),
            neighbors.end()
        );
    }

    int maximum_degree = 0;
    std::vector<int> degree(n);
    for (int vertex = 0; vertex < n; vertex++) {
        degree[vertex] = int(adjacency[vertex].size());
        maximum_degree = std::max(maximum_degree, degree[vertex]);
    }

    // Compute a degeneracy ordering in linear time. A clique is assigned to
    // its first vertex in this ordering, and all its other vertices are among
    // that vertex's forward neighbors.
    std::vector<std::vector<int>> bucket(maximum_degree + 1);
    for (int vertex = 0; vertex < n; vertex++) {
        bucket[degree[vertex]].push_back(vertex);
    }
    std::vector<char> active(n, true);
    std::vector<std::vector<int>> forward(n);
    int minimum_degree = 0;
    int degeneracy = 0;
    for (int removed = 0; removed < n; removed++) {
        while (true) {
            while (bucket[minimum_degree].empty()) minimum_degree++;
            int vertex = bucket[minimum_degree].back();
            if (active[vertex] && degree[vertex] == minimum_degree) break;
            bucket[minimum_degree].pop_back();
        }

        int vertex = bucket[minimum_degree].back();
        bucket[minimum_degree].pop_back();
        active[vertex] = false;
        degeneracy = std::max(degeneracy, minimum_degree);
        forward[vertex].reserve(minimum_degree);
        for (int to : adjacency[vertex]) {
            if (!active[to]) continue;
            forward[vertex].push_back(to);
            degree[to]--;
            bucket[degree[to]].push_back(to);
            minimum_degree = std::min(minimum_degree, degree[to]);
        }
    }

    std::vector<int> clique;
    clique.reserve(degeneracy + 1);
    std::vector<std::vector<int>> candidates(degeneracy + 1);
    for (int vertex = 0; vertex < n; vertex++) {
        const std::vector<int>& neighbors = forward[vertex];
        const int neighbor_count = int(neighbors.size());

        clique.clear();
        clique.push_back(vertex);
        callback(std::as_const(clique));
        if (neighbor_count == 0) continue;

        std::vector<char> connected(
            std::size_t(neighbor_count) * neighbor_count,
            false
        );
        for (int first = 0; first < neighbor_count; first++) {
            for (int second = first + 1; second < neighbor_count; second++) {
                bool adjacent = std::binary_search(
                    adjacency[neighbors[first]].begin(),
                    adjacency[neighbors[first]].end(),
                    neighbors[second]
                );
                connected[std::size_t(first) * neighbor_count + second] =
                    adjacent;
                connected[std::size_t(second) * neighbor_count + first] =
                    adjacent;
            }
        }

        candidates[0].resize(neighbor_count);
        for (int i = 0; i < neighbor_count; i++) candidates[0][i] = i;
        auto enumerate = [&](auto&& self, int depth) -> void {
            const std::vector<int>& current = candidates[depth];
            for (int position = 0; position < int(current.size()); position++) {
                int chosen = current[position];
                clique.push_back(neighbors[chosen]);
                callback(std::as_const(clique));

                std::vector<int>& next = candidates[depth + 1];
                next.clear();
                for (int next_position = position + 1;
                     next_position < int(current.size());
                     next_position++) {
                    int candidate = current[next_position];
                    if (connected[
                            std::size_t(chosen) * neighbor_count + candidate
                        ]) {
                        next.push_back(candidate);
                    }
                }
                if (!next.empty()) self(self, depth + 1);
                clique.pop_back();
            }
        };
        enumerate(enumerate, 0);
    }
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/enumerate_triangles.hpp"



#line 7 "graph/enumerate_triangles.hpp"

#line 9 "graph/enumerate_triangles.hpp"

namespace m1une {
namespace graph {

template <class T, class Callback>
void enumerate_triangles(const Graph<T>& graph, Callback&& callback) {
    const int n = graph.size();
    const std::vector<Edge<T>> edges = graph.edges();

    std::vector<int> degree(n, 0);
    for (const Edge<T>& edge : edges) {
        assert(edge.from != edge.to);
        degree[edge.from]++;
        degree[edge.to]++;
    }

    std::vector<std::vector<int>> oriented(n);
    for (const Edge<T>& edge : edges) {
        int from = edge.from;
        int to = edge.to;
        if (degree[from] > degree[to] ||
            (degree[from] == degree[to] && from > to)) {
            std::swap(from, to);
        }
        oriented[from].push_back(to);
    }

    std::vector<int> marked(n, -1);
    for (int vertex = 0; vertex < n; vertex++) {
        for (int to : oriented[vertex]) marked[to] = vertex;
        for (int middle : oriented[vertex]) {
            for (int to : oriented[middle]) {
                if (marked[to] != vertex) continue;
                int first = vertex;
                int second = middle;
                int third = to;
                if (first > second) std::swap(first, second);
                if (second > third) std::swap(second, third);
                if (first > second) std::swap(first, second);
                callback(first, second, third);
            }
        }
    }
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/eulerian_trail.hpp"



#line 9 "graph/eulerian_trail.hpp"

#line 11 "graph/eulerian_trail.hpp"

namespace m1une {
namespace graph {

struct EulerianTrail {
    std::vector<int> vertices;
    std::vector<int> edge_ids;

    int edge_count() const {
        return int(edge_ids.size());
    }

    bool is_circuit() const {
        return vertices.empty() || vertices.front() == vertices.back();
    }
};

namespace internal {

template <class T>
std::optional<EulerianTrail> hierholzer(
    const Graph<T>& graph,
    int start,
    int active_edge_count
) {
    EulerianTrail result;
    if (active_edge_count == 0) {
        if (start != -1) result.vertices.push_back(start);
        return result;
    }

    assert(0 <= start && start < graph.size());
    std::vector<char> used(graph.edge_count(), false);
    std::vector<int> cursor(graph.size(), 0);
    std::vector<int> vertex_stack(1, start);
    std::vector<int> incoming_edge_stack(1, -1);
    std::vector<int> reversed_vertices;
    std::vector<int> reversed_edges;
    reversed_vertices.reserve(active_edge_count + 1);
    reversed_edges.reserve(active_edge_count);

    while (!vertex_stack.empty()) {
        const int vertex = vertex_stack.back();
        while (cursor[vertex] < int(graph[vertex].size())) {
            const Edge<T>& edge = graph[vertex][cursor[vertex]];
            if (edge.alive && !used[edge.id]) break;
            cursor[vertex]++;
        }

        if (cursor[vertex] < int(graph[vertex].size())) {
            const Edge<T>& edge = graph[vertex][cursor[vertex]++];
            used[edge.id] = true;
            vertex_stack.push_back(edge.to);
            incoming_edge_stack.push_back(edge.id);
            continue;
        }

        reversed_vertices.push_back(vertex);
        const int incoming_edge = incoming_edge_stack.back();
        if (incoming_edge != -1) reversed_edges.push_back(incoming_edge);
        vertex_stack.pop_back();
        incoming_edge_stack.pop_back();
    }

    if (int(reversed_edges.size()) != active_edge_count) return std::nullopt;
    std::reverse(reversed_vertices.begin(), reversed_vertices.end());
    std::reverse(reversed_edges.begin(), reversed_edges.end());
    result.vertices = std::move(reversed_vertices);
    result.edge_ids = std::move(reversed_edges);
    return result;
}

template <class T>
std::vector<int> edge_incidence_count(const Graph<T>& graph) {
    std::vector<int> count(graph.edge_count(), 0);
    for (int vertex = 0; vertex < graph.size(); vertex++) {
        for (const Edge<T>& edge : graph[vertex]) {
            if (!edge.alive) continue;
            assert(0 <= edge.id && edge.id < graph.edge_count());
            count[edge.id]++;
        }
    }
    return count;
}

}  // namespace internal

template <class T>
std::optional<EulerianTrail> directed_eulerian_trail(
    const Graph<T>& graph,
    int start = -1
) {
    assert(start == -1 || (0 <= start && start < graph.size()));
    const int n = graph.size();
    std::vector<int> incidence = internal::edge_incidence_count(graph);
    std::vector<int> in_degree(n, 0);
    std::vector<int> out_degree(n, 0);
    int active_edge_count = 0;
    for (int vertex = 0; vertex < n; vertex++) {
        for (const Edge<T>& edge : graph[vertex]) {
            if (!edge.alive) continue;
            out_degree[vertex]++;
            in_degree[edge.to]++;
        }
    }
    for (int count : incidence) {
        if (count == 0) continue;
        assert(count == 1);
        active_edge_count++;
    }

    int required_start = -1;
    int required_end = -1;
    for (int vertex = 0; vertex < n; vertex++) {
        const int difference = out_degree[vertex] - in_degree[vertex];
        if (difference == 1) {
            if (required_start != -1) return std::nullopt;
            required_start = vertex;
        } else if (difference == -1) {
            if (required_end != -1) return std::nullopt;
            required_end = vertex;
        } else if (difference != 0) {
            return std::nullopt;
        }
    }
    if ((required_start == -1) != (required_end == -1)) return std::nullopt;

    int chosen_start = start;
    if (active_edge_count == 0) {
        if (chosen_start == -1 && n > 0) chosen_start = 0;
        return internal::hierholzer(graph, chosen_start, 0);
    }
    if (required_start != -1) {
        if (chosen_start != -1 && chosen_start != required_start) return std::nullopt;
        chosen_start = required_start;
    } else if (chosen_start == -1) {
        for (int vertex = 0; vertex < n; vertex++) {
            if (out_degree[vertex] > 0) {
                chosen_start = vertex;
                break;
            }
        }
    } else if (out_degree[chosen_start] == 0) {
        return std::nullopt;
    }
    return internal::hierholzer(graph, chosen_start, active_edge_count);
}

template <class T>
std::optional<EulerianTrail> undirected_eulerian_trail(
    const Graph<T>& graph,
    int start = -1
) {
    assert(start == -1 || (0 <= start && start < graph.size()));
    const int n = graph.size();
    std::vector<int> incidence = internal::edge_incidence_count(graph);
    std::vector<int> degree(n, 0);
    int active_edge_count = 0;
    for (int vertex = 0; vertex < n; vertex++) {
        for (const Edge<T>& edge : graph[vertex]) {
            if (edge.alive) degree[vertex]++;
        }
    }
    for (int count : incidence) {
        if (count == 0) continue;
        assert(count == 2);
        active_edge_count++;
    }

    std::vector<int> odd;
    for (int vertex = 0; vertex < n; vertex++) {
        if (degree[vertex] & 1) odd.push_back(vertex);
    }
    if (!odd.empty() && odd.size() != 2) return std::nullopt;

    int chosen_start = start;
    if (active_edge_count == 0) {
        if (chosen_start == -1 && n > 0) chosen_start = 0;
        return internal::hierholzer(graph, chosen_start, 0);
    }
    if (odd.size() == 2) {
        if (chosen_start != -1 && chosen_start != odd[0] && chosen_start != odd[1]) {
            return std::nullopt;
        }
        if (chosen_start == -1) chosen_start = odd[0];
    } else if (chosen_start == -1) {
        for (int vertex = 0; vertex < n; vertex++) {
            if (degree[vertex] > 0) {
                chosen_start = vertex;
                break;
            }
        }
    } else if (degree[chosen_start] == 0) {
        return std::nullopt;
    }
    return internal::hierholzer(graph, chosen_start, active_edge_count);
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/general_matching.hpp"



#line 9 "graph/general_matching.hpp"

#line 11 "graph/general_matching.hpp"

namespace m1une {
namespace graph {

struct GeneralMatching {
    struct Edge {
        int from;
        int to;
        int id;
        bool alive;

        int other(int v) const {
            assert(v == from || v == to);
            return from ^ to ^ v;
        }
    };

    struct Pair {
        int from;
        int to;
        int edge_id;
    };

   private:
    int _n;
    std::vector<Edge> _edges;
    std::vector<std::vector<int>> _adj;
    std::vector<int> _mate;
    std::vector<int> _mate_edge;
    bool _calculated;

    void invalidate() {
        _calculated = false;
    }

    void ensure_matching() {
        if (!_calculated) max_matching();
    }

    bool is_matched_edge(int id) const {
        const auto& e = _edges[id];
        return _mate[e.from] == e.to && _mate_edge[e.from] == id;
    }

    enum MatchingLabel : char {
        even_label,
        odd_label,
        unlabeled
    };

    struct MutablePartition {
        std::vector<int> parent;
        std::vector<int> rank;
        std::vector<int> representative;

        MutablePartition() = default;

        explicit MutablePartition(int n) {
            reset(n);
        }

        void reset(int n) {
            parent.resize(n);
            rank.assign(n, 0);
            representative.resize(n);
            for (int i = 0; i < n; i++) {
                parent[i] = i;
                representative[i] = i;
            }
        }

        int root(int v) {
            if (parent[v] == v) return v;
            return parent[v] = root(parent[v]);
        }

        int operator()(int v) {
            return representative[root(v)];
        }

        void unite(int a, int b) {
            int ra = root(a);
            int rb = root(b);
            if (ra == rb) return;
            if (rank[ra] < rank[rb]) std::swap(ra, rb);
            parent[rb] = ra;
            if (rank[ra] == rank[rb]) rank[ra]++;
        }

        void make_rep(int v) {
            representative[root(v)] = v;
        }
    };

    struct EdgeBucketQueue {
        std::vector<std::vector<int>> bucket;
        std::vector<int> head;

        void reset(int n) {
            bucket.assign(n + 3, {});
            head.assign(n + 3, 0);
        }

        void insert(int edge_id, int key) {
            if (key < 0 || int(bucket.size()) <= key) return;
            bucket[key].push_back(edge_id);
        }

        int pop(int key) {
            if (key < 0 || int(bucket.size()) <= key) return -1;
            if (head[key] == int(bucket[key].size())) return -1;
            return bucket[key][head[key]++];
        }
    };

    struct NewMatchingPair {
        int from;
        int to;
        int edge_id;
    };

    // General-graph shortest augmenting path phase solver.
    struct MicaliVaziraniSolver {
        GeneralMatching& graph;
        int n;
        int matching_size;
        int delta;
        int visit_token;
        int even_time_token;
        MutablePartition base;
        MutablePartition delayed_base;
        EdgeBucketQueue queue;
        std::vector<MatchingLabel> label;
        std::vector<MatchingLabel> h_label;
        std::vector<int> parent;
        std::vector<int> parent_edge;
        std::vector<int> source_bridge;
        std::vector<int> target_bridge;
        std::vector<int> bridge_edge;
        std::vector<int> lcp;
        std::vector<int> path_mark_1;
        std::vector<int> path_mark_2;
        std::vector<int> restore_vertex;
        std::vector<int> restore_value;
        std::vector<int> rep;
        std::vector<int> h_mate;
        std::vector<char> is_h_edge;
        std::vector<std::vector<int>> contracted_into;
        std::vector<int> h_parent_edge;
        std::vector<int> h_even_time;
        std::vector<int> h_bridge_edge;
        std::vector<int> h_bridge_dir;

        explicit MicaliVaziraniSolver(GeneralMatching& graph_)
            : graph(graph_),
              n(graph_._n),
              matching_size(0),
              delta(0),
              visit_token(0),
              even_time_token(0),
              base(n),
              delayed_base(n),
              label(n, unlabeled),
              h_label(n, unlabeled),
              parent(n, -1),
              parent_edge(n, -1),
              source_bridge(n, -1),
              target_bridge(n, -1),
              bridge_edge(n, -1),
              lcp(n, 0),
              path_mark_1(n, 0),
              path_mark_2(n, 0),
              rep(n, -1),
              h_mate(n, -1),
              is_h_edge(graph_._edges.size(), false),
              contracted_into(n),
              h_parent_edge(n, -1),
              h_even_time(n, 0),
              h_bridge_edge(n, -1),
              h_bridge_dir(n, 0) {}

        bool active(int edge_id) const {
            return graph._edges[edge_id].alive;
        }

        int other(int edge_id, int v) const {
            return graph._edges[edge_id].other(v);
        }

        int edge_weight(int edge_id) const {
            return graph.is_matched_edge(edge_id) ? 2 : 0;
        }

        void set_match(int edge_id) {
            const auto& e = graph._edges[edge_id];
            graph._mate[e.from] = e.to;
            graph._mate[e.to] = e.from;
            graph._mate_edge[e.from] = edge_id;
            graph._mate_edge[e.to] = edge_id;
        }

        void initialize_greedy_matching() {
            graph._mate.assign(n, -1);
            graph._mate_edge.assign(n, -1);
            matching_size = 0;
            for (const auto& e : graph._edges) {
                if (!e.alive) continue;
                if (graph._mate[e.from] != -1 || graph._mate[e.to] != -1) continue;
                set_match(e.id);
                matching_size++;
            }
        }

        void scan_edge(int edge_id, int from) {
            if (!active(edge_id)) return;
            int to = other(edge_id, from);
            if (to == from || graph._mate[to] == from || label[base(to)] == odd_label) return;
            if (label[to] == unlabeled) {
                queue.insert(edge_id, lcp[from] + 2);
            } else {
                queue.insert(edge_id, (lcp[from] + lcp[to]) / 2 + 1);
            }
        }

        void shrink_path(int blossom_base, int x, int y, int edge_id,
                         std::vector<std::pair<int, int>>& delayed_unions) {
            int v = base(x);
            while (v != blossom_base) {
                base.unite(v, blossom_base);
                delayed_unions.push_back({v, blossom_base});

                v = graph._mate[v];
                assert(v != -1);
                base.unite(v, blossom_base);
                delayed_unions.push_back({v, blossom_base});
                base.make_rep(blossom_base);

                source_bridge[v] = x;
                target_bridge[v] = y;
                bridge_edge[v] = edge_id;
                restore_vertex.push_back(v);
                restore_value.push_back(lcp[v]);
                lcp[v] = lcp[x] + lcp[y] - lcp[graph._mate[v]] + 2;

                for (int id : graph._adj[v]) scan_edge(id, v);
                assert(parent[v] != -1);
                v = base(parent[v]);
            }
            delayed_unions.push_back({blossom_base, blossom_base});
        }

        void build_phase_graph() {
            std::fill(h_mate.begin(), h_mate.end(), -1);
            std::fill(is_h_edge.begin(), is_h_edge.end(), false);
            for (auto& vertices : contracted_into) vertices.clear();

            for (int v = 0; v < n; v++) contracted_into[delayed_base(v)].push_back(v);

            for (const auto& e : graph._edges) {
                if (!e.alive) continue;
                int u = e.from;
                int v = e.to;
                int uh = delayed_base(u);
                int vh = delayed_base(v);
                if (uh == vh) continue;
                if (label[uh] == odd_label && label[vh] == odd_label) continue;

                int w = edge_weight(e.id);
                bool even_odd =
                    (label[uh] == even_label && label[vh] == odd_label && lcp[v] == lcp[u] + 1 - w) ||
                    (label[vh] == even_label && label[uh] == odd_label && lcp[u] == lcp[v] + 1 - w);
                bool unlabeled_unlabeled = label[uh] == unlabeled && label[vh] == unlabeled && w == 2;
                bool even_unlabeled =
                    (label[uh] == even_label && label[vh] == unlabeled && lcp[u] == delta - 2) ||
                    (label[vh] == even_label && label[uh] == unlabeled && lcp[v] == delta - 2);
                bool even_even = label[uh] == even_label && label[vh] == even_label;
                bool tight_even_even = even_even && lcp[u] + lcp[v] == 2 * delta + w - 2;

                if (even_odd || unlabeled_unlabeled || even_unlabeled || tight_even_even) {
                    is_h_edge[e.id] = true;
                    if (w == 2) {
                        h_mate[uh] = vh;
                        h_mate[vh] = uh;
                    }
                }
            }
        }

        bool phase_one() {
            delta = 0;
            base.reset(n);
            delayed_base.reset(n);
            queue.reset(n);
            std::fill(label.begin(), label.end(), unlabeled);
            std::fill(parent.begin(), parent.end(), -1);
            std::fill(parent_edge.begin(), parent_edge.end(), -1);
            std::fill(source_bridge.begin(), source_bridge.end(), -1);
            std::fill(target_bridge.begin(), target_bridge.end(), -1);
            std::fill(bridge_edge.begin(), bridge_edge.end(), -1);
            std::fill(lcp.begin(), lcp.end(), 0);

            for (int v = 0; v < n; v++) {
                if (graph._mate[v] == -1) label[v] = even_label;
            }
            for (int v = 0; v < n; v++) {
                if (label[v] != even_label) continue;
                for (int id : graph._adj[v]) scan_edge(id, v);
            }

            std::vector<std::pair<int, int>> delayed_unions;
            while (delta <= n + 1) {
                restore_vertex.clear();
                restore_value.clear();

                while (true) {
                    int edge_id = queue.pop(delta);
                    if (edge_id == -1) break;
                    if (!active(edge_id)) continue;

                    int x = graph._edges[edge_id].from;
                    int y = graph._edges[edge_id].to;
                    if (label[base(x)] != even_label) std::swap(x, y);
                    if (label[base(x)] != even_label) continue;
                    if (graph._mate[x] == y || base(x) == base(y) || label[base(y)] == odd_label) continue;

                    if (label[base(y)] == unlabeled) {
                        int z = graph._mate[y];
                        assert(z != -1);
                        lcp[y] = lcp[x] + 1;
                        lcp[z] = lcp[x] + 2;
                        parent[y] = x;
                        parent_edge[y] = edge_id;
                        parent[z] = y;
                        parent_edge[z] = graph._mate_edge[z];
                        label[y] = odd_label;
                        label[z] = even_label;
                        for (int id : graph._adj[z]) scan_edge(id, z);
                        continue;
                    }

                    if (label[base(y)] != even_label || lcp[x] + lcp[y] != 2 * delta - 2) continue;

                    ++visit_token;
                    int hx = base(x);
                    int hy = base(y);
                    path_mark_1[hx] = visit_token;
                    path_mark_2[hy] = visit_token;
                    while (path_mark_1[hy] != visit_token && path_mark_2[hx] != visit_token &&
                           (graph._mate[hx] != -1 || graph._mate[hy] != -1)) {
                        if (graph._mate[hx] != -1) {
                            assert(parent[graph._mate[hx]] != -1);
                            hx = base(parent[graph._mate[hx]]);
                            path_mark_1[hx] = visit_token;
                        }
                        if (graph._mate[hy] != -1) {
                            assert(parent[graph._mate[hy]] != -1);
                            hy = base(parent[graph._mate[hy]]);
                            path_mark_2[hy] = visit_token;
                        }
                    }

                    if (path_mark_1[hy] == visit_token || path_mark_2[hx] == visit_token) {
                        int blossom_base = path_mark_1[hy] == visit_token ? hy : hx;
                        shrink_path(blossom_base, x, y, edge_id, delayed_unions);
                        shrink_path(blossom_base, y, x, edge_id, delayed_unions);
                    } else {
                        for (int i = int(restore_vertex.size()) - 1; i >= 0; i--) {
                            lcp[restore_vertex[i]] = restore_value[i];
                        }
                        build_phase_graph();
                        return true;
                    }
                }

                for (auto [a, b] : delayed_unions) {
                    if (a == b) {
                        delayed_base.make_rep(a);
                    } else {
                        delayed_base.unite(a, b);
                    }
                }
                delayed_unions.clear();
                delta++;
            }
            return false;
        }

        int next_h_vertex_through_edge(int edge_id, int current_h) const {
            const auto& e = graph._edges[edge_id];
            return rep[rep[e.from] == current_h ? e.to : e.from];
        }

        int find_path_in_h(int h_vertex) {
            for (int v : contracted_into[h_vertex]) {
                for (int edge_id : graph._adj[v]) {
                    if (!is_h_edge[edge_id]) continue;
                    int uh = rep[other(edge_id, v)];
                    if (h_mate[h_vertex] == uh) continue;

                    if (h_label[uh] == unlabeled) {
                        int mate_uh = h_mate[uh];
                        h_label[uh] = odd_label;
                        h_parent_edge[uh] = edge_id;
                        if (mate_uh == -1) return uh;

                        h_label[mate_uh] = even_label;
                        h_even_time[mate_uh] = even_time_token++;
                        int found = find_path_in_h(mate_uh);
                        if (found != -1) return found;
                    } else {
                        int bh = delayed_base(h_vertex);
                        int zh = delayed_base(uh);
                        if (h_even_time[bh] >= h_even_time[zh]) continue;

                        std::vector<int> blossom_path;
                        std::vector<int> blossom_vertices;
                        while (zh != bh) {
                            blossom_vertices.push_back(zh);
                            zh = h_mate[zh];
                            assert(zh != -1);
                            blossom_vertices.push_back(zh);
                            blossom_path.push_back(zh);
                            assert(h_parent_edge[zh] != -1);
                            zh = delayed_base(next_h_vertex_through_edge(h_parent_edge[zh], zh));
                        }

                        for (int x : blossom_vertices) delayed_base.unite(x, bh);
                        delayed_base.make_rep(bh);

                        std::reverse(blossom_path.begin(), blossom_path.end());
                        for (int x : blossom_path) {
                            h_bridge_edge[x] = edge_id;
                            h_bridge_dir[x] = graph._edges[edge_id].to == v ? 1 : -1;
                        }
                        for (int x : blossom_path) {
                            int found = find_path_in_h(x);
                            if (found != -1) return found;
                        }
                    }
                }
            }
            return -1;
        }

        void collect_path_in_h(std::vector<int>& path, int from_h, int to_h) {
            if (from_h == to_h) return;
            if (h_label[from_h] == even_label) {
                int mate_from = h_mate[from_h];
                assert(mate_from != -1);
                int edge_id = h_parent_edge[mate_from];
                assert(edge_id != -1);
                path.push_back(edge_id);
                collect_path_in_h(path, next_h_vertex_through_edge(edge_id, mate_from), to_h);
            } else {
                int edge_id = h_bridge_edge[from_h];
                assert(edge_id != -1);
                const auto& e = graph._edges[edge_id];
                int first = rep[h_bridge_dir[from_h] == 1 ? e.from : e.to];
                int second = rep[h_bridge_dir[from_h] == 1 ? e.to : e.from];
                collect_path_in_h(path, first, rep[h_mate[from_h]]);
                path.push_back(edge_id);
                collect_path_in_h(path, second, to_h);
            }
        }

        void add_new_pair(std::vector<NewMatchingPair>& pairs, int from, int to, int edge_id) const {
            const auto& e = graph._edges[edge_id];
            assert(e.alive);
            assert((e.from == from && e.to == to) || (e.from == to && e.to == from));
            pairs.push_back(NewMatchingPair{from, to, edge_id});
        }

        void collect_path_in_graph(std::vector<NewMatchingPair>& pairs, int from, int to) {
            if (from == to) return;
            if (label[from] == even_label) {
                int mate_from = graph._mate[from];
                assert(mate_from != -1);
                int parent_of_mate = parent[mate_from];
                int edge_id = parent_edge[mate_from];
                assert(parent_of_mate != -1 && edge_id != -1);
                add_new_pair(pairs, mate_from, parent_of_mate, edge_id);
                collect_path_in_graph(pairs, parent_of_mate, to);
            } else {
                assert(source_bridge[from] != -1 && target_bridge[from] != -1 && bridge_edge[from] != -1);
                collect_path_in_graph(pairs, source_bridge[from], graph._mate[from]);
                add_new_pair(pairs, source_bridge[from], target_bridge[from], bridge_edge[from]);
                collect_path_in_graph(pairs, target_bridge[from], to);
            }
        }

        void augment_path(const std::vector<int>& h_path) {
            std::vector<NewMatchingPair> pairs;
            for (int edge_id : h_path) {
                const auto& e = graph._edges[edge_id];
                add_new_pair(pairs, e.from, e.to, edge_id);
                collect_path_in_graph(pairs, e.from, rep[e.from]);
                collect_path_in_graph(pairs, e.to, rep[e.to]);
            }

            for (const auto& p : pairs) {
                if (graph._mate[p.from] != -1) {
                    int old = graph._mate[p.from];
                    graph._mate[old] = -1;
                    graph._mate_edge[old] = -1;
                }
                if (graph._mate[p.to] != -1) {
                    int old = graph._mate[p.to];
                    graph._mate[old] = -1;
                    graph._mate_edge[old] = -1;
                }
                graph._mate[p.from] = graph._mate[p.to] = -1;
                graph._mate_edge[p.from] = graph._mate_edge[p.to] = -1;
            }
            for (const auto& p : pairs) {
                assert(graph._mate[p.from] == -1 && graph._mate[p.to] == -1);
                graph._mate[p.from] = p.to;
                graph._mate[p.to] = p.from;
                graph._mate_edge[p.from] = p.edge_id;
                graph._mate_edge[p.to] = p.edge_id;
            }
            matching_size++;
        }

        void phase_two() {
            std::fill(h_label.begin(), h_label.end(), unlabeled);
            std::fill(h_parent_edge.begin(), h_parent_edge.end(), -1);
            std::fill(h_bridge_edge.begin(), h_bridge_edge.end(), -1);
            std::fill(h_bridge_dir.begin(), h_bridge_dir.end(), 0);
            for (int v = 0; v < n; v++) rep[v] = delayed_base(v);

            std::vector<std::vector<int>> paths;
            for (int h_vertex = 0; h_vertex < n; h_vertex++) {
                if (rep[h_vertex] != h_vertex) continue;
                if (h_label[h_vertex] != unlabeled || h_mate[h_vertex] != -1) continue;

                h_label[h_vertex] = even_label;
                h_even_time[h_vertex] = even_time_token++;
                int free_h = find_path_in_h(h_vertex);
                if (free_h == -1) continue;

                std::vector<int> path;
                int edge_id = h_parent_edge[free_h];
                assert(edge_id != -1);
                path.push_back(edge_id);
                collect_path_in_h(path, next_h_vertex_through_edge(edge_id, free_h), h_vertex);
                paths.push_back(path);
            }

            assert(!paths.empty());
            for (const auto& path : paths) augment_path(path);
            for (auto& vertices : contracted_into) vertices.clear();
        }

        int solve() {
            initialize_greedy_matching();
            while (phase_one()) phase_two();
            return matching_size;
        }
    };

   public:
    GeneralMatching() : GeneralMatching(0) {}

    explicit GeneralMatching(int n) : _n(n), _adj(n), _mate(n, -1), _mate_edge(n, -1), _calculated(false) {
        assert(0 <= n);
    }

    int size() const {
        return _n;
    }

    int edge_count() const {
        return int(_edges.size());
    }

    int add_edge(int from, int to) {
        assert(0 <= from && from < _n);
        assert(0 <= to && to < _n);
        assert(from != to);
        int id = int(_edges.size());
        _edges.push_back(Edge{from, to, id, true});
        _adj[from].push_back(id);
        _adj[to].push_back(id);
        invalidate();
        return id;
    }

    Edge get_edge(int i) const {
        assert(0 <= i && i < int(_edges.size()));
        return _edges[i];
    }

    std::vector<Edge> edges(bool include_inactive = false) const {
        std::vector<Edge> result;
        result.reserve(_edges.size());
        for (const auto& e : _edges) {
            if (include_inactive || e.alive) result.push_back(e);
        }
        return result;
    }

    void set_edge_alive(int id, bool alive) {
        assert(0 <= id && id < int(_edges.size()));
        _edges[id].alive = alive;
        invalidate();
    }

    void erase_edge(int id) {
        set_edge_alive(id, false);
    }

    void revive_edge(int id) {
        set_edge_alive(id, true);
    }

    bool is_edge_alive(int id) const {
        assert(0 <= id && id < int(_edges.size()));
        return _edges[id].alive;
    }

    int max_matching() {
        MicaliVaziraniSolver solver(*this);
        int result = solver.solve();

        _calculated = true;
        return result;
    }

    int matching_size() {
        ensure_matching();
        int result = 0;
        for (int v = 0; v < _n; v++) {
            if (v < _mate[v]) result++;
        }
        return result;
    }

    std::vector<int> mate() {
        ensure_matching();
        return _mate;
    }

    std::vector<int> mate_edge() {
        ensure_matching();
        return _mate_edge;
    }

    std::vector<Pair> matching() {
        ensure_matching();
        std::vector<Pair> result;
        for (int v = 0; v < _n; v++) {
            if (v < _mate[v]) result.push_back(Pair{v, _mate[v], _mate_edge[v]});
        }
        return result;
    }

    std::optional<std::vector<int>> minimum_edge_cover() {
        ensure_matching();

        std::vector<int> result;
        std::vector<char> covered(_n, false), used_edge(_edges.size(), false);

        auto use_edge = [&](int id) {
            if (used_edge[id]) return;
            used_edge[id] = true;
            result.push_back(id);
            covered[_edges[id].from] = true;
            covered[_edges[id].to] = true;
        };

        for (int v = 0; v < _n; v++) {
            if (v < _mate[v]) use_edge(_mate_edge[v]);
        }

        for (int v = 0; v < _n; v++) {
            if (covered[v]) continue;
            int id = -1;
            for (int edge_id : _adj[v]) {
                if (_edges[edge_id].alive) {
                    id = edge_id;
                    break;
                }
            }
            if (id == -1) return std::nullopt;
            use_edge(id);
        }

        return result;
    }
};

struct GeneralMatchingGraph {
    GeneralMatching matching;
    std::vector<int> original_edge_id;

    int original_edge(int edge_id) const {
        assert(0 <= edge_id && edge_id < int(original_edge_id.size()));
        return original_edge_id[edge_id];
    }
};

template <class T>
GeneralMatchingGraph make_general_matching(const Graph<T>& g) {
    GeneralMatchingGraph result;
    result.matching = GeneralMatching(g.size());
    for (const auto& e : g.edges()) {
        int id = result.matching.add_edge(e.from, e.to);
        if (int(result.original_edge_id.size()) <= id) result.original_edge_id.resize(id + 1);
        result.original_edge_id[id] = e.id;
    }
    return result;
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/general_weighted_matching.hpp"



#line 10 "graph/general_weighted_matching.hpp"
#include <type_traits>
#line 13 "graph/general_weighted_matching.hpp"

#line 15 "graph/general_weighted_matching.hpp"

namespace m1une {
namespace graph {
namespace internal {

// Primal-dual weighted blossom algorithm using Gabow's event queues.
// Reference: H. N. Gabow, "Data Structures for Weighted Matching and
// Extensions to b-matching and f-factors", 2016.
// Vertices 1..n are atoms; larger indices represent contracted blossoms.
template <class Cost, class TotalCost>
class WeightedBlossomSolver {
   public:
    using cost_type = Cost;
    using total_type = TotalCost;

   private:
    enum BlossomLabel : int {
        separated_label = -2,
        inner_label = -1,
        free_label = 0,
        outer_label = 1
    };
    static constexpr cost_type infinity = cost_type(1) << (sizeof(cost_type) * 8 - 2);

    template <class T>
    class MutableBinaryHeap {
       public:
        struct Node {
            bool operator<(const Node& rhs) const {
                if (value < rhs.value) {
                    return true;
                }
                if (rhs.value < value) {
                    return false;
                }
                return id < rhs.id;
            }
            T value;
            int id;
        };

        MutableBinaryHeap() = default;
        explicit MutableBinaryHeap(int capacity) : _size(0), _nodes(capacity + 1), _position(capacity, 0) {
        }

        bool empty() const {
            return _size == 0;
        }
        void clear() {
            while (_size > 0) {
                _position[_nodes[_size].id] = 0;
                _size--;
            }
        }
        T min() const {
            return _nodes[1].value;
        }
        int argmin() const {
            return _nodes[1].id;
        }
        void pop() {
            if (_size > 0) pop(1);
        }
        void erase(int id) {
            if (_position[id]) pop(_position[id]);
        }
        bool has(int id) const {
            return _position[id] != 0;
        }
        void update(int id, T v) {
            if (!has(id)) return push(id, v);
            bool up = (v < _nodes[_position[id]].value);
            _nodes[_position[id]].value = v;
            if (up) {
                up_heap(_position[id]);
            } else {
                down_heap(_position[id]);
            }
        }
        void decrease_key(int id, T v) {
            if (!has(id)) return push(id, v);
            if (v < _nodes[_position[id]].value) {
                _nodes[_position[id]].value = v;
                up_heap(_position[id]);
            }
        }
        void push(int id, T v) {
            _position[id] = ++_size;
            _nodes[_size] = {v, id};
            up_heap(_size);
        }

       private:
        void pop(int pos) {
            _position[_nodes[pos].id] = 0;
            if (pos == _size) {
                --_size;
                return;
            }
            bool up = (_nodes[_size].value < _nodes[pos].value);
            _nodes[pos] = _nodes[_size--];
            _position[_nodes[pos].id] = pos;
            if (up) {
                up_heap(pos);
            } else {
                down_heap(pos);
            }
        }
        void swap_node(int a, int b) {
            std::swap(_nodes[a], _nodes[b]);
            _position[_nodes[a].id] = a;
            _position[_nodes[b].id] = b;
        }
        void down_heap(int pos) {
            for (int current = pos;;) {
                int next = current;
                if (2 * current <= _size && _nodes[2 * current] < _nodes[next]) {
                    next = 2 * current;
                }
                if (2 * current + 1 <= _size && _nodes[2 * current + 1] < _nodes[next]) {
                    next = 2 * current + 1;
                }
                if (next == current) break;
                swap_node(current, next);
                current = next;
            }
        }
        void up_heap(int pos) {
            for (int current = pos; current > 1 && _nodes[current] < _nodes[current >> 1]; current >>= 1) {
                swap_node(current, current >> 1);
            }
        }
        int _size;
        std::vector<Node> _nodes;
        std::vector<int> _position;
    };

    template <class Key>
    class DisjointPairingHeaps {
       private:
        struct Node {
            Node() : key(), child(0), next(0), prev(-1) {
            }
            explicit Node(Key value) : key(value), child(0), next(0), prev(0) {
            }
            Key key;
            int child;
            int next;
            int prev;
        };

       public:
        DisjointPairingHeaps(int heap_count, int node_count) : _roots(heap_count), _nodes(node_count) {
        }

        void clear(int h) {
            if (_roots[h]) {
                clear_rec(_roots[h]);
                _roots[h] = 0;
            }
        }
        bool empty(int h) const {
            return !_roots[h];
        }
        bool used(int v) const {
            return _nodes[v].prev >= 0;
        }
        Key min(int h) const {
            return _nodes[_roots[h]].key;
        }
        void push(int h, int v, Key key) {
            _nodes[v] = Node(key);
            _roots[h] = merge(_roots[h], v);
        }
        void erase(int h, int v) {
            if (!used(v)) return;
            int w = two_pass_pairing(_nodes[v].child);
            if (!_nodes[v].prev) {
                _roots[h] = w;
            } else {
                cut(v);
                _roots[h] = merge(_roots[h], w);
            }
            _nodes[v].prev = -1;
        }
        void decrease_key(int h, int v, Key key) {
            if (!used(v)) return push(h, v, key);
            if (!_nodes[v].prev) {
                _nodes[v].key = key;
            } else {
                cut(v);
                _nodes[v].key = key;
                _roots[h] = merge(_roots[h], v);
            }
        }

       private:
        void clear_rec(int v) {
            for (; v; v = _nodes[v].next) {
                if (_nodes[v].child) clear_rec(_nodes[v].child);
                _nodes[v].prev = -1;
            }
        }

        inline void cut(int v) {
            auto& n = _nodes[v];
            int previous = n.prev;
            int next = n.next;
            auto& previous_node = _nodes[previous];
            if (previous_node.child == v) {
                previous_node.child = next;
            } else {
                previous_node.next = next;
            }
            _nodes[next].prev = previous;
            n.next = n.prev = 0;
        }

        int merge(int l, int r) {
            if (!l) return r;
            if (!r) return l;
            if (_nodes[l].key > _nodes[r].key) std::swap(l, r);
            int lc = _nodes[r].next = _nodes[l].child;
            _nodes[l].child = _nodes[lc].prev = r;
            return _nodes[r].prev = l;
        }

        int two_pass_pairing(int root) {
            if (!root) return 0;
            int a = root;
            root = 0;
            while (a) {
                int b = _nodes[a].next;
                int next_a = 0;
                _nodes[a].prev = _nodes[a].next = 0;
                if (b) {
                    next_a = _nodes[b].next;
                    _nodes[b].prev = _nodes[b].next = 0;
                }
                a = merge(a, b);
                _nodes[a].next = root;
                root = a;
                a = next_a;
            }
            int s = _nodes[root].next;
            _nodes[root].next = 0;
            while (s) {
                int t = _nodes[s].next;
                _nodes[s].next = 0;
                root = merge(root, s);
                s = t;
            }
            return root;
        }

       private:
        std::vector<int> _roots;
        std::vector<Node> _nodes;
    };

    template <class T>
    struct ReservablePriorityQueue : public std::priority_queue<T, std::vector<T>, std::greater<T>> {
        ReservablePriorityQueue() = default;
        explicit ReservablePriorityQueue(int capacity) {
            this->c.reserve(capacity);
        }
        T min() const {
            return this->top();
        }
        void clear() {
            this->c.clear();
        }
    };

    template <class T>
    struct FixedQueue {
        FixedQueue() = default;
        explicit FixedQueue(int capacity) : _head(0), _tail(0), _data(capacity) {
        }
        void enqueue(int u) {
            _data[_tail++] = u;
        }
        int dequeue() {
            return _data[_head++];
        }
        bool empty() const {
            return _head == _tail;
        }
        void clear() {
            _head = _tail = 0;
        }
        int _head = 0;
        int _tail = 0;
        std::vector<T> _data;
    };

   public:
    struct InputEdge {
        int from;
        int to;
        cost_type cost;
    };

   private:
    struct SolverEdge {
        int to;
        cost_type cost;
    };

    struct BlossomLink {
        int from;
        int to;
    };

    struct BlossomNode {
        struct CycleLink {
            int blossom;
            int vertex;
        };

        BlossomNode() = default;
        explicit BlossomNode(int vertex) : parent(0), size(1) {
            cycle[0] = cycle[1] = CycleLink{vertex, vertex};
        }

        int next_v() const {
            return cycle[0].vertex;
        }
        int next_b() const {
            return cycle[0].blossom;
        }
        int prev_v() const {
            return cycle[1].vertex;
        }
        int prev_b() const {
            return cycle[1].blossom;
        }

        int parent = 0;
        int size = 0;
        CycleLink cycle[2];
    };

    struct VertexEvent {
        VertexEvent() = default;
        VertexEvent(cost_type event_time, int vertex) : time(event_time), id(vertex) {
        }

        bool operator<(const VertexEvent& rhs) const {
            if (time < rhs.time) {
                return true;
            }
            if (rhs.time < time) {
                return false;
            }
            return id < rhs.id;
        }
        bool operator>(const VertexEvent& rhs) const {
            return rhs < *this;
        }

        cost_type time = cost_type();
        int id = 0;
    };

    struct EdgeEvent {
        EdgeEvent() = default;
        EdgeEvent(cost_type event_time, int from_, int to_) : time(event_time), from(from_), to(to_) {
        }

        bool operator<(const EdgeEvent& rhs) const {
            if (time < rhs.time) {
                return true;
            }
            if (time > rhs.time) {
                return false;
            }
            return std::make_pair(from, to) < std::make_pair(rhs.from, rhs.to);
        }
        bool operator>(const EdgeEvent& rhs) const {
            return rhs < *this;
        }

        cost_type time = cost_type();
        int from = 0;
        int to = 0;
    };

   public:
    WeightedBlossomSolver(int n, const std::vector<InputEdge>& input_edges)
        : _vertex_count(n),
          _blossom_count((n - 1) / 2),
          _state_count(n + _blossom_count + 1),
          _offset(n + 2),
          _edges(input_edges.size() * 2),
          _grow_heap(_state_count),
          _blossom_grow_heaps(_state_count, _state_count),
          _contract_heap(int(_edges.size())),
          _expand_heap(_state_count) {
        for (const InputEdge& edge : input_edges) {
            _offset[edge.from + 1]++;
            _offset[edge.to + 1]++;
        }
        for (int i = 1; i <= _vertex_count + 1; i++) _offset[i] += _offset[i - 1];
        for (const InputEdge& edge : input_edges) {
            _edges[_offset[edge.from]++] = SolverEdge{edge.to, edge.cost * 2};
            _edges[_offset[edge.to]++] = SolverEdge{edge.from, edge.cost * 2};
        }
        for (int i = _vertex_count + 1; i > 0; i--) _offset[i] = _offset[i - 1];
        _offset[0] = 0;
    }

    total_type solve(std::vector<std::pair<int, int>>& matching) {
        initialize_state();
        initialize_potentials();
        for (int vertex = 1; vertex <= _vertex_count; vertex++) {
            if (_mate[vertex] == 0) augment_from(vertex);
        }

        matching.clear();
        for (int vertex = 1; vertex <= _vertex_count; vertex++) {
            if (_mate[vertex] > vertex) matching.emplace_back(vertex, _mate[vertex]);
        }
        return compute_matching_weight();
    }

   private:
    total_type compute_matching_weight() const {
        total_type result = 0;
        for (int vertex = 1; vertex <= _vertex_count; vertex++) {
            if (_mate[vertex] > vertex) {
                cost_type best_cost = 0;
                for (int edge_id = _offset[vertex]; edge_id < _offset[vertex + 1]; edge_id++) {
                    if (_edges[edge_id].to == _mate[vertex]) {
                        best_cost = std::max(best_cost, _edges[edge_id].cost);
                    }
                }
                result += best_cost;
            }
        }
        return result >> 1;
    }

    total_type reduced_cost(int from, int to, const SolverEdge& edge) const {
        return total_type(_potential[from]) + _potential[to] - edge.cost;
    }

    void rematch(int vertex, int new_mate) {
        int old_mate = _mate[vertex];
        _mate[vertex] = new_mate;
        if (_mate[old_mate] != vertex) return;
        if (_tree_link[vertex].to == _surface[_tree_link[vertex].to]) {
            _mate[old_mate] = _tree_link[vertex].from;
            rematch(_mate[old_mate], old_mate);
        } else {
            int from = _tree_link[vertex].from;
            int to = _tree_link[vertex].to;
            rematch(from, to);
            rematch(to, from);
        }
    }

    void repair_matching(int blossom) {
        if (blossom <= _vertex_count) return;
        int child = _base[blossom];
        int first_vertex = _nodes[child].cycle[0].vertex;
        int first_neighbor = _nodes[child].cycle[0].blossom;
        int direction = (_nodes[first_neighbor].cycle[1].vertex == _mate[first_vertex]) ? 0 : 1;
        while (true) {
            int matched_vertex = _nodes[child].cycle[direction].vertex;
            int matched_child = _nodes[child].cycle[direction].blossom;
            if (_nodes[matched_child].cycle[1 ^ direction].vertex != _mate[matched_vertex]) break;
            repair_matching(child);
            repair_matching(matched_child);
            child = _nodes[matched_child].cycle[direction].blossom;
        }
        _base[blossom] = child;
        repair_matching(child);
        _mate[blossom] = _mate[child];
    }

    void reset_clock() {
        _time = 0;
        _vertex_event = {infinity, 0};
    }

    void reset_blossom(int blossom) {
        _label[blossom] = free_label;
        _tree_link[blossom].from = 0;
        _slack[blossom] = infinity;
        _lazy[blossom] = 0;
    }

    void reset_search_state() {
        _label[0] = free_label;
        _tree_link[0].from = 0;
        for (int vertex = 1; vertex <= _vertex_count; vertex++) {
            if (_label[vertex] == outer_label) {
                _potential[vertex] -= _time;
            } else {
                int blossom = _surface[vertex];
                _potential[vertex] += _lazy[blossom];
                if (_label[blossom] == inner_label) {
                    _potential[vertex] += _time - _created_at[blossom];
                }
            }
            reset_blossom(vertex);
        }
        int remaining_blossoms = _blossom_count - _unused_count;
        for (int blossom = _vertex_count + 1; remaining_blossoms > 0 && blossom < _state_count; blossom++) {
            if (_base[blossom] != blossom) {
                if (_surface[blossom] == blossom) {
                    repair_matching(blossom);
                    if (_label[blossom] == outer_label) {
                        _potential[blossom] += (_time - _created_at[blossom]) << 1;
                    } else if (_label[blossom] == inner_label) {
                        materialize_potential<inner_label>(blossom);
                    } else {
                        materialize_potential<free_label>(blossom);
                    }
                }
                _blossom_grow_heaps.clear(blossom);
                reset_blossom(blossom);
                remaining_blossoms--;
            }
        }

        _queue.clear();
        reset_clock();
        _grow_heap.clear();
        _contract_heap.clear();
        _expand_heap.clear();
    }

    void augment_from(int root) {
        if (_potential[root] == 0) return;
        link_blossom(_surface[root], {0, 0});
        make_outer(_surface[root], 0);
        for (bool augmented = false; !augmented;) {
            augmented = scan_tight_edges(root);
            if (augmented) break;
            augmented = advance_dual(root);
        }
        reset_search_state();
    }

    template <BlossomLabel target_label>
    cost_type materialize_potential(int blossom) {
        cost_type delta = _lazy[blossom];
        _lazy[blossom] = 0;
        if (target_label == inner_label) {
            cost_type elapsed = _time - _created_at[blossom];
            if (blossom > _vertex_count) _potential[blossom] -= elapsed << 1;
            delta += elapsed;
        }
        return delta;
    }

    template <BlossomLabel target_label>
    void update_grow_event(int from, int to, int to_blossom, cost_type slack) {
        if (slack >= _slack[to]) return;
        _slack[to] = slack;
        _best_from[to] = from;
        if (to == to_blossom) {
            if (target_label != inner_label) {
                _grow_heap.decrease_key(to, EdgeEvent(slack + _lazy[to], from, to));
            }
        } else {
            int to_group = _group[to];
            if (to_group != to) {
                if (slack >= _slack[to_group]) return;
                _slack[to_group] = slack;
            }
            _blossom_grow_heaps.decrease_key(to_blossom, to_group, EdgeEvent(slack, from, to));
            if (target_label == inner_label) return;
            EdgeEvent event = _blossom_grow_heaps.min(to_blossom);
            _grow_heap.decrease_key(to_blossom,
                                    EdgeEvent(event.time + _lazy[to_blossom], event.from, event.to));
        }
    }

    void activate_grow_event(int blossom) {
        if (blossom <= _vertex_count) {
            if (_slack[blossom] < infinity) {
                _grow_heap.push(blossom,
                                EdgeEvent(_slack[blossom] + _lazy[blossom], _best_from[blossom], blossom));
            }
        } else {
            if (_blossom_grow_heaps.empty(blossom)) return;
            EdgeEvent event = _blossom_grow_heaps.min(blossom);
            _grow_heap.push(blossom, EdgeEvent(event.time + _lazy[blossom], event.from, event.to));
        }
    }

    void swap_blossoms(int a, int b) {
        // b is a maximal blossom.
        std::swap(_base[a], _base[b]);
        if (_base[a] == a) _base[a] = b;
        std::swap(_heavy[a], _heavy[b]);
        if (_heavy[a] == a) _heavy[a] = b;
        std::swap(_tree_link[a], _tree_link[b]);
        std::swap(_mate[a], _mate[b]);
        std::swap(_potential[a], _potential[b]);
        std::swap(_lazy[a], _lazy[b]);
        std::swap(_created_at[a], _created_at[b]);
        for (int direction = 0; direction < 2; direction++) {
            int child = _nodes[a].cycle[direction].blossom;
            _nodes[child].cycle[1 ^ direction].blossom = b;
        }
        std::swap(_nodes[a], _nodes[b]);
    }

    void assign_surface(int blossom, int surface, int group) {
        _surface[blossom] = surface;
        _group[blossom] = group;
        if (blossom <= _vertex_count) return;
        for (int child = _base[blossom]; _surface[child] != surface; child = _nodes[child].next_b()) {
            assign_surface(child, surface, group);
        }
    }

    void merge_blossom_children(int blossom) {
        int largest_child = blossom;
        int largest_size = 1;
        int first_child = _base[blossom];
        for (int child = first_child;; child = _nodes[child].next_b()) {
            if (_nodes[child].size > largest_size) {
                largest_size = _nodes[child].size;
                largest_child = child;
            }
            if (_nodes[child].next_b() == first_child) break;
        }
        for (int child = first_child;; child = _nodes[child].next_b()) {
            if (child != largest_child) assign_surface(child, largest_child, child);
            if (_nodes[child].next_b() == first_child) break;
        }
        _group[largest_child] = largest_child;
        if (largest_size > 1) {
            _surface[blossom] = _heavy[blossom] = largest_child;
            swap_blossoms(largest_child, blossom);
        } else {
            _heavy[blossom] = 0;
        }
    }

    void contract_blossom(int x, int y, int edge_id) {
        int x_blossom = _surface[x];
        int y_blossom = _surface[y];
        assert(x_blossom != y_blossom);
        const int visit_mark = -(edge_id + 1);
        _tree_link[_surface[_mate[x_blossom]]].from = visit_mark;
        _tree_link[_surface[_mate[y_blossom]]].from = visit_mark;

        int lca = -1;
        while (true) {
            if (_mate[y_blossom] != 0) std::swap(x_blossom, y_blossom);
            x_blossom = lca = _surface[_tree_link[x_blossom].from];
            if (_tree_link[_surface[_mate[x_blossom]]].from == visit_mark) break;
            _tree_link[_surface[_mate[x_blossom]]].from = visit_mark;
        }

        const int blossom = _unused_blossoms[--_unused_count];
        assert(_unused_count >= 0);
        int tree_size = 0;
        for (int direction = 0; direction < 2; direction++) {
            for (int child = _surface[x]; child != lca;) {
                int matched_vertex = _mate[child];
                int matched_child = _surface[matched_vertex];
                int vertex = _mate[matched_vertex];
                int link_from = _tree_link[vertex].from;
                int link_to = _tree_link[vertex].to;
                tree_size += _nodes[child].size + _nodes[matched_child].size;
                _tree_link[matched_vertex] = {x, y};

                if (child > _vertex_count) {
                    _potential[child] += (_time - _created_at[child]) << 1;
                }
                if (matched_child > _vertex_count) _expand_heap.erase(matched_child);
                make_outer(matched_child, materialize_potential<inner_label>(matched_child));

                _nodes[child].cycle[direction] = {matched_child, matched_vertex};
                _nodes[matched_child].cycle[1 ^ direction] = {child, vertex};
                child = _surface[link_from];
                _nodes[matched_child].cycle[direction] = {child, link_from};
                _nodes[child].cycle[1 ^ direction] = {matched_child, link_to};
            }
            _nodes[_surface[x]].cycle[1 ^ direction] = {_surface[y], y};
            std::swap(x, y);
        }
        if (lca > _vertex_count) _potential[lca] += (_time - _created_at[lca]) << 1;
        _nodes[blossom].size = tree_size + _nodes[lca].size;
        _base[blossom] = lca;
        _tree_link[blossom] = _tree_link[lca];
        _mate[blossom] = _mate[lca];
        _label[blossom] = outer_label;
        _surface[blossom] = blossom;
        _created_at[blossom] = _time;
        _potential[blossom] = 0;
        _lazy[blossom] = 0;

        merge_blossom_children(blossom);
    }

    void link_blossom(int blossom, BlossomLink link) {
        _tree_link[blossom] = link;
        if (blossom <= _vertex_count) return;
        int first_child = _base[blossom];
        link_blossom(first_child, link);
        int previous_child = _nodes[first_child].prev_b();
        link = {_nodes[previous_child].next_v(), _nodes[first_child].prev_v()};
        for (int child = first_child;;) {
            int next_child = _nodes[child].next_b();
            if (next_child == first_child) break;
            link_blossom(next_child, link);
            BlossomLink next_link = {_nodes[next_child].prev_v(), _nodes[child].next_v()};
            child = _nodes[next_child].next_b();
            link_blossom(child, next_link);
        }
    }

    void make_outer(int blossom, cost_type delta) {
        _label[blossom] = outer_label;
        if (blossom > _vertex_count) {
            for (int child = _base[blossom]; _label[child] != outer_label; child = _nodes[child].next_b()) {
                make_outer(child, delta);
            }
        } else {
            _potential[blossom] += _time + delta;
            if (_potential[blossom] < _vertex_event.time) {
                _vertex_event = {_potential[blossom], blossom};
            }
            _queue.enqueue(blossom);
        }
    }

    bool grow_tree(int from, int to) {
        int inner_blossom = _surface[to];
        bool visited = (_label[inner_blossom] != free_label);
        if (!visited) link_blossom(inner_blossom, {0, 0});
        _label[inner_blossom] = inner_label;
        _created_at[inner_blossom] = _time;
        _grow_heap.erase(inner_blossom);
        if (to != inner_blossom) {
            _expand_heap.update(inner_blossom, _time + (_potential[inner_blossom] >> 1));
        }
        int matched_vertex = _mate[inner_blossom];
        if (matched_vertex == 0) {
            rematch(from, to);
            rematch(to, from);
            return true;
        }
        int outer_blossom = _surface[matched_vertex];
        if (!visited) {
            link_blossom(outer_blossom, {from, to});
        } else {
            _tree_link[outer_blossom] = _tree_link[matched_vertex] = {from, to};
        }
        make_outer(outer_blossom, materialize_potential<free_label>(outer_blossom));
        _created_at[outer_blossom] = _time;
        _grow_heap.erase(outer_blossom);
        return false;
    }

    void release_blossom(int blossom) {
        _unused_blossoms[_unused_count++] = blossom;
        _base[blossom] = blossom;
    }

    int recompute_slack(int blossom, int group) {
        if (blossom <= _vertex_count) {
            if (_slack[blossom] >= _slack[group]) return 0;
            _slack[group] = _slack[blossom];
            _best_from[group] = _best_from[blossom];
            return blossom;
        }
        int destination = 0;
        int first_child = _base[blossom];
        for (int child = first_child;; child = _nodes[child].next_b()) {
            int candidate = recompute_slack(child, group);
            if (candidate != 0) destination = candidate;
            if (_nodes[child].next_b() == first_child) break;
        }
        return destination;
    }

    void rebuild_components(int blossom, int surface, int group) {
        _surface[blossom] = surface;
        _group[blossom] = group;
        if (blossom <= _vertex_count) return;
        for (int child = _base[blossom]; _surface[child] != surface; child = _nodes[child].next_b()) {
            if (child == _heavy[blossom]) {
                rebuild_components(child, surface, group);
            } else {
                assign_surface(child, surface, child);
                int destination = 0;
                if (child > _vertex_count) {
                    _slack[child] = infinity;
                    destination = recompute_slack(child, child);
                } else if (_slack[child] < infinity) {
                    destination = child;
                }
                if (destination > 0) {
                    _blossom_grow_heaps.push(surface, child,
                                             EdgeEvent(_slack[child], _best_from[child], destination));
                }
            }
        }
    }

    void promote_largest_child(int blossom) {
        int largest_child = _heavy[blossom];
        cost_type delta = (_time - _created_at[blossom]) + _lazy[blossom];
        _lazy[blossom] = 0;
        int first_child = _base[blossom];
        for (int child = first_child;; child = _nodes[child].next_b()) {
            _created_at[child] = _time;
            _lazy[child] = delta;
            if (child != largest_child) {
                rebuild_components(child, child, child);
                _blossom_grow_heaps.erase(blossom, child);
            }
            if (_nodes[child].next_b() == first_child) break;
        }
        if (largest_child > 0) {
            swap_blossoms(largest_child, blossom);
            blossom = largest_child;
        }
        release_blossom(blossom);
    }

    void expand_blossom(int blossom) {
        int matched_vertex = _mate[_base[blossom]];
        promote_largest_child(blossom);
        BlossomLink old_link = _tree_link[matched_vertex];
        int old_base = _surface[_mate[matched_vertex]];
        int root = _surface[old_link.to];
        int direction = (_mate[root] == _nodes[root].cycle[0].vertex) ? 1 : 0;
        for (int child = _nodes[old_base].cycle[direction ^ 1].blossom; child != root;) {
            _label[child] = separated_label;
            activate_grow_event(child);
            child = _nodes[child].cycle[direction ^ 1].blossom;
            _label[child] = separated_label;
            activate_grow_event(child);
            child = _nodes[child].cycle[direction ^ 1].blossom;
        }
        for (int child = old_base;; child = _nodes[child].cycle[direction].blossom) {
            _label[child] = inner_label;
            int next_child = _nodes[child].cycle[direction].blossom;
            if (child == root) {
                _tree_link[_mate[child]] = old_link;
            } else {
                _tree_link[_mate[child]] = {_nodes[child].cycle[direction].vertex,
                                            _nodes[next_child].cycle[direction ^ 1].vertex};
            }
            _tree_link[_surface[_mate[child]]] = _tree_link[_mate[child]];
            if (child > _vertex_count) {
                if (_potential[child] == 0) {
                    expand_blossom(child);
                } else {
                    _expand_heap.push(child, _time + (_potential[child] >> 1));
                }
            }
            if (child == root) break;
            child = next_child;
            make_outer(next_child, materialize_potential<inner_label>(next_child));
        }
    }

    bool scan_tight_edges(int root) {
        while (!_queue.empty()) {
            int from = _queue.dequeue();
            int from_blossom = _surface[from];
            if (_potential[from] == _time) {
                if (from != root) rematch(from, 0);
                return true;
            }
            for (int edge_id = _offset[from]; edge_id < _offset[from + 1]; edge_id++) {
                const SolverEdge& edge = _edges[edge_id];
                int to = edge.to;
                int to_blossom = _surface[to];
                if (from_blossom == to_blossom) continue;
                BlossomLabel to_label = _label[to_blossom];
                if (to_label == outer_label) {
                    cost_type event_time = cost_type(reduced_cost(from, to, edge) >> 1);
                    if (event_time == _time) {
                        contract_blossom(from, to, edge_id);
                        from_blossom = _surface[from];
                    } else if (event_time < _vertex_event.time) {
                        _contract_heap.emplace(event_time, from, edge_id);
                    }
                } else {
                    total_type event_time = reduced_cost(from, to, edge);
                    if (event_time >= infinity) continue;
                    if (to_label != inner_label) {
                        if (cost_type(event_time) + _lazy[to_blossom] == _time) {
                            if (grow_tree(from, to)) return true;
                        } else {
                            update_grow_event<free_label>(from, to, to_blossom, cost_type(event_time));
                        }
                    } else if (_mate[from] != to) {
                        update_grow_event<inner_label>(from, to, to_blossom, cost_type(event_time));
                    }
                }
            }
        }
        return false;
    }

    bool advance_dual(int root) {
        cost_type rematch_time = _vertex_event.time;
        cost_type grow_time = infinity;
        if (!_grow_heap.empty()) grow_time = _grow_heap.min().time;

        cost_type contract_time = infinity;
        while (!_contract_heap.empty()) {
            EdgeEvent event = _contract_heap.min();
            int from = event.from;
            int to = _edges[event.to].to;
            if (_surface[from] != _surface[to]) {
                contract_time = event.time;
                break;
            } else {
                _contract_heap.pop();
            }
        }

        cost_type expand_time = infinity;
        if (!_expand_heap.empty()) expand_time = _expand_heap.min();

        cost_type next_time =
            std::min(std::min(rematch_time, grow_time), std::min(contract_time, expand_time));
        assert(_time <= next_time && next_time < infinity);
        _time = next_time;

        if (_time == _vertex_event.time) {
            int x = _vertex_event.id;
            if (x != root) rematch(x, 0);
            return true;
        }
        while (!_grow_heap.empty() && _grow_heap.min().time == _time) {
            int from = _grow_heap.min().from;
            int to = _grow_heap.min().to;
            if (grow_tree(from, to)) return true;
        }
        while (!_contract_heap.empty() && _contract_heap.min().time == _time) {
            int from = _contract_heap.min().from;
            int edge_id = _contract_heap.min().to;
            int to = _edges[edge_id].to;
            _contract_heap.pop();
            if (_surface[from] == _surface[to]) continue;
            contract_blossom(from, to, edge_id);
        }
        while (!_expand_heap.empty() && _expand_heap.min() == _time) {
            int blossom = _expand_heap.argmin();
            _expand_heap.pop();
            expand_blossom(blossom);
        }
        return false;
    }

   private:
    void initialize_state() {
        _queue = FixedQueue<int>(_vertex_count);
        _mate.assign(_state_count, 0);
        _tree_link.assign(_state_count, {0, 0});
        _label.assign(_state_count, free_label);
        _base.resize(_state_count);
        for (int state = 1; state < _state_count; state++) _base[state] = state;
        _surface.resize(_state_count);
        for (int state = 1; state < _state_count; state++) _surface[state] = state;

        _potential.resize(_state_count);
        _nodes.resize(_state_count);
        for (int state = 1; state < _state_count; state++) {
            _nodes[state] = BlossomNode(state);
        }

        _unused_blossoms.resize(_blossom_count);
        for (int i = 0; i < _blossom_count; i++) {
            _unused_blossoms[i] = _vertex_count + _blossom_count - i;
        }
        _unused_count = _blossom_count;

        reset_clock();
        _created_at.resize(_state_count);
        _slack.assign(_state_count, infinity);
        _best_from.assign(_state_count, 0);
        _heavy.assign(_state_count, 0);
        _lazy.assign(_state_count, 0);
        _group.resize(_state_count);
        for (int state = 0; state < _state_count; state++) _group[state] = state;
    }

    void initialize_potentials() {
        for (int vertex = 1; vertex <= _vertex_count; vertex++) {
            cost_type maximum_cost = 0;
            for (int edge_id = _offset[vertex]; edge_id < _offset[vertex + 1]; edge_id++) {
                maximum_cost = std::max(maximum_cost, _edges[edge_id].cost);
            }
            _potential[vertex] = maximum_cost >> 1;
        }
    }

    const int _vertex_count;
    const int _blossom_count;
    const int _state_count;
    std::vector<int> _offset;
    std::vector<SolverEdge> _edges;

    FixedQueue<int> _queue;
    std::vector<int> _mate;
    std::vector<int> _surface;
    std::vector<int> _base;
    std::vector<BlossomLink> _tree_link;
    std::vector<BlossomLabel> _label;
    std::vector<cost_type> _potential;

    std::vector<int> _unused_blossoms;
    int _unused_count;
    std::vector<BlossomNode> _nodes;

    // Heavy children and event queues keep each search phase at O(m log n).
    std::vector<int> _heavy;
    std::vector<int> _group;
    std::vector<cost_type> _created_at;
    std::vector<cost_type> _lazy;
    std::vector<cost_type> _slack;
    std::vector<int> _best_from;

    cost_type _time;
    VertexEvent _vertex_event;
    MutableBinaryHeap<EdgeEvent> _grow_heap;
    DisjointPairingHeaps<EdgeEvent> _blossom_grow_heaps;
    ReservablePriorityQueue<EdgeEvent> _contract_heap;
    MutableBinaryHeap<cost_type> _expand_heap;
};

}  // namespace internal

template <class Cost, class TotalCost = Cost>
struct GeneralWeightedMatching {
    static_assert(std::is_integral_v<Cost> && std::is_signed_v<Cost>);
    static_assert(std::is_integral_v<TotalCost> && std::is_signed_v<TotalCost>);

    struct Edge {
        int from;
        int to;
        Cost cost;
        int id;
        bool alive;

        int other(int vertex) const {
            assert(vertex == from || vertex == to);
            return from ^ to ^ vertex;
        }
    };

    struct Pair {
        int from;
        int to;
        Cost cost;
        int edge_id;
    };

   private:
    int _n;
    std::vector<Edge> _edges;
    std::vector<std::vector<int>> _adj;
    std::vector<int> _mate;
    std::vector<int> _mate_edge;
    TotalCost _matching_weight;
    bool _calculated;

    void invalidate() {
        _calculated = false;
    }

    void ensure_matching() {
        if (!_calculated) max_weight_matching();
    }

   public:
    GeneralWeightedMatching() : GeneralWeightedMatching(0) {
    }

    explicit GeneralWeightedMatching(int n)
        : _n(n), _adj(n), _mate(n, -1), _mate_edge(n, -1), _matching_weight(), _calculated(false) {
        assert(0 <= n);
    }

    int size() const {
        return _n;
    }

    int edge_count() const {
        return int(_edges.size());
    }

    int add_edge(int from, int to, Cost cost) {
        assert(0 <= from && from < _n);
        assert(0 <= to && to < _n);
        assert(from != to);
        assert(cost <= std::numeric_limits<Cost>::max() / Cost(2));
        int id = int(_edges.size());
        _edges.push_back(Edge{from, to, cost, id, true});
        _adj[from].push_back(id);
        _adj[to].push_back(id);
        invalidate();
        return id;
    }

    Edge get_edge(int id) const {
        assert(0 <= id && id < int(_edges.size()));
        return _edges[id];
    }

    std::vector<Edge> edges(bool include_inactive = false) const {
        std::vector<Edge> result;
        result.reserve(_edges.size());
        for (const Edge& edge : _edges) {
            if (include_inactive || edge.alive) result.push_back(edge);
        }
        return result;
    }

    void set_edge_alive(int id, bool alive) {
        assert(0 <= id && id < int(_edges.size()));
        _edges[id].alive = alive;
        invalidate();
    }

    void erase_edge(int id) {
        set_edge_alive(id, false);
    }

    void revive_edge(int id) {
        set_edge_alive(id, true);
    }

    bool is_edge_alive(int id) const {
        assert(0 <= id && id < int(_edges.size()));
        return _edges[id].alive;
    }

    TotalCost max_weight_matching() {
        using Solver = internal::WeightedBlossomSolver<Cost, TotalCost>;
        std::vector<typename Solver::InputEdge> input;
        input.reserve(_edges.size());
        for (const Edge& edge : _edges) {
            if (!edge.alive || edge.cost <= Cost()) continue;
            input.push_back(typename Solver::InputEdge{edge.from + 1, edge.to + 1, edge.cost});
        }

        Solver solver(_n, input);
        std::vector<std::pair<int, int>> vertex_pairs;
        solver.solve(vertex_pairs);

        _mate.assign(_n, -1);
        _mate_edge.assign(_n, -1);
        _matching_weight = TotalCost();
        for (auto [one_based_from, one_based_to] : vertex_pairs) {
            int from = one_based_from - 1;
            int to = one_based_to - 1;
            int best_edge = -1;
            for (int id : _adj[from]) {
                const Edge& edge = _edges[id];
                if (!edge.alive || edge.other(from) != to || edge.cost <= Cost()) continue;
                if (best_edge == -1 || _edges[best_edge].cost < edge.cost) best_edge = id;
            }
            assert(best_edge != -1);
            _mate[from] = to;
            _mate[to] = from;
            _mate_edge[from] = best_edge;
            _mate_edge[to] = best_edge;
            _matching_weight += static_cast<TotalCost>(_edges[best_edge].cost);
        }

        _calculated = true;
        return _matching_weight;
    }

    TotalCost matching_weight() {
        ensure_matching();
        return _matching_weight;
    }

    int matching_size() {
        ensure_matching();
        int result = 0;
        for (int vertex = 0; vertex < _n; vertex++) {
            if (vertex < _mate[vertex]) result++;
        }
        return result;
    }

    std::vector<int> mate() {
        ensure_matching();
        return _mate;
    }

    std::vector<int> mate_edge() {
        ensure_matching();
        return _mate_edge;
    }

    std::vector<Pair> matching() {
        ensure_matching();
        std::vector<Pair> result;
        for (int vertex = 0; vertex < _n; vertex++) {
            if (vertex < _mate[vertex]) {
                int id = _mate_edge[vertex];
                result.push_back(Pair{vertex, _mate[vertex], _edges[id].cost, id});
            }
        }
        return result;
    }
};

template <class Cost, class TotalCost = Cost>
struct GeneralWeightedMatchingGraph {
    GeneralWeightedMatching<Cost, TotalCost> matching;
    std::vector<int> original_edge_id;

    int original_edge(int edge_id) const {
        assert(0 <= edge_id && edge_id < int(original_edge_id.size()));
        return original_edge_id[edge_id];
    }
};

template <class T>
GeneralWeightedMatchingGraph<T> make_general_weighted_matching(const Graph<T>& graph) {
    GeneralWeightedMatchingGraph<T> result;
    result.matching = GeneralWeightedMatching<T>(graph.size());
    for (const auto& edge : graph.edges()) {
        int id = result.matching.add_edge(edge.from, edge.to, edge.cost);
        if (int(result.original_edge_id.size()) <= id) {
            result.original_edge_id.resize(id + 1);
        }
        result.original_edge_id[id] = edge.id;
    }
    return result;
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/grid.hpp"



#line 8 "graph/grid.hpp"

#line 10 "graph/grid.hpp"

namespace m1une {
namespace graph {

struct Grid {
   private:
    int _h;
    int _w;

   public:
    static constexpr std::array<int, 4> di4 = {-1, 0, 1, 0};
    static constexpr std::array<int, 4> dj4 = {0, 1, 0, -1};
    static constexpr std::array<int, 8> di8 = {-1, -1, -1, 0, 0, 1, 1, 1};
    static constexpr std::array<int, 8> dj8 = {-1, 0, 1, -1, 1, -1, 0, 1};

    Grid() : _h(0), _w(0) {}
    Grid(int h, int w) : _h(h), _w(w) {
        assert(0 <= h);
        assert(0 <= w);
    }

    int height() const {
        return _h;
    }

    int width() const {
        return _w;
    }

    int size() const {
        return _h * _w;
    }

    bool empty() const {
        return size() == 0;
    }

    bool inside(int i, int j) const {
        return 0 <= i && i < _h && 0 <= j && j < _w;
    }

    int id(int i, int j) const {
        assert(inside(i, j));
        return i * _w + j;
    }

    std::pair<int, int> pos(int v) const {
        assert(0 <= v && v < size());
        return {v / _w, v % _w};
    }

    std::vector<std::pair<int, int>> adj4(int i, int j) const {
        assert(inside(i, j));
        std::vector<std::pair<int, int>> result;
        result.reserve(4);
        for (int k = 0; k < 4; k++) {
            int ni = i + di4[k], nj = j + dj4[k];
            if (inside(ni, nj)) result.emplace_back(ni, nj);
        }
        return result;
    }

    std::vector<std::pair<int, int>> adj8(int i, int j) const {
        assert(inside(i, j));
        std::vector<std::pair<int, int>> result;
        result.reserve(8);
        for (int k = 0; k < 8; k++) {
            int ni = i + di8[k], nj = j + dj8[k];
            if (inside(ni, nj)) result.emplace_back(ni, nj);
        }
        return result;
    }

    std::vector<int> adj4_ids(int v) const {
        auto [i, j] = pos(v);
        std::vector<int> result;
        result.reserve(4);
        for (auto [ni, nj] : adj4(i, j)) result.push_back(id(ni, nj));
        return result;
    }

    std::vector<int> adj8_ids(int v) const {
        auto [i, j] = pos(v);
        std::vector<int> result;
        result.reserve(8);
        for (auto [ni, nj] : adj8(i, j)) result.push_back(id(ni, nj));
        return result;
    }

    Graph<int> graph4() const {
        return graph4([](int, int) { return true; });
    }

    Graph<int> graph8() const {
        return graph8([](int, int) { return true; });
    }

    template <class Passable>
    Graph<int> graph4(Passable passable) const {
        Graph<int> g(size());
        for (int i = 0; i < _h; i++) {
            for (int j = 0; j < _w; j++) {
                if (!passable(i, j)) continue;
                int v = id(i, j);
                for (auto [ni, nj] : adj4(i, j)) {
                    if (!passable(ni, nj)) continue;
                    int to = id(ni, nj);
                    if (v < to) g.add_edge(v, to);
                }
            }
        }
        return g;
    }

    template <class Passable>
    Graph<int> graph8(Passable passable) const {
        Graph<int> g(size());
        for (int i = 0; i < _h; i++) {
            for (int j = 0; j < _w; j++) {
                if (!passable(i, j)) continue;
                int v = id(i, j);
                for (auto [ni, nj] : adj8(i, j)) {
                    if (!passable(ni, nj)) continue;
                    int to = id(ni, nj);
                    if (v < to) g.add_edge(v, to);
                }
            }
        }
        return g;
    }
};

}  // namespace graph
}  // namespace m1une


#line 1 "graph/kruskal.hpp"



#line 6 "graph/kruskal.hpp"

#line 9 "graph/kruskal.hpp"

namespace m1une {
namespace graph {

template <class T>
struct MinimumSpanningForest {
    T cost;
    std::vector<Edge<T>> edges;
    int components;

    bool is_spanning_tree(int n) const {
        return components <= 1 && int(edges.size()) == std::max(0, n - 1);
    }
};

template <class T>
MinimumSpanningForest<T> kruskal(const Graph<T>& g) {
    int n = g.size();
    auto edges = g.edges();
    std::sort(edges.begin(), edges.end(), [](const auto& a, const auto& b) {
        return a.cost < b.cost;
    });

    m1une::ds::Dsu dsu(n);
    MinimumSpanningForest<T> result;
    result.cost = T(0);
    result.components = n;

    for (const auto& e : edges) {
        if (dsu.same(e.from, e.to)) continue;
        dsu.merge(e.from, e.to);
        result.cost += e.cost;
        result.edges.push_back(e);
        result.components--;
    }

    return result;
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/lowlink.hpp"



#line 6 "graph/lowlink.hpp"

#line 8 "graph/lowlink.hpp"

namespace m1une {
namespace graph {

template <class T>
struct LowLinkResult {
    std::vector<int> ord;
    std::vector<int> low;
    std::vector<int> articulation;
    std::vector<Edge<T>> bridges;
    std::vector<int> bridge_ids;
};

template <class T>
LowLinkResult<T> lowlink(const Graph<T>& g) {
    int n = g.size();
    LowLinkResult<T> result;
    result.ord.assign(n, -1);
    result.low.assign(n, -1);
    int now = 0;

    auto dfs = [&](auto self, int v, int parent_edge) -> void {
        result.ord[v] = result.low[v] = now++;
        int child_count = 0;
        bool is_articulation = false;

        for (const auto& e : g[v]) {
            if (!e.alive) continue;
            if (e.id == parent_edge) continue;
            int to = e.to;
            if (result.ord[to] == -1) {
                child_count++;
                self(self, to, e.id);
                result.low[v] = std::min(result.low[v], result.low[to]);
                if (parent_edge != -1 && result.ord[v] <= result.low[to]) is_articulation = true;
                if (result.ord[v] < result.low[to]) {
                    result.bridges.push_back(e);
                    result.bridge_ids.push_back(e.id);
                }
            } else {
                result.low[v] = std::min(result.low[v], result.ord[to]);
            }
        }

        if (parent_edge == -1 && child_count >= 2) is_articulation = true;
        if (is_articulation) result.articulation.push_back(v);
    };

    for (int v = 0; v < n; v++) {
        if (result.ord[v] == -1) dfs(dfs, v, -1);
    }
    std::sort(result.articulation.begin(), result.articulation.end());
    std::sort(result.bridge_ids.begin(), result.bridge_ids.end());
    return result;
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/maximum_clique.hpp"



#line 7 "graph/maximum_clique.hpp"

#line 9 "graph/maximum_clique.hpp"

namespace m1une {
namespace graph {

struct MaximumCliqueResult {
    std::vector<int> vertices;

    int size() const {
        return int(vertices.size());
    }

    bool empty() const {
        return vertices.empty();
    }
};

struct MaximumIndependentSetResult {
    std::vector<int> vertices;

    int size() const {
        return int(vertices.size());
    }

    bool empty() const {
        return vertices.empty();
    }
};

struct MinimumVertexCoverResult {
    std::vector<int> vertices;

    int size() const {
        return int(vertices.size());
    }

    bool empty() const {
        return vertices.empty();
    }
};

namespace detail {

struct MaximumIndependentSetBranching {
    int n;
    std::vector<std::vector<char>> adjacent;
    std::vector<std::vector<int>> graph;

    explicit MaximumIndependentSetBranching(const std::vector<std::vector<char>>& adjacent_)
        : n(int(adjacent_.size())), adjacent(adjacent_), graph(n) {
        for (int v = 0; v < n; v++) {
            for (int to = 0; to < n; to++) {
                if (adjacent[v][to]) graph[v].push_back(to);
            }
        }
    }

    std::vector<int> solve_path(const std::vector<int>& order) const {
        int m = int(order.size());
        if (m == 0) return {};

        std::vector<int> dp0(m, 0), dp1(m, 0);
        dp1[0] = 1;
        for (int i = 1; i < m; i++) {
            dp0[i] = std::max(dp0[i - 1], dp1[i - 1]);
            dp1[i] = dp0[i - 1] + 1;
        }

        std::vector<int> result;
        int state = (dp1[m - 1] > dp0[m - 1] ? 1 : 0);
        for (int i = m - 1; i >= 0; i--) {
            if (state == 1) {
                result.push_back(order[i]);
                state = 0;
            } else if (i > 0) {
                state = (dp1[i - 1] > dp0[i - 1] ? 1 : 0);
            }
        }
        return result;
    }

    std::vector<int> solve_cycle(const std::vector<int>& order) const {
        int m = int(order.size());
        if (m == 0) return {};
        if (m == 1) return {order[0]};

        std::vector<int> without_first(order.begin() + 1, order.end());
        auto result_without = solve_path(without_first);

        std::vector<int> result_with = {order[0]};
        if (m >= 4) {
            std::vector<int> middle(order.begin() + 2, order.end() - 1);
            auto middle_result = solve_path(middle);
            result_with.insert(result_with.end(), middle_result.begin(), middle_result.end());
        }

        return (result_with.size() > result_without.size() ? result_with : result_without);
    }

    std::vector<int> solve_degree_at_most_two(const std::vector<char>& active,
                                              const std::vector<int>& degree) const {
        std::vector<int> result;
        std::vector<char> visited(n, false);

        for (int s = 0; s < n; s++) {
            if (!active[s] || visited[s]) continue;

            std::vector<int> component;
            std::vector<int> stack = {s};
            visited[s] = true;
            for (int it = 0; it < int(stack.size()); it++) {
                int v = stack[it];
                component.push_back(v);
                for (int to : graph[v]) {
                    if (!active[to] || visited[to]) continue;
                    visited[to] = true;
                    stack.push_back(to);
                }
            }

            if (component.size() == 1) {
                result.push_back(component[0]);
                continue;
            }

            int endpoint = -1;
            for (int v : component) {
                if (degree[v] <= 1) {
                    endpoint = v;
                    break;
                }
            }

            std::vector<int> order;
            if (endpoint != -1) {
                int prev = -1, cur = endpoint;
                while (cur != -1) {
                    order.push_back(cur);
                    int next = -1;
                    for (int to : graph[cur]) {
                        if (active[to] && to != prev) {
                            next = to;
                            break;
                        }
                    }
                    prev = cur;
                    cur = next;
                }
                auto part = solve_path(order);
                result.insert(result.end(), part.begin(), part.end());
            } else {
                int start = component[0];
                int first = -1;
                for (int to : graph[start]) {
                    if (active[to]) {
                        first = to;
                        break;
                    }
                }
                assert(first != -1);

                order.push_back(start);
                int prev = start, cur = first;
                while (cur != start) {
                    order.push_back(cur);
                    int next = -1;
                    for (int to : graph[cur]) {
                        if (active[to] && to != prev) {
                            next = to;
                            break;
                        }
                    }
                    assert(next != -1);
                    prev = cur;
                    cur = next;
                }
                auto part = solve_cycle(order);
                result.insert(result.end(), part.begin(), part.end());
            }
        }

        return result;
    }

    std::vector<int> solve(std::vector<char> active) const {
        int active_count = 0;
        int max_degree = -1;
        int branch_vertex = -1;
        std::vector<int> degree(n, 0);

        for (int v = 0; v < n; v++) {
            if (!active[v]) continue;
            active_count++;
            for (int to : graph[v]) {
                if (active[to]) degree[v]++;
            }
            if (degree[v] > max_degree) {
                max_degree = degree[v];
                branch_vertex = v;
            }
        }

        if (active_count == 0) return {};
        if (max_degree <= 2) {
            auto result = solve_degree_at_most_two(active, degree);
            std::sort(result.begin(), result.end());
            return result;
        }

        auto without = active;
        without[branch_vertex] = false;
        auto result_without = solve(without);

        auto with = active;
        with[branch_vertex] = false;
        for (int to : graph[branch_vertex]) with[to] = false;
        auto result_with = solve(with);
        result_with.push_back(branch_vertex);

        auto result = (result_with.size() > result_without.size() ? result_with : result_without);
        std::sort(result.begin(), result.end());
        return result;
    }

    std::vector<int> solve() const {
        std::vector<char> active(n, true);
        return solve(active);
    }
};

template <class T>
std::vector<std::vector<char>> undirected_adjacency_matrix(const Graph<T>& g) {
    int n = g.size();
    std::vector<std::vector<char>> adjacent(n, std::vector<char>(n, false));
    for (const auto& e : g.edges()) {
        if (e.from == e.to) continue;
        adjacent[e.from][e.to] = true;
        adjacent[e.to][e.from] = true;
    }
    return adjacent;
}

std::vector<std::vector<char>> complement_adjacency_matrix(const std::vector<std::vector<char>>& adjacent) {
    int n = int(adjacent.size());
    std::vector<std::vector<char>> complement(n, std::vector<char>(n, false));
    for (int i = 0; i < n; i++) {
        for (int j = i + 1; j < n; j++) {
            if (adjacent[i][j]) continue;
            complement[i][j] = true;
            complement[j][i] = true;
        }
    }
    return complement;
}

}  // namespace detail

template <class T>
bool is_clique(const Graph<T>& g, const std::vector<int>& vertices) {
    auto adjacent = detail::undirected_adjacency_matrix(g);
    for (int v : vertices) {
        assert(0 <= v && v < g.size());
    }
    for (int i = 0; i < int(vertices.size()); i++) {
        for (int j = i + 1; j < int(vertices.size()); j++) {
            if (!adjacent[vertices[i]][vertices[j]]) return false;
        }
    }
    return true;
}

template <class T>
bool is_independent_set(const Graph<T>& g, const std::vector<int>& vertices) {
    auto adjacent = detail::undirected_adjacency_matrix(g);
    for (int v : vertices) {
        assert(0 <= v && v < g.size());
    }
    for (int i = 0; i < int(vertices.size()); i++) {
        for (int j = i + 1; j < int(vertices.size()); j++) {
            if (adjacent[vertices[i]][vertices[j]]) return false;
        }
    }
    return true;
}

template <class T>
bool is_vertex_cover(const Graph<T>& g, const std::vector<int>& vertices) {
    std::vector<char> selected(g.size(), false);
    for (int v : vertices) {
        assert(0 <= v && v < g.size());
        selected[v] = true;
    }
    for (const auto& e : g.edges()) {
        if (e.from == e.to) continue;
        if (!selected[e.from] && !selected[e.to]) return false;
    }
    return true;
}

template <class T>
MaximumCliqueResult maximum_clique(const Graph<T>& g) {
    auto adjacent = detail::undirected_adjacency_matrix(g);
    auto complement = detail::complement_adjacency_matrix(adjacent);
    detail::MaximumIndependentSetBranching solver(complement);
    return MaximumCliqueResult{solver.solve()};
}

template <class T>
int maximum_clique_size(const Graph<T>& g) {
    return maximum_clique(g).size();
}

template <class T>
MaximumIndependentSetResult maximum_independent_set(const Graph<T>& g) {
    auto adjacent = detail::undirected_adjacency_matrix(g);
    detail::MaximumIndependentSetBranching solver(adjacent);
    return MaximumIndependentSetResult{solver.solve()};
}

template <class T>
int maximum_independent_set_size(const Graph<T>& g) {
    return maximum_independent_set(g).size();
}

template <class T>
MinimumVertexCoverResult minimum_vertex_cover(const Graph<T>& g) {
    auto independent = maximum_independent_set(g);
    std::vector<char> in_independent(g.size(), false);
    for (int v : independent.vertices) in_independent[v] = true;

    MinimumVertexCoverResult result;
    for (int v = 0; v < g.size(); v++) {
        if (!in_independent[v]) result.vertices.push_back(v);
    }
    return result;
}

template <class T>
int minimum_vertex_cover_size(const Graph<T>& g) {
    return minimum_vertex_cover(g).size();
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/matrix_tree_theorem.hpp"



#line 7 "graph/matrix_tree_theorem.hpp"

#line 1 "math/matrix/linear_algebra.hpp"



#line 7 "math/matrix/linear_algebra.hpp"

#line 1 "math/matrix/matrix.hpp"



#line 9 "math/matrix/matrix.hpp"

namespace m1une {
namespace matrix {

template <class T>
class Matrix {
   private:
    int _rows;
    int _cols;
    std::vector<T> _data;

    static std::size_t storage_size(int rows, int cols) {
        assert(rows >= 0);
        assert(cols >= 0);
        return std::size_t(rows) * std::size_t(cols);
    }

   public:
    using value_type = T;

    Matrix() : _rows(0), _cols(0) {}

    Matrix(int rows, int cols, const T& value = T())
        : _rows(rows), _cols(cols), _data(storage_size(rows, cols), value) {}

    Matrix(int rows, int cols, std::vector<T> values)
        : _rows(rows), _cols(cols), _data(std::move(values)) {
        assert(rows >= 0);
        assert(cols >= 0);
        assert(_data.size() == std::size_t(rows) * std::size_t(cols));
    }

    explicit Matrix(const std::vector<std::vector<T>>& values)
        : _rows(int(values.size())), _cols(values.empty() ? 0 : int(values[0].size())),
          _data(storage_size(_rows, _cols)) {
        for (int row = 0; row < _rows; row++) {
            assert(int(values[std::size_t(row)].size()) == _cols);
            for (int col = 0; col < _cols; col++) {
                (*this)[row][col] = values[std::size_t(row)][std::size_t(col)];
            }
        }
    }

    int rows() const {
        return _rows;
    }

    int cols() const {
        return _cols;
    }

    bool empty() const {
        return _rows == 0 || _cols == 0;
    }

    std::vector<T>& data() {
        return _data;
    }

    const std::vector<T>& data() const {
        return _data;
    }

    T* operator[](int row) {
        assert(0 <= row && row < _rows);
        return _data.data() + std::size_t(row) * std::size_t(_cols);
    }

    const T* operator[](int row) const {
        assert(0 <= row && row < _rows);
        return _data.data() + std::size_t(row) * std::size_t(_cols);
    }

    T& operator()(int row, int col) {
        assert(0 <= col && col < _cols);
        return (*this)[row][col];
    }

    const T& operator()(int row, int col) const {
        assert(0 <= col && col < _cols);
        return (*this)[row][col];
    }

    static Matrix identity(int size) {
        assert(size >= 0);
        Matrix result(size, size);
        for (int i = 0; i < size; i++) result[i][i] = T(1);
        return result;
    }

    Matrix transposed() const {
        Matrix result(_cols, _rows);
        for (int row = 0; row < _rows; row++) {
            for (int col = 0; col < _cols; col++) {
                result[col][row] = (*this)[row][col];
            }
        }
        return result;
    }

    void swap_rows(int first, int second) {
        assert(0 <= first && first < _rows);
        assert(0 <= second && second < _rows);
        if (first == second) return;
        for (int col = 0; col < _cols; col++) {
            std::swap((*this)[first][col], (*this)[second][col]);
        }
    }

    Matrix& operator+=(const Matrix& rhs) {
        assert(_rows == rhs._rows && _cols == rhs._cols);
        for (std::size_t i = 0; i < _data.size(); i++) _data[i] += rhs._data[i];
        return *this;
    }

    Matrix& operator-=(const Matrix& rhs) {
        assert(_rows == rhs._rows && _cols == rhs._cols);
        for (std::size_t i = 0; i < _data.size(); i++) _data[i] -= rhs._data[i];
        return *this;
    }

    Matrix& operator*=(const T& scalar) {
        for (T& value : _data) value *= scalar;
        return *this;
    }

    Matrix& operator/=(const T& scalar) {
        for (T& value : _data) value /= scalar;
        return *this;
    }

    Matrix& operator*=(const Matrix& rhs) {
        return *this = *this * rhs;
    }

    Matrix operator+() const {
        return *this;
    }

    Matrix operator-() const {
        Matrix result = *this;
        for (T& value : result._data) value = T() - value;
        return result;
    }

    friend Matrix operator+(Matrix lhs, const Matrix& rhs) {
        return lhs += rhs;
    }

    friend Matrix operator-(Matrix lhs, const Matrix& rhs) {
        return lhs -= rhs;
    }

    friend Matrix operator*(Matrix lhs, const T& rhs) {
        return lhs *= rhs;
    }

    friend Matrix operator*(const T& lhs, Matrix rhs) {
        return rhs *= lhs;
    }

    friend Matrix operator/(Matrix lhs, const T& rhs) {
        return lhs /= rhs;
    }

    friend Matrix operator*(const Matrix& lhs, const Matrix& rhs) {
        assert(lhs._cols == rhs._rows);
        Matrix result(lhs._rows, rhs._cols);
        for (int row = 0; row < lhs._rows; row++) {
            T* output = result[row];
            for (int middle = 0; middle < lhs._cols; middle++) {
                const T coefficient = lhs[row][middle];
                if (coefficient == T()) continue;
                const T* input = rhs[middle];
                for (int col = 0; col < rhs._cols; col++) {
                    output[col] += coefficient * input[col];
                }
            }
        }
        return result;
    }

    friend std::vector<T> operator*(const Matrix& lhs, const std::vector<T>& rhs) {
        assert(lhs._cols == int(rhs.size()));
        std::vector<T> result(std::size_t(lhs._rows));
        for (int row = 0; row < lhs._rows; row++) {
            T value = T();
            for (int col = 0; col < lhs._cols; col++) {
                value += lhs[row][col] * rhs[std::size_t(col)];
            }
            result[std::size_t(row)] = value;
        }
        return result;
    }

    friend std::vector<T> operator*(const std::vector<T>& lhs, const Matrix& rhs) {
        assert(int(lhs.size()) == rhs._rows);
        std::vector<T> result(std::size_t(rhs._cols));
        for (int row = 0; row < rhs._rows; row++) {
            if (lhs[std::size_t(row)] == T()) continue;
            for (int col = 0; col < rhs._cols; col++) {
                result[std::size_t(col)] += lhs[std::size_t(row)] * rhs[row][col];
            }
        }
        return result;
    }

    bool operator==(const Matrix& rhs) const {
        return _rows == rhs._rows && _cols == rhs._cols && _data == rhs._data;
    }

    bool operator!=(const Matrix& rhs) const {
        return !(*this == rhs);
    }

    Matrix pow(std::uint64_t exponent) const {
        assert(_rows == _cols);
        Matrix result = identity(_rows);
        Matrix base = *this;
        while (exponent > 0) {
            if (exponent & 1) result *= base;
            exponent >>= 1;
            if (exponent > 0) base *= base;
        }
        return result;
    }
};

}  // namespace matrix
}  // namespace m1une


#line 9 "math/matrix/linear_algebra.hpp"

namespace m1une {
namespace matrix {

template <class T>
constexpr T default_epsilon() {
    if constexpr (std::is_floating_point_v<T>) {
        return T(1e-10);
    } else {
        return T();
    }
}

namespace detail {

template <class T>
T matrix_abs(T value) {
    return value < T() ? T() - value : value;
}

template <class T>
bool is_zero(const T& value, const T& eps) {
    if constexpr (std::is_floating_point_v<T>) {
        return matrix_abs(value) <= eps;
    } else {
        (void)eps;
        return value == T();
    }
}

template <class T>
int choose_pivot(const Matrix<T>& matrix, int first_row, int col, const T& eps) {
    int pivot = -1;
    if constexpr (std::is_floating_point_v<T>) {
        for (int row = first_row; row < matrix.rows(); row++) {
            if (is_zero(matrix[row][col], eps)) continue;
            if (pivot == -1 || matrix_abs(matrix[pivot][col]) < matrix_abs(matrix[row][col])) {
                pivot = row;
            }
        }
    } else {
        for (int row = first_row; row < matrix.rows(); row++) {
            if (!is_zero(matrix[row][col], eps)) {
                pivot = row;
                break;
            }
        }
    }
    return pivot;
}

template <class T>
std::vector<int> row_reduce(Matrix<T>& matrix, int pivot_col_limit, const T& eps,
                            bool reduced) {
    std::vector<int> pivot_columns;
    int pivot_row = 0;
    for (int col = 0; col < pivot_col_limit && pivot_row < matrix.rows(); col++) {
        int pivot = choose_pivot(matrix, pivot_row, col, eps);
        if (pivot == -1) continue;
        matrix.swap_rows(pivot_row, pivot);

        const T pivot_value = matrix[pivot_row][col];
        if (reduced) {
            for (int j = col; j < matrix.cols(); j++) matrix[pivot_row][j] /= pivot_value;
        }

        const int first_row = reduced ? 0 : pivot_row + 1;
        for (int row = first_row; row < matrix.rows(); row++) {
            if (row == pivot_row || is_zero(matrix[row][col], eps)) continue;
            T factor = matrix[row][col];
            if (!reduced) factor /= pivot_value;
            matrix[row][col] = T();
            for (int j = col + 1; j < matrix.cols(); j++) {
                matrix[row][j] -= factor * matrix[pivot_row][j];
            }
        }

        pivot_columns.push_back(col);
        pivot_row++;
    }

    if constexpr (std::is_floating_point_v<T>) {
        for (T& value : matrix.data()) {
            if (is_zero(value, eps)) value = T();
        }
    }
    return pivot_columns;
}

}  // namespace detail

template <class T>
struct RowReduction {
    Matrix<T> matrix;
    std::vector<int> pivot_columns;

    int rank() const {
        return int(pivot_columns.size());
    }
};

template <class T>
RowReduction<T> reduced_row_echelon_form(Matrix<T> matrix,
                                         T eps = default_epsilon<T>()) {
    RowReduction<T> result;
    result.pivot_columns = detail::row_reduce(matrix, matrix.cols(), eps, true);
    result.matrix = std::move(matrix);
    return result;
}

template <class T>
int matrix_rank(Matrix<T> matrix, T eps = default_epsilon<T>()) {
    return int(detail::row_reduce(matrix, matrix.cols(), eps, false).size());
}

template <class T>
T determinant(Matrix<T> matrix, T eps = default_epsilon<T>()) {
    assert(matrix.rows() == matrix.cols());
    const int size = matrix.rows();
    T result = T(1);
    bool negate = false;

    for (int col = 0; col < size; col++) {
        int pivot = detail::choose_pivot(matrix, col, col, eps);
        if (pivot == -1) return T();
        if (pivot != col) {
            matrix.swap_rows(pivot, col);
            negate = !negate;
        }

        const T pivot_value = matrix[col][col];
        result *= pivot_value;
        for (int row = col + 1; row < size; row++) {
            if (detail::is_zero(matrix[row][col], eps)) continue;
            const T factor = matrix[row][col] / pivot_value;
            matrix[row][col] = T();
            for (int j = col + 1; j < size; j++) {
                matrix[row][j] -= factor * matrix[col][j];
            }
        }
    }
    return negate ? T() - result : result;
}

template <class T>
std::optional<Matrix<T>> inverse(const Matrix<T>& matrix,
                                 T eps = default_epsilon<T>()) {
    assert(matrix.rows() == matrix.cols());
    const int size = matrix.rows();
    Matrix<T> augmented(size, size * 2);
    for (int row = 0; row < size; row++) {
        for (int col = 0; col < size; col++) {
            augmented[row][col] = matrix[row][col];
        }
        augmented[row][size + row] = T(1);
    }

    const std::vector<int> pivots = detail::row_reduce(augmented, size, eps, true);
    if (int(pivots.size()) != size) return std::nullopt;

    Matrix<T> result(size, size);
    for (int row = 0; row < size; row++) {
        for (int col = 0; col < size; col++) {
            result[row][col] = augmented[row][size + col];
        }
    }
    return result;
}

template <class T>
struct LinearSystemResult {
    bool consistent = false;
    std::vector<T> particular_solution;
    std::vector<std::vector<T>> nullspace_basis;
    std::vector<int> pivot_columns;

    int rank() const {
        return int(pivot_columns.size());
    }

    int nullity() const {
        return consistent ? int(nullspace_basis.size()) : 0;
    }

    bool has_unique_solution() const {
        return consistent && nullspace_basis.empty();
    }
};

template <class T>
LinearSystemResult<T> solve_linear_system(const Matrix<T>& coefficients,
                                          const std::vector<T>& constants,
                                          T eps = default_epsilon<T>()) {
    assert(coefficients.rows() == int(constants.size()));
    const int equation_count = coefficients.rows();
    const int variable_count = coefficients.cols();
    Matrix<T> augmented(equation_count, variable_count + 1);
    for (int row = 0; row < equation_count; row++) {
        for (int col = 0; col < variable_count; col++) {
            augmented[row][col] = coefficients[row][col];
        }
        augmented[row][variable_count] = constants[std::size_t(row)];
    }

    LinearSystemResult<T> result;
    result.pivot_columns =
        detail::row_reduce(augmented, variable_count, eps, true);

    for (int row = result.rank(); row < equation_count; row++) {
        bool zero_left = true;
        for (int col = 0; col < variable_count; col++) {
            if (!detail::is_zero(augmented[row][col], eps)) {
                zero_left = false;
                break;
            }
        }
        if (zero_left && !detail::is_zero(augmented[row][variable_count], eps)) {
            return result;
        }
    }

    result.consistent = true;
    result.particular_solution.assign(std::size_t(variable_count), T());
    std::vector<bool> is_pivot(std::size_t(variable_count), false);
    for (int row = 0; row < result.rank(); row++) {
        const int col = result.pivot_columns[std::size_t(row)];
        is_pivot[std::size_t(col)] = true;
        result.particular_solution[std::size_t(col)] = augmented[row][variable_count];
    }

    for (int free_col = 0; free_col < variable_count; free_col++) {
        if (is_pivot[std::size_t(free_col)]) continue;
        std::vector<T> direction(static_cast<std::size_t>(variable_count));
        direction[std::size_t(free_col)] = T(1);
        for (int row = 0; row < result.rank(); row++) {
            const int pivot_col = result.pivot_columns[std::size_t(row)];
            direction[std::size_t(pivot_col)] = T() - augmented[row][free_col];
        }
        result.nullspace_basis.push_back(std::move(direction));
    }
    return result;
}

}  // namespace matrix
}  // namespace m1une


#line 10 "graph/matrix_tree_theorem.hpp"

namespace m1une {
namespace graph {

namespace matrix_tree_detail {

inline int minor_index(int vertex, int removed) {
    assert(vertex != removed);
    return vertex < removed ? vertex : vertex - 1;
}

template <class Weight>
void assert_edge_incidence(const Graph<Weight>& graph, int expected) {
#ifndef NDEBUG
    std::vector<int> incidence(graph.edge_count(), 0);
    for (int vertex = 0; vertex < graph.size(); vertex++) {
        for (const Edge<Weight>& edge : graph[vertex]) {
            if (!edge.alive) continue;
            assert(0 <= edge.id && edge.id < graph.edge_count());
            incidence[edge.id]++;
        }
    }
    for (int count : incidence) {
        if (count != 0) assert(count == expected);
    }
#else
    (void)graph;
    (void)expected;
#endif
}

template <class Field, class Weight>
Field count_arborescences(
    const Graph<Weight>& graph,
    int root,
    bool outward
) {
    const int n = graph.size();
    assert(0 <= root && root < n);
    assert_edge_incidence(graph, 1);

    matrix::Matrix<Field> minor(n - 1, n - 1);
    for (int vertex = 0; vertex < n; vertex++) {
        for (const Edge<Weight>& edge : graph[vertex]) {
            if (!edge.alive || edge.from == edge.to) continue;
            const int row = outward ? edge.to : edge.from;
            const int col = outward ? edge.from : edge.to;
            if (row == root) continue;

            const Field weight(edge.cost);
            const int reduced_row = minor_index(row, root);
            minor[reduced_row][reduced_row] += weight;
            if (col != root) {
                minor[reduced_row][minor_index(col, root)] -= weight;
            }
        }
    }
    return matrix::determinant(std::move(minor));
}

}  // namespace matrix_tree_detail

// Returns the total weight of all undirected spanning trees. The weight of a
// tree is the product of its edge costs.
template <class Field, class Weight>
Field count_spanning_trees(const Graph<Weight>& graph) {
    const int n = graph.size();
    assert(n > 0);
    matrix_tree_detail::assert_edge_incidence(graph, 2);

    const int removed = n - 1;
    matrix::Matrix<Field> minor(n - 1, n - 1);
    for (int vertex = 0; vertex < n; vertex++) {
        for (const Edge<Weight>& edge : graph[vertex]) {
            if (!edge.alive || edge.from >= edge.to) continue;
            const int from = edge.from;
            const int to = edge.to;
            const Field weight(edge.cost);

            if (from != removed) {
                const int reduced_from = matrix_tree_detail::minor_index(from, removed);
                minor[reduced_from][reduced_from] += weight;
            }
            if (to != removed) {
                const int reduced_to = matrix_tree_detail::minor_index(to, removed);
                minor[reduced_to][reduced_to] += weight;
            }
            if (from != removed && to != removed) {
                const int reduced_from = matrix_tree_detail::minor_index(from, removed);
                const int reduced_to = matrix_tree_detail::minor_index(to, removed);
                minor[reduced_from][reduced_to] -= weight;
                minor[reduced_to][reduced_from] -= weight;
            }
        }
    }
    return matrix::determinant(std::move(minor));
}

// Counts directed spanning trees whose edges point away from root, so every
// vertex is reachable from root.
template <class Field, class Weight>
Field count_out_arborescences(const Graph<Weight>& graph, int root) {
    return matrix_tree_detail::count_arborescences<Field>(graph, root, true);
}

// Counts directed spanning trees whose edges point toward root, so root is
// reachable from every vertex.
template <class Field, class Weight>
Field count_in_arborescences(const Graph<Weight>& graph, int root) {
    return matrix_tree_detail::count_arborescences<Field>(graph, root, false);
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/minimum_steiner_tree.hpp"



#line 15 "graph/minimum_steiner_tree.hpp"

#line 17 "graph/minimum_steiner_tree.hpp"

namespace m1une {
namespace graph {

template <class Cost>
struct MinimumSteinerTreeResult {
    Cost cost;
    std::vector<int> edge_ids;
    std::vector<int> vertices;
};

namespace internal {

inline std::vector<int> steiner_terminals(int n, std::vector<int> terminals) {
    for (int v : terminals) assert(0 <= v && v < n);
    std::sort(terminals.begin(), terminals.end());
    terminals.erase(std::unique(terminals.begin(), terminals.end()), terminals.end());
    assert(terminals.size() < std::numeric_limits<std::size_t>::digits);
    return terminals;
}

template <class Cost>
struct MinimumSteinerTreeDp {
    Cost cost;
    Cost inf;
    std::size_t states;
    std::size_t width;
    std::vector<Cost> dp;
    std::vector<int> terminals;
};

template <class Cost, class GraphCost, class EdgeCost>
std::optional<MinimumSteinerTreeDp<Cost>> minimum_steiner_tree_dp(
    const Graph<GraphCost>& g,
    std::vector<int> terminals,
    const std::vector<Cost>& vertex_cost,
    EdgeCost edge_cost,
    Cost inf
) {
    const int n = g.size();
    assert(vertex_cost.size() == std::size_t(n));
    for (Cost cost : vertex_cost) assert(Cost(0) <= cost);
    terminals = steiner_terminals(n, std::move(terminals));
    const int k = int(terminals.size());
    if (k == 0) return MinimumSteinerTreeDp<Cost>{Cost(0), inf, 1, std::size_t(n), {}, {}};

    assert(Cost(0) < inf);
    for (int v = 0; v < n; v++) {
        for (const auto& edge : g[v]) {
            if (edge.alive) assert(Cost(0) <= edge_cost(edge));
        }
    }

    const std::size_t states = std::size_t(1) << k;
    const std::size_t width = std::size_t(n);
    assert(width <= std::numeric_limits<std::size_t>::max() / states);
    std::vector<Cost> dp(states * width, inf);
    for (int i = 0; i < k; i++) {
        const int terminal = terminals[i];
        if (vertex_cost[terminal] < inf) {
            dp[(std::size_t(1) << i) * width + std::size_t(terminal)] = vertex_cost[terminal];
        }
    }

    using QueueEntry = std::pair<Cost, int>;
    for (std::size_t mask = 1; mask < states; mask++) {
        const std::size_t mask_offset = mask * width;
        for (std::size_t sub = (mask - 1) & mask; sub != 0; sub = (sub - 1) & mask) {
            const std::size_t other = mask ^ sub;
            if (sub > other) continue;
            const std::size_t sub_offset = sub * width;
            const std::size_t other_offset = other * width;
            for (int v = 0; v < n; v++) {
                const std::size_t vertex = std::size_t(v);
                const Cost left = dp[sub_offset + vertex];
                const Cost right = dp[other_offset + vertex];
                if (left == inf || right == inf) continue;
                assert(vertex_cost[v] <= right);
                const Cost extra = right - vertex_cost[v];
                if (left > inf - extra) continue;
                const Cost candidate = left + extra;
                Cost& current = dp[mask_offset + vertex];
                if (candidate < current) current = candidate;
            }
        }

        std::priority_queue<QueueEntry, std::vector<QueueEntry>, std::greater<QueueEntry>> queue;
        for (int v = 0; v < n; v++) {
            const Cost distance = dp[mask_offset + std::size_t(v)];
            if (distance != inf) queue.emplace(distance, v);
        }
        while (!queue.empty()) {
            auto [distance, v] = queue.top();
            queue.pop();
            if (distance != dp[mask_offset + std::size_t(v)]) continue;
            for (const auto& edge : g[v]) {
                if (!edge.alive) continue;
                const Cost cost = edge_cost(edge);
                if (cost >= inf || vertex_cost[edge.to] > inf - cost) continue;
                const Cost extra = cost + vertex_cost[edge.to];
                if (distance > inf - extra) continue;
                const Cost candidate = distance + extra;
                Cost& current = dp[mask_offset + std::size_t(edge.to)];
                if (current <= candidate) continue;
                current = candidate;
                queue.emplace(candidate, edge.to);
            }
        }
    }

    const auto answer_begin = dp.begin() + (states - 1) * width;
    const Cost answer = *std::min_element(answer_begin, dp.end());
    if (answer == inf) return std::nullopt;
    return MinimumSteinerTreeDp<Cost>{
        answer,
        inf,
        states,
        width,
        std::move(dp),
        std::move(terminals)
    };
}

template <class T>
std::optional<MinimumSteinerTreeDp<int>> minimum_steiner_tree_unweighted_dp(
    const Graph<T>& g,
    std::vector<int> terminals
) {
    const int n = g.size();
    terminals = steiner_terminals(n, std::move(terminals));
    const int k = int(terminals.size());
    if (k == 0) return MinimumSteinerTreeDp<int>{0, n, 1, std::size_t(n), {}, {}};

    const std::size_t states = std::size_t(1) << k;
    const std::size_t width = std::size_t(n);
    assert(width <= std::numeric_limits<std::size_t>::max() / states);
    const int inf = n;
    std::vector<int> dp(states * width, inf);
    for (int i = 0; i < k; i++) {
        dp[(std::size_t(1) << i) * width + std::size_t(terminals[i])] = 0;
    }

    for (std::size_t mask = 1; mask < states; mask++) {
        const std::size_t mask_offset = mask * width;
        for (std::size_t sub = (mask - 1) & mask; sub != 0; sub = (sub - 1) & mask) {
            const std::size_t other = mask ^ sub;
            if (sub > other) continue;
            const std::size_t sub_offset = sub * width;
            const std::size_t other_offset = other * width;
            for (int v = 0; v < n; v++) {
                const std::size_t vertex = std::size_t(v);
                const int candidate = dp[sub_offset + vertex] + dp[other_offset + vertex];
                int& current = dp[mask_offset + vertex];
                if (candidate < current) current = candidate;
            }
        }

        std::vector<int> bucket_head(n, -1);
        std::vector<int> entry_vertex;
        std::vector<int> entry_next;
        entry_vertex.reserve(2 * width);
        entry_next.reserve(2 * width);
        auto push = [&](int distance, int v) {
            entry_vertex.push_back(v);
            entry_next.push_back(bucket_head[distance]);
            bucket_head[distance] = int(entry_vertex.size()) - 1;
        };
        for (int v = 0; v < n; v++) {
            const int distance = dp[mask_offset + std::size_t(v)];
            if (distance != inf) push(distance, v);
        }
        for (int distance = 0; distance < n; distance++) {
            for (int entry = bucket_head[distance]; entry != -1; entry = entry_next[entry]) {
                const int v = entry_vertex[entry];
                if (dp[mask_offset + std::size_t(v)] != distance) continue;
                for (const auto& edge : g[v]) {
                    if (!edge.alive) continue;
                    int& current = dp[mask_offset + std::size_t(edge.to)];
                    if (distance + 1 >= current) continue;
                    current = distance + 1;
                    push(current, edge.to);
                }
            }
        }
    }

    const auto answer_begin = dp.begin() + (states - 1) * width;
    const int answer = *std::min_element(answer_begin, dp.end());
    if (answer == inf) return std::nullopt;
    return MinimumSteinerTreeDp<int>{
        answer,
        inf,
        states,
        width,
        std::move(dp),
        std::move(terminals)
    };
}

template <class Cost, class GraphCost, class EdgeCost>
MinimumSteinerTreeResult<Cost> restore_minimum_steiner_tree(
    const Graph<GraphCost>& g,
    const MinimumSteinerTreeDp<Cost>& data,
    const std::vector<Cost>& vertex_cost,
    EdgeCost edge_cost
) {
    MinimumSteinerTreeResult<Cost> result;
    result.cost = data.cost;
    if (data.terminals.empty()) return result;

    const int n = g.size();
    const std::size_t cells = data.states * data.width;
    std::vector<char> state(cells, 0);
    std::vector<char> selected_edge(g.edge_count(), false);

    std::function<bool(std::size_t, int)> restore = [&](std::size_t mask, int start) {
        const std::size_t position = mask * data.width + std::size_t(start);
        if (state[position] == 2) return true;
        if (state[position] == 1) return false;
        state[position] = 1;

        std::vector<int> search_parent(n, -2), search_edge(n, -1), stack;
        search_parent[start] = -1;
        stack.push_back(start);
        int seed = -1;
        std::size_t seed_split = 0;

        while (!stack.empty() && seed == -1) {
            const int v = stack.back();
            stack.pop_back();
            const std::size_t vertex_position = mask * data.width + std::size_t(v);
            const Cost current = data.dp[vertex_position];

            if (v != start && state[vertex_position] == 2) {
                seed = v;
                break;
            }
            if ((mask & (mask - 1)) == 0) {
                const int terminal_index = int(std::countr_zero(mask));
                if (v == data.terminals[terminal_index] && current == vertex_cost[v]) {
                    seed = v;
                    break;
                }
            }
            for (std::size_t sub = (mask - 1) & mask; sub != 0; sub = (sub - 1) & mask) {
                const std::size_t other = mask ^ sub;
                if (sub > other) continue;
                const Cost left = data.dp[sub * data.width + std::size_t(v)];
                const Cost right = data.dp[other * data.width + std::size_t(v)];
                if (left == data.inf || right == data.inf || right < vertex_cost[v]) continue;
                const Cost extra = right - vertex_cost[v];
                if (left > data.inf - extra || left + extra != current) continue;
                seed = v;
                seed_split = sub;
                break;
            }
            if (seed != -1) break;

            for (const auto& edge : g[v]) {
                if (!edge.alive || search_parent[edge.to] != -2) continue;
                const Cost cost = edge_cost(edge);
                if (cost >= data.inf || vertex_cost[v] > data.inf - cost) continue;
                const Cost extra = cost + vertex_cost[v];
                const Cost previous = data.dp[mask * data.width + std::size_t(edge.to)];
                if (previous == data.inf || previous > data.inf - extra) continue;
                if (previous + extra != current) continue;
                search_parent[edge.to] = v;
                search_edge[edge.to] = edge.id;
                stack.push_back(edge.to);
            }
        }

        if (seed == -1) {
            state[position] = 0;
            return false;
        }
        if (seed_split != 0) {
            const bool restored_left = restore(seed_split, seed);
            const bool restored_right = restore(mask ^ seed_split, seed);
            assert(restored_left && restored_right);
            if (!restored_left || !restored_right) {
                state[position] = 0;
                return false;
            }
        }

        for (int v = seed; v != -1; v = search_parent[v]) {
            state[mask * data.width + std::size_t(v)] = 2;
            if (search_parent[v] == -1) continue;
            const int id = search_edge[v];
            assert(0 <= id && id < g.edge_count());
            selected_edge[id] = true;
        }
        return true;
    };

    int root = -1;
    const std::size_t full_mask = data.states - 1;
    for (int v = 0; v < n; v++) {
        if (data.dp[full_mask * data.width + std::size_t(v)] == data.cost) {
            root = v;
            break;
        }
    }
    assert(root != -1);
    const bool restored = restore(full_mask, root);
    assert(restored);
    (void)restored;

    std::vector<Edge<GraphCost>> edge_by_id(g.edge_count());
    std::vector<char> has_edge(g.edge_count(), false);
    for (const auto& edge : g.edges()) {
        edge_by_id[edge.id] = edge;
        has_edge[edge.id] = true;
    }

    std::vector<int> parent(n), component_size(n, 1);
    for (int v = 0; v < n; v++) parent[v] = v;
    auto leader = [&](auto&& self, int v) -> int {
        if (parent[v] == v) return v;
        return parent[v] = self(self, parent[v]);
    };

    std::vector<char> tree_edge(g.edge_count(), false);
    for (int id = 0; id < g.edge_count(); id++) {
        if (!selected_edge[id]) continue;
        assert(has_edge[id]);
        const auto& edge = edge_by_id[id];
        int u = leader(leader, edge.from);
        int v = leader(leader, edge.to);
        if (u == v) continue;
        if (component_size[u] < component_size[v]) std::swap(u, v);
        parent[v] = u;
        component_size[u] += component_size[v];
        tree_edge[id] = true;
    }

    std::vector<std::vector<std::pair<int, int>>> tree(n);
    std::vector<int> degree(n, 0);
    std::vector<char> in_tree(n, false), is_terminal(n, false);
    for (int terminal : data.terminals) {
        in_tree[terminal] = true;
        is_terminal[terminal] = true;
    }
    for (int id = 0; id < g.edge_count(); id++) {
        if (!tree_edge[id]) continue;
        const auto& edge = edge_by_id[id];
        tree[edge.from].emplace_back(edge.to, id);
        tree[edge.to].emplace_back(edge.from, id);
        degree[edge.from]++;
        degree[edge.to]++;
        in_tree[edge.from] = true;
        in_tree[edge.to] = true;
    }

    std::queue<int> leaves;
    for (int v = 0; v < n; v++) {
        if (in_tree[v] && !is_terminal[v] && degree[v] <= 1) leaves.push(v);
    }
    std::vector<char> removed_vertex(n, false), removed_edge(g.edge_count(), false);
    while (!leaves.empty()) {
        const int v = leaves.front();
        leaves.pop();
        if (removed_vertex[v] || is_terminal[v] || degree[v] > 1) continue;
        removed_vertex[v] = true;
        for (auto [to, id] : tree[v]) {
            if (removed_edge[id]) continue;
            removed_edge[id] = true;
            degree[v]--;
            degree[to]--;
            if (!is_terminal[to] && degree[to] <= 1) leaves.push(to);
            break;
        }
    }

    Cost restored_cost = Cost(0);
    for (int id = 0; id < g.edge_count(); id++) {
        if (!tree_edge[id] || removed_edge[id]) continue;
        result.edge_ids.push_back(id);
        restored_cost += edge_cost(edge_by_id[id]);
    }
    for (int v = 0; v < n; v++) {
        if (!in_tree[v] || removed_vertex[v]) continue;
        result.vertices.push_back(v);
        restored_cost += vertex_cost[v];
    }
    if constexpr (std::is_integral_v<Cost>) assert(restored_cost == result.cost);
    result.cost = restored_cost;
    return result;
}

}  // namespace internal

template <class T>
std::optional<T> minimum_steiner_tree(
    const Graph<T>& g,
    std::vector<int> terminals,
    const std::vector<T>& vertex_cost,
    T inf = std::numeric_limits<T>::max() / T(4)
) {
    auto result = internal::minimum_steiner_tree_dp(
        g,
        std::move(terminals),
        vertex_cost,
        [](const Edge<T>& edge) { return edge.cost; },
        inf
    );
    if (!result) return std::nullopt;
    return result->cost;
}

template <class T>
std::optional<T> minimum_steiner_tree(
    const Graph<T>& g,
    std::vector<int> terminals,
    T inf = std::numeric_limits<T>::max() / T(4)
) {
    return minimum_steiner_tree(g, std::move(terminals), std::vector<T>(g.size(), T(0)), inf);
}

template <class GraphCost, class Cost>
std::optional<Cost> minimum_steiner_tree_unweighted(
    const Graph<GraphCost>& g,
    std::vector<int> terminals,
    const std::vector<Cost>& vertex_cost,
    Cost inf = std::numeric_limits<Cost>::max() / Cost(4)
) {
    auto result = internal::minimum_steiner_tree_dp(
        g,
        std::move(terminals),
        vertex_cost,
        [](const Edge<GraphCost>&) { return Cost(1); },
        inf
    );
    if (!result) return std::nullopt;
    return result->cost;
}

template <class T>
std::optional<MinimumSteinerTreeResult<T>> build_minimum_steiner_tree(
    const Graph<T>& g,
    std::vector<int> terminals,
    const std::vector<T>& vertex_cost,
    T inf = std::numeric_limits<T>::max() / T(4)
) {
    auto data = internal::minimum_steiner_tree_dp(
        g,
        std::move(terminals),
        vertex_cost,
        [](const Edge<T>& edge) { return edge.cost; },
        inf
    );
    if (!data) return std::nullopt;
    return internal::restore_minimum_steiner_tree(
        g,
        *data,
        vertex_cost,
        [](const Edge<T>& edge) { return edge.cost; }
    );
}

template <class T>
std::optional<MinimumSteinerTreeResult<T>> build_minimum_steiner_tree(
    const Graph<T>& g,
    std::vector<int> terminals,
    T inf = std::numeric_limits<T>::max() / T(4)
) {
    std::vector<T> vertex_cost(g.size(), T(0));
    return build_minimum_steiner_tree(g, std::move(terminals), vertex_cost, inf);
}

template <class GraphCost, class Cost>
std::optional<MinimumSteinerTreeResult<Cost>> build_minimum_steiner_tree_unweighted(
    const Graph<GraphCost>& g,
    std::vector<int> terminals,
    const std::vector<Cost>& vertex_cost,
    Cost inf = std::numeric_limits<Cost>::max() / Cost(4)
) {
    auto data = internal::minimum_steiner_tree_dp(
        g,
        std::move(terminals),
        vertex_cost,
        [](const Edge<GraphCost>&) { return Cost(1); },
        inf
    );
    if (!data) return std::nullopt;
    return internal::restore_minimum_steiner_tree(
        g,
        *data,
        vertex_cost,
        [](const Edge<GraphCost>&) { return Cost(1); }
    );
}

template <class T>
std::optional<MinimumSteinerTreeResult<int>> build_minimum_steiner_tree_unweighted(
    const Graph<T>& g,
    std::vector<int> terminals
) {
    auto data = internal::minimum_steiner_tree_unweighted_dp(g, std::move(terminals));
    if (!data) return std::nullopt;
    std::vector<int> vertex_cost(g.size(), 0);
    return internal::restore_minimum_steiner_tree(
        g,
        *data,
        vertex_cost,
        [](const Edge<T>&) { return 1; }
    );
}

template <class T>
std::optional<int> minimum_steiner_tree_unweighted(
    const Graph<T>& g,
    std::vector<int> terminals
) {
    auto result = internal::minimum_steiner_tree_unweighted_dp(g, std::move(terminals));
    if (!result) return std::nullopt;
    return result->cost;
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/namori.hpp"



#line 9 "graph/namori.hpp"

#line 11 "graph/namori.hpp"

namespace m1une {
namespace graph {

template <class T>
struct NamoriDecomposition {
    int component_count;
    std::vector<std::vector<int>> cycles;
    std::vector<std::vector<int>> cycle_edge_ids;
    std::vector<std::vector<T>> cycle_edge_costs;

    std::vector<bool> on_cycle;
    std::vector<int> component;
    std::vector<int> cycle_root;
    std::vector<int> cycle_position;
    std::vector<int> parent;
    std::vector<int> parent_edge;
    std::vector<int> depth;
    std::vector<T> dist_to_cycle;
    std::vector<std::vector<int>> children;

    bool same_component(int u, int v) const {
        assert(0 <= u && u < int(component.size()));
        assert(0 <= v && v < int(component.size()));
        return component[u] == component[v];
    }

    bool same_tree(int u, int v) const {
        assert(0 <= u && u < int(cycle_root.size()));
        assert(0 <= v && v < int(cycle_root.size()));
        return cycle_root[u] == cycle_root[v];
    }
};

template <class T>
std::optional<NamoriDecomposition<T>> namori_decomposition(const Graph<T>& graph) {
    int n = graph.size();
    NamoriDecomposition<T> result;
    result.component_count = 0;
    result.on_cycle.assign(n, false);
    result.component.assign(n, -1);
    result.cycle_root.assign(n, -1);
    result.cycle_position.assign(n, -1);
    result.parent.assign(n, -1);
    result.parent_edge.assign(n, -1);
    result.depth.assign(n, 0);
    result.dist_to_cycle.assign(n, T(0));
    result.children.assign(n, {});
    if (n == 0) return result;

    std::vector<int> degree(n, 0);
    for (int v = 0; v < n; v++) {
        for (const auto& edge : graph[v]) {
            if (edge.alive) degree[v]++;
        }
    }

    std::queue<int> queue;
    std::vector<bool> removed(n, false);
    for (int v = 0; v < n; v++) {
        if (degree[v] <= 1) queue.push(v);
    }
    while (!queue.empty()) {
        int v = queue.front();
        queue.pop();
        if (removed[v] || degree[v] > 1) continue;
        removed[v] = true;
        for (const auto& edge : graph[v]) {
            if (!edge.alive || removed[edge.to]) continue;
            degree[edge.to]--;
            if (degree[edge.to] == 1) queue.push(edge.to);
        }
    }

    for (int v = 0; v < n; v++) {
        result.on_cycle[v] = !removed[v];
    }
    for (int v = 0; v < n; v++) {
        if (!result.on_cycle[v]) continue;
        int cycle_degree = 0;
        for (const auto& edge : graph[v]) {
            if (edge.alive && result.on_cycle[edge.to]) cycle_degree++;
        }
        if (cycle_degree != 2) return std::nullopt;
    }

    std::vector<bool> cycle_visited(n, false);
    for (int start = 0; start < n; start++) {
        if (!result.on_cycle[start] || cycle_visited[start]) continue;
        int component_id = int(result.cycles.size());
        std::vector<int> vertices;
        std::vector<int> edge_ids;
        std::vector<T> edge_costs;

        int current = start;
        int previous_edge = -1;
        while (true) {
            if (cycle_visited[current]) return std::nullopt;
            cycle_visited[current] = true;
            vertices.push_back(current);

            int next_vertex = -1;
            int next_edge = -1;
            T next_cost = T(0);
            for (const auto& edge : graph[current]) {
                if (!edge.alive || !result.on_cycle[edge.to] || edge.id == previous_edge) continue;
                next_vertex = edge.to;
                next_edge = edge.id;
                next_cost = edge.cost;
                break;
            }
            if (next_edge == -1) return std::nullopt;
            edge_ids.push_back(next_edge);
            edge_costs.push_back(next_cost);
            if (next_vertex == start) break;
            previous_edge = next_edge;
            current = next_vertex;
            if (int(vertices.size()) > n) return std::nullopt;
        }

        for (int position = 0; position < int(vertices.size()); position++) {
            int v = vertices[position];
            result.component[v] = component_id;
            result.cycle_root[v] = v;
            result.cycle_position[v] = position;
        }
        result.cycles.push_back(std::move(vertices));
        result.cycle_edge_ids.push_back(std::move(edge_ids));
        result.cycle_edge_costs.push_back(std::move(edge_costs));
    }
    if (result.cycles.empty()) return std::nullopt;

    std::vector<int> stack;
    stack.reserve(n);
    for (const auto& cycle : result.cycles) {
        for (int v : cycle) stack.push_back(v);
    }
    while (!stack.empty()) {
        int v = stack.back();
        stack.pop_back();
        for (const auto& edge : graph[v]) {
            if (!edge.alive || result.on_cycle[edge.to] || edge.id == result.parent_edge[v]) continue;
            int to = edge.to;
            if (result.component[to] != -1) continue;
            result.component[to] = result.component[v];
            result.cycle_root[to] = result.cycle_root[v];
            result.cycle_position[to] = result.cycle_position[v];
            result.parent[to] = v;
            result.parent_edge[to] = edge.id;
            result.depth[to] = result.depth[v] + 1;
            result.dist_to_cycle[to] = result.dist_to_cycle[v] + edge.cost;
            result.children[v].push_back(to);
            stack.push_back(to);
        }
    }
    for (int v = 0; v < n; v++) {
        if (result.component[v] == -1) return std::nullopt;
    }

    result.component_count = int(result.cycles.size());
    return result;
}

template <class T>
std::optional<NamoriDecomposition<T>> decompose_namori(const Graph<T>& graph) {
    return namori_decomposition(graph);
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/replacement_paths.hpp"



#line 11 "graph/replacement_paths.hpp"

#line 1 "graph/dijkstra.hpp"



#line 8 "graph/dijkstra.hpp"

#line 10 "graph/dijkstra.hpp"

namespace m1une {
namespace graph {

template <class T>
struct DijkstraResult {
    std::vector<T> dist;
    std::vector<char> reached;
    std::vector<int> parent;
    std::vector<int> parent_edge;
    T inf = T();

    bool reachable(int v) const {
        assert(0 <= v && v < int(dist.size()));
        return reached[v];
    }

    std::vector<int> path(int t) const {
        assert(reachable(t));
        std::vector<int> result;
        for (int v = t; v != -1; v = parent[v]) result.push_back(v);
        std::reverse(result.begin(), result.end());
        return result;
    }
};

namespace internal {

template <class T>
class DijkstraHeap {
   private:
    const std::vector<T>& dist_;
    std::vector<int> heap_;
    std::vector<int> position_;

    bool less(int first, int second) const {
        return dist_[heap_[first]] < dist_[heap_[second]];
    }

    void swap_nodes(int first, int second) {
        std::swap(heap_[first], heap_[second]);
        position_[heap_[first]] = first;
        position_[heap_[second]] = second;
    }

    void sift_up(int index) {
        while (index != 0) {
            const int parent = (index - 1) / 2;
            if (!less(index, parent)) break;
            swap_nodes(index, parent);
            index = parent;
        }
    }

    void sift_down(int index) {
        while (2 * index + 1 < int(heap_.size())) {
            int child = 2 * index + 1;
            if (child + 1 < int(heap_.size()) && less(child + 1, child)) {
                ++child;
            }
            if (!less(child, index)) break;
            swap_nodes(index, child);
            index = child;
        }
    }

   public:
    DijkstraHeap(const std::vector<T>& dist, int size)
        : dist_(dist), position_(size, -1) {
        heap_.reserve(size);
    }

    bool empty() const {
        return heap_.empty();
    }

    void push_or_decrease(int vertex) {
        int& position = position_[vertex];
        if (position == -1) {
            position = int(heap_.size());
            heap_.push_back(vertex);
        }
        sift_up(position);
    }

    int pop_min() {
        const int result = heap_.front();
        position_[result] = -1;
        if (heap_.size() == 1) {
            heap_.pop_back();
            return result;
        }
        heap_.front() = heap_.back();
        position_[heap_.front()] = 0;
        heap_.pop_back();
        sift_down(0);
        return result;
    }
};

}  // namespace internal

template <class T>
DijkstraResult<T> dijkstra(const Graph<T>& g,
                           const std::vector<int>& sources) {
    int n = g.size();
    DijkstraResult<T> result;
    result.dist.resize(n);
    result.reached.assign(n, false);
    result.parent.assign(n, -1);
    result.parent_edge.assign(n, -1);

    internal::DijkstraHeap<T> que(result.dist, n);
    for (int s : sources) {
        assert(0 <= s && s < n);
        if (result.reached[s]) continue;
        result.reached[s] = true;
        result.dist[s] = T();
        que.push_or_decrease(s);
    }

    while (!que.empty()) {
        const int current = que.pop_min();
        for (const auto& e : g[current]) {
            if (!e.alive) continue;
            T nd = result.dist[current] + e.cost;
            if (result.reached[e.to] && !(nd < result.dist[e.to])) continue;
            result.reached[e.to] = true;
            result.dist[e.to] = std::move(nd);
            result.parent[e.to] = current;
            result.parent_edge[e.to] = e.id;
            que.push_or_decrease(e.to);
        }
    }

    return result;
}

template <class T>
DijkstraResult<T> dijkstra(const Graph<T>& g, int s) {
    return dijkstra(g, std::vector<int>{s});
}

// Compatibility overload: unreachable distances are replaced by inf after the
// search. Reachability itself never depends on this sentinel.
template <class T>
DijkstraResult<T> dijkstra(const Graph<T>& g,
                           const std::vector<int>& sources, const T& inf) {
    DijkstraResult<T> result = dijkstra(g, sources);
    result.inf = inf;
    for (int v = 0; v < int(result.dist.size()); v++) {
        if (!result.reachable(v)) result.dist[v] = inf;
    }
    return result;
}

template <class T>
DijkstraResult<T> dijkstra(const Graph<T>& g, int s, const T& inf) {
    return dijkstra(g, std::vector<int>{s}, inf);
}

}  // namespace graph
}  // namespace m1une


#line 14 "graph/replacement_paths.hpp"

namespace m1une {
namespace graph {

struct GraphPath {
    std::vector<int> vertices;
    std::vector<int> edges;
};

template <class T>
struct EdgeReplacementPathsResult {
    GraphPath path;
    std::vector<T> replacement_dist;
    T inf;

    bool reachable(int path_edge_index) const {
        assert(0 <= path_edge_index && path_edge_index < int(replacement_dist.size()));
        return replacement_dist[path_edge_index] != inf;
    }
};

template <class T>
struct VertexReplacementPathsResult {
    GraphPath path;
    std::vector<T> replacement_dist;
    T inf;

    bool reachable(int path_vertex_index) const {
        assert(0 <= path_vertex_index && path_vertex_index < int(replacement_dist.size()));
        return replacement_dist[path_vertex_index] != inf;
    }
};

namespace internal {

template <class T>
T replacement_paths_safe_add(T a, T b, T inf) {
    if (a >= inf || b >= inf) return inf;
    if (a > inf - b) return inf;
    return a + b;
}

template <class T>
DijkstraResult<T> replacement_paths_dijkstra(const Graph<T>& g, int s, T inf) {
    int n = g.size();
    assert(0 <= s && s < n);
    DijkstraResult<T> result;
    result.dist.assign(n, inf);
    result.reached.assign(n, false);
    result.parent.assign(n, -1);
    result.parent_edge.assign(n, -1);
    result.inf = inf;

    using P = std::pair<T, int>;
    std::priority_queue<P, std::vector<P>, std::greater<P>> que;
    result.dist[s] = T(0);
    result.reached[s] = true;
    que.emplace(T(0), s);
    while (!que.empty()) {
        auto [d, v] = que.top();
        que.pop();
        if (result.dist[v] != d) continue;
        for (const auto& e : g[v]) {
            if (!e.alive) continue;
            T nd = replacement_paths_safe_add(d, e.cost, inf);
            if (result.dist[e.to] <= nd) continue;
            result.reached[e.to] = true;
            result.dist[e.to] = nd;
            result.parent[e.to] = v;
            result.parent_edge[e.to] = e.id;
            que.emplace(nd, e.to);
        }
    }
    return result;
}

template <class T>
std::vector<Edge<T>> replacement_paths_validate_graph(const Graph<T>& g, T inf) {
    assert(T(0) < inf);
    std::vector<int> occurrence(g.edge_count(), 0);
    std::vector<Edge<T>> edge_by_id(g.edge_count());
    for (int v = 0; v < g.size(); v++) {
        for (const auto& e : g[v]) {
            assert(e.from == v);
            assert(0 <= e.to && e.to < g.size());
            assert(0 <= e.id && e.id < g.edge_count());
            if (e.alive) assert(T(0) < e.cost);
            if (occurrence[e.id] == 0) {
                edge_by_id[e.id] = e;
            } else {
                assert(occurrence[e.id] == 1);
                const auto& other = edge_by_id[e.id];
                assert(e.from == other.to && e.to == other.from);
                assert(e.cost == other.cost && e.alive == other.alive);
            }
            occurrence[e.id]++;
        }
    }

    for (int id = 0; id < g.edge_count(); id++) {
        // add_edge creates exactly two mutually reversed arcs with one logical id.
        assert(occurrence[id] == 2);
    }
    return edge_by_id;
}

template <class T>
void replacement_paths_validate_path(const Graph<T>& g, const GraphPath& path,
                                     const std::vector<Edge<T>>& edge_by_id,
                                     const DijkstraResult<T>& from_s, T inf) {
    assert(!path.vertices.empty());
    assert(path.edges.size() + 1 == path.vertices.size());
    std::vector<char> used_vertex(g.size(), false);
    for (int v : path.vertices) {
        assert(0 <= v && v < g.size());
        assert(!used_vertex[v]);
        used_vertex[v] = true;
    }

    T path_cost = T(0);
    for (int i = 0; i < int(path.edges.size()); i++) {
        int id = path.edges[i];
        assert(0 <= id && id < g.edge_count());
        assert(g.is_edge_alive(id));
        const auto& e = edge_by_id[id];
        int u = path.vertices[i];
        int v = path.vertices[i + 1];
        assert((e.from == u && e.to == v) || (e.from == v && e.to == u));
        assert(T(0) < e.cost);
        path_cost = replacement_paths_safe_add(path_cost, e.cost, inf);
    }
    assert(from_s.reachable(path.vertices.back()));
    assert(path_cost == from_s.dist[path.vertices.back()]);
}

template <class T>
GraphPath replacement_paths_restore_path(const DijkstraResult<T>& result, int s, int t) {
    assert(result.reachable(t));
    GraphPath path;
    for (int v = t; v != s; v = result.parent[v]) {
        assert(v != -1 && result.parent[v] != -1 && result.parent_edge[v] != -1);
        path.vertices.push_back(v);
        path.edges.push_back(result.parent_edge[v]);
    }
    path.vertices.push_back(s);
    std::reverse(path.vertices.begin(), path.vertices.end());
    std::reverse(path.edges.begin(), path.edges.end());
    return path;
}

template <class T>
struct ReplacementPathsData {
    GraphPath path;
    std::vector<T> dist_s;
    std::vector<T> dist_t;
    std::vector<int> block;
    std::vector<char> is_path_edge;
    std::vector<Edge<T>> edge_by_id;
    T inf;
};

template <class T>
ReplacementPathsData<T> replacement_paths_prepare(const Graph<T>& g, const GraphPath& path,
                                                   T inf, const DijkstraResult<T>* known_from_s) {
    auto edge_by_id = replacement_paths_validate_graph(g, inf);
    int s = path.vertices.front();
    int t = path.vertices.back();
    auto computed_from_s = known_from_s == nullptr
                               ? replacement_paths_dijkstra(g, s, inf)
                               : DijkstraResult<T>();
    const auto& from_s = known_from_s == nullptr ? computed_from_s : *known_from_s;
    replacement_paths_validate_path(g, path, edge_by_id, from_s, inf);
    auto from_t = replacement_paths_dijkstra(g, t, inf);

    int n = g.size();
    std::vector<int> path_position(n, -1);
    std::vector<char> is_path_edge(g.edge_count(), false);
    for (int i = 0; i < int(path.vertices.size()); i++) path_position[path.vertices[i]] = i;
    for (int id : path.edges) is_path_edge[id] = true;

    std::vector<int> parent(n, -1);
    for (int i = 0; i < int(path.edges.size()); i++) {
        int v = path.vertices[i + 1];
        parent[v] = path.vertices[i];
        const auto& e = edge_by_id[path.edges[i]];
        assert(replacement_paths_safe_add(from_s.dist[parent[v]], e.cost, inf) == from_s.dist[v]);
    }
    for (int v = 0; v < n; v++) {
        if (!from_s.reachable(v) || v == s || path_position[v] != -1) continue;
        for (const auto& e : g[v]) {
            if (!e.alive || !from_s.reachable(e.to)) continue;
            if (replacement_paths_safe_add(from_s.dist[e.to], e.cost, inf) != from_s.dist[v]) {
                continue;
            }
            parent[v] = e.to;
            break;
        }
        assert(parent[v] != -1);
        assert(from_s.dist[parent[v]] < from_s.dist[v]);
    }

    std::vector<std::vector<int>> children(n);
    for (int v = 0; v < n; v++) {
        if (parent[v] != -1) children[parent[v]].push_back(v);
    }
    std::vector<int> block(n, -1);
    block[s] = 0;
    std::vector<int> stack = {s};
    while (!stack.empty()) {
        int v = stack.back();
        stack.pop_back();
        for (int to : children[v]) {
            block[to] = path_position[to] == -1 ? block[v] : path_position[to];
            stack.push_back(to);
        }
    }
    for (int v = 0; v < n; v++) assert(!from_s.reachable(v) || block[v] != -1);

    return {path, from_s.dist, from_t.dist, block, is_path_edge, edge_by_id, inf};
}

template <class T>
class ReplacementPathsRangeChmin {
   private:
    int _size;
    std::vector<T> _lazy;

   public:
    ReplacementPathsRangeChmin(int n, T inf) : _size(1) {
        while (_size < n) _size <<= 1;
        _lazy.assign(2 * _size, inf);
    }

    void apply(int l, int r, T value) {
        assert(0 <= l && l <= r && r <= _size);
        for (l += _size, r += _size; l < r; l >>= 1, r >>= 1) {
            if (l & 1) {
                _lazy[l] = std::min(_lazy[l], value);
                l++;
            }
            if (r & 1) {
                --r;
                _lazy[r] = std::min(_lazy[r], value);
            }
        }
    }

    std::vector<T> values(int n) {
        for (int v = 1; v < _size; v++) {
            _lazy[2 * v] = std::min(_lazy[2 * v], _lazy[v]);
            _lazy[2 * v + 1] = std::min(_lazy[2 * v + 1], _lazy[v]);
        }
        return std::vector<T>(_lazy.begin() + _size, _lazy.begin() + _size + n);
    }
};

template <class T>
std::vector<T> replacement_paths_solve_edges(const ReplacementPathsData<T>& data) {
    int answer_size = int(data.path.edges.size());
    ReplacementPathsRangeChmin<T> range_chmin(answer_size, data.inf);
    for (const auto& e : data.edge_by_id) {
        if (!e.alive || data.is_path_edge[e.id]) continue;
        int u = e.from;
        int v = e.to;
        if (data.block[u] == -1 || data.block[v] == -1 || data.block[u] == data.block[v]) continue;
        if (data.block[u] > data.block[v]) std::swap(u, v);
        int a = data.block[u];
        int b = data.block[v];
        T candidate = replacement_paths_safe_add(data.dist_s[u], e.cost, data.inf);
        candidate = replacement_paths_safe_add(candidate, data.dist_t[v], data.inf);
        if (candidate == data.inf) continue;
        range_chmin.apply(a, b, candidate);
    }
    return range_chmin.values(answer_size);
}

template <class T>
T replacement_paths_without_vertex(const Graph<T>& g, int s, int t, int removed, T inf) {
    if (s == removed || t == removed) return inf;
    std::vector<T> dist(g.size(), inf);
    using P = std::pair<T, int>;
    std::priority_queue<P, std::vector<P>, std::greater<P>> que;
    dist[s] = T(0);
    que.emplace(T(0), s);
    while (!que.empty()) {
        auto [d, v] = que.top();
        que.pop();
        if (dist[v] != d) continue;
        for (const auto& e : g[v]) {
            if (!e.alive || e.to == removed) continue;
            T nd = replacement_paths_safe_add(d, e.cost, inf);
            if (dist[e.to] <= nd) continue;
            dist[e.to] = nd;
            que.emplace(nd, e.to);
        }
    }
    return dist[t];
}

template <class T>
std::vector<T> replacement_paths_solve_vertices(const Graph<T>& g,
                                                const ReplacementPathsData<T>& data) {
    // One edge can cross an edge cut, but a vertex-avoiding path may enter and
    // leave the failed vertex's tree block through two different detour edges.
    int path_size = int(data.path.vertices.size());
    std::vector<T> answer(path_size, data.inf);
    int s = data.path.vertices.front();
    int t = data.path.vertices.back();
    for (int i = 1; i + 1 < path_size; i++) {
        answer[i] = replacement_paths_without_vertex(g, s, t, data.path.vertices[i], data.inf);
    }
    return answer;
}

}  // namespace internal

template <class T>
EdgeReplacementPathsResult<T> edge_replacement_paths(
    const Graph<T>& g, const GraphPath& path, T inf = std::numeric_limits<T>::max() / T(4)) {
    assert(!path.vertices.empty());
    auto data = internal::replacement_paths_prepare(
        g, path, inf, static_cast<const DijkstraResult<T>*>(nullptr));
    auto replacement_dist = internal::replacement_paths_solve_edges(data);
    return {path, replacement_dist, inf};
}

template <class T>
EdgeReplacementPathsResult<T> edge_replacement_paths(
    const Graph<T>& g, int s, int t, T inf = std::numeric_limits<T>::max() / T(4)) {
    assert(0 <= s && s < g.size());
    assert(0 <= t && t < g.size());
    auto from_s = internal::replacement_paths_dijkstra(g, s, inf);
    assert(from_s.reachable(t));
    auto path = internal::replacement_paths_restore_path(from_s, s, t);
    auto data = internal::replacement_paths_prepare(g, path, inf, &from_s);
    auto replacement_dist = internal::replacement_paths_solve_edges(data);
    return {path, replacement_dist, inf};
}

template <class T>
VertexReplacementPathsResult<T> vertex_replacement_paths(
    const Graph<T>& g, const GraphPath& path, T inf = std::numeric_limits<T>::max() / T(4)) {
    assert(!path.vertices.empty());
    auto data = internal::replacement_paths_prepare(
        g, path, inf, static_cast<const DijkstraResult<T>*>(nullptr));
    auto replacement_dist = internal::replacement_paths_solve_vertices(g, data);
    return {path, replacement_dist, inf};
}

template <class T>
VertexReplacementPathsResult<T> vertex_replacement_paths(
    const Graph<T>& g, int s, int t, T inf = std::numeric_limits<T>::max() / T(4)) {
    assert(0 <= s && s < g.size());
    assert(0 <= t && t < g.size());
    auto from_s = internal::replacement_paths_dijkstra(g, s, inf);
    assert(from_s.reachable(t));
    auto path = internal::replacement_paths_restore_path(from_s, s, t);
    auto data = internal::replacement_paths_prepare(g, path, inf, &from_s);
    auto replacement_dist = internal::replacement_paths_solve_vertices(g, data);
    return {path, replacement_dist, inf};
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/shortest_path.hpp"



#line 1 "graph/bellman_ford.hpp"



#line 9 "graph/bellman_ford.hpp"

#line 11 "graph/bellman_ford.hpp"

namespace m1une {
namespace graph {

template <class T>
struct BellmanFordResult {
    std::vector<T> dist;
    std::vector<int> parent;
    std::vector<int> parent_edge;
    std::vector<bool> negative;
    T inf;
    bool has_negative_cycle;

    bool reachable(int v) const {
        assert(0 <= v && v < int(dist.size()));
        return dist[v] != inf;
    }

    bool affected_by_negative_cycle(int v) const {
        assert(0 <= v && v < int(negative.size()));
        return negative[v];
    }

    std::vector<int> path(int t) const {
        assert(reachable(t));
        assert(!affected_by_negative_cycle(t));
        std::vector<int> result;
        for (int v = t; v != -1; v = parent[v]) result.push_back(v);
        std::reverse(result.begin(), result.end());
        return result;
    }
};

template <class T>
BellmanFordResult<T> bellman_ford(const Graph<T>& g, const std::vector<int>& sources,
                                  T inf = std::numeric_limits<T>::max() / T(4)) {
    int n = g.size();
    BellmanFordResult<T> result;
    result.dist.assign(n, inf);
    result.parent.assign(n, -1);
    result.parent_edge.assign(n, -1);
    result.negative.assign(n, false);
    result.inf = inf;
    result.has_negative_cycle = false;

    for (int s : sources) {
        assert(0 <= s && s < n);
        result.dist[s] = T(0);
    }

    std::vector<int> relaxed_vertices;
    for (int iter = 0; iter < n; iter++) {
        bool updated = false;
        for (int v = 0; v < n; v++) {
            if (result.dist[v] == inf) continue;
            for (const auto& e : g[v]) {
                if (!e.alive) continue;
                T nd = result.dist[v] + e.cost;
                if (result.dist[e.to] <= nd) continue;
                result.dist[e.to] = nd;
                result.parent[e.to] = v;
                result.parent_edge[e.to] = e.id;
                updated = true;
                if (iter == n - 1) relaxed_vertices.push_back(e.to);
            }
        }
        if (!updated) break;
    }

    std::queue<int> que;
    for (int v : relaxed_vertices) {
        if (result.negative[v]) continue;
        result.negative[v] = true;
        que.push(v);
    }
    while (!que.empty()) {
        int v = que.front();
        que.pop();
        for (const auto& e : g[v]) {
            if (!e.alive) continue;
            if (result.negative[e.to]) continue;
            result.negative[e.to] = true;
            que.push(e.to);
        }
    }

    for (bool x : result.negative) result.has_negative_cycle = result.has_negative_cycle || x;
    return result;
}

template <class T>
BellmanFordResult<T> bellman_ford(const Graph<T>& g, int s, T inf = std::numeric_limits<T>::max() / T(4)) {
    return bellman_ford(g, std::vector<int>{s}, inf);
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/bfs.hpp"



#line 11 "graph/bfs.hpp"

#line 13 "graph/bfs.hpp"

namespace m1une {
namespace graph {

struct BfsResult {
    std::vector<int> dist;
    std::vector<int> parent;
    std::vector<int> parent_edge;

    bool reachable(int v) const {
        assert(0 <= v && v < int(dist.size()));
        return dist[v] != -1;
    }

    std::vector<int> path(int t) const {
        assert(reachable(t));
        std::vector<int> result;
        for (int v = t; v != -1; v = parent[v]) result.push_back(v);
        std::reverse(result.begin(), result.end());
        return result;
    }
};

namespace bfs_detail {

template <class Callback>
concept BfsCallback =
    std::invocable<Callback&, int, int> ||
    std::invocable<Callback&, int>;

template <BfsCallback Callback>
void invoke_callback(Callback& callback, int vertex, int parent) {
    if constexpr (std::invocable<Callback&, int, int>) {
        std::invoke(callback, vertex, parent);
    } else {
        std::invoke(callback, vertex);
    }
}

template <class T, class Callback>
BfsResult run_bfs(
    const Graph<T>& g,
    const std::vector<int>& sources,
    Callback& callback
) {
    int n = g.size();
    BfsResult result;
    result.dist.assign(n, -1);
    result.parent.assign(n, -1);
    result.parent_edge.assign(n, -1);

    std::queue<int> que;
    for (int s : sources) {
        assert(0 <= s && s < n);
        if (result.dist[s] != -1) continue;
        result.dist[s] = 0;
        invoke_callback(callback, s, -1);
        que.push(s);
    }

    while (!que.empty()) {
        int v = que.front();
        que.pop();
        for (const auto& e : g[v]) {
            if (!e.alive) continue;
            if (result.dist[e.to] != -1) continue;
            result.dist[e.to] = result.dist[v] + 1;
            result.parent[e.to] = v;
            result.parent_edge[e.to] = e.id;
            invoke_callback(callback, e.to, v);
            que.push(e.to);
        }
    }

    return result;
}

}  // namespace bfs_detail

template <class T>
BfsResult bfs(const Graph<T>& g, const std::vector<int>& sources) {
    auto callback = [](int) {};
    return bfs_detail::run_bfs(g, sources, callback);
}

template <class T>
BfsResult bfs(const Graph<T>& g, int s) {
    return bfs(g, std::vector<int>{s});
}

template <class T, class Callback>
requires bfs_detail::BfsCallback<Callback>
BfsResult bfs(
    const Graph<T>& g,
    const std::vector<int>& sources,
    Callback&& callback
) {
    return bfs_detail::run_bfs(g, sources, callback);
}

template <class T, class Callback>
requires bfs_detail::BfsCallback<Callback>
BfsResult bfs(const Graph<T>& g, int source, Callback&& callback) {
    return bfs(
        g,
        std::vector<int>{source},
        std::forward<Callback>(callback)
    );
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/cow_game.hpp"



#line 10 "graph/cow_game.hpp"

namespace m1une {
namespace graph {

template <class T>
struct CowGameConstraint {
    int a;
    int b;
    T upper_bound;
};

template <class T>
struct CowGameSolution {
    bool feasible = false;
    std::vector<T> value;

    bool is_feasible() const {
        return feasible;
    }
};

template <class T>
struct CowGameUpperBounds {
    bool feasible;
    std::vector<T> upper_bound;
    T inf;

    bool is_feasible() const {
        return feasible;
    }

    bool bounded(int variable) const {
        assert(0 <= variable && variable < int(upper_bound.size()));
        return feasible && upper_bound[variable] != inf;
    }
};

template <class T>
struct CowGameDifferenceBounds {
    bool feasible;
    std::optional<T> lower_bound;
    std::optional<T> upper_bound;

    bool is_feasible() const {
        return feasible;
    }

    bool bounded_below() const {
        return feasible && lower_bound.has_value();
    }

    bool bounded_above() const {
        return feasible && upper_bound.has_value();
    }
};

template <class T>
class CowGame {
    static_assert(std::is_arithmetic_v<T> && std::is_signed_v<T>);

    struct RelaxationResult {
        bool has_negative_cycle;
        std::vector<T> dist;
    };

    int _n;
    std::vector<CowGameConstraint<T>> _constraints;
    std::vector<std::vector<int>> _outgoing_constraints;
    bool _has_negative_upper_bound = false;
    mutable bool _solution_cached = false;
    mutable CowGameSolution<T> _cached_solution;

    void assert_variable(int variable) const {
        (void)variable;
        assert(0 <= variable && variable < _n);
    }

    T negate(T value) const {
        assert(value != std::numeric_limits<T>::lowest());
        return -value;
    }

    RelaxationResult check_feasibility() const {
        std::vector<T> dist(_n, T());
        for (int iteration = 0; iteration < _n; iteration++) {
            bool updated = false;
            for (const auto& constraint : _constraints) {
                T candidate = dist[constraint.b] + constraint.upper_bound;
                if (dist[constraint.a] <= candidate) continue;
                dist[constraint.a] = candidate;
                updated = true;
                if (iteration == _n - 1) return RelaxationResult{true, std::move(dist)};
            }
            if (!updated) break;
        }
        return RelaxationResult{false, std::move(dist)};
    }

    std::vector<T> shortest_paths(int source, T inf) const {
        const auto& potential = _cached_solution.value;
        std::vector<T> dist(_n, inf);
        std::vector<int> heap;
        // -1 is unseen, -2 is fixed, and every other value is a heap index.
        std::vector<int> position(_n, -1);
        heap.reserve(_n);

        auto swap_heap = [&](int i, int j) {
            std::swap(heap[i], heap[j]);
            position[heap[i]] = i;
            position[heap[j]] = j;
        };
        auto sift_up = [&](int i) {
            while (i > 0) {
                int parent = (i - 1) / 2;
                if (dist[heap[parent]] <= dist[heap[i]]) break;
                swap_heap(parent, i);
                i = parent;
            }
        };
        auto sift_down = [&](int i) {
            while (2 * i + 1 < int(heap.size())) {
                int child = 2 * i + 1;
                if (child + 1 < int(heap.size()) &&
                    dist[heap[child + 1]] < dist[heap[child]]) {
                    child++;
                }
                if (dist[heap[i]] <= dist[heap[child]]) break;
                swap_heap(i, child);
                i = child;
            }
        };

        dist[source] = T();
        position[source] = 0;
        heap.push_back(source);

        while (!heap.empty()) {
            int b = heap[0];
            position[b] = -2;
            int last = heap.back();
            heap.pop_back();
            if (!heap.empty()) {
                heap[0] = last;
                position[last] = 0;
                sift_down(0);
            }

            for (int id : _outgoing_constraints[b]) {
                const auto& constraint = _constraints[id];
                T cost = constraint.upper_bound + potential[b] -
                         potential[constraint.a];
                assert(cost >= T());
                T candidate = dist[b] + cost;
                if (dist[constraint.a] <= candidate) continue;
                dist[constraint.a] = candidate;
                assert(position[constraint.a] != -2);
                if (position[constraint.a] == -1) {
                    position[constraint.a] = int(heap.size());
                    heap.push_back(constraint.a);
                }
                sift_up(position[constraint.a]);
            }
        }

        for (int v = 0; v < _n; v++) {
            if (dist[v] == inf) continue;
            dist[v] = dist[v] - potential[source] + potential[v];
        }
        return dist;
    }

   public:
    CowGame() : CowGame(0) {}

    explicit CowGame(int variable_count)
        : _n(variable_count),
          _outgoing_constraints(variable_count < 0 ? 0 : variable_count) {
        assert(variable_count >= 0);
    }

    int size() const {
        return _n;
    }

    int constraint_count() const {
        return int(_constraints.size());
    }

    const CowGameConstraint<T>& get_constraint(int id) const {
        assert(0 <= id && id < int(_constraints.size()));
        return _constraints[id];
    }

    const std::vector<CowGameConstraint<T>>& constraints() const {
        return _constraints;
    }

    bool can_use_dijkstra() const {
        return !_has_negative_upper_bound ||
               (_solution_cached && _cached_solution.feasible);
    }

    int add_upper_bound(int a, int b, T upper_bound) {
        assert_variable(a);
        assert_variable(b);
        int id = int(_constraints.size());
        _constraints.push_back(CowGameConstraint<T>{a, b, upper_bound});
        _outgoing_constraints[b].push_back(id);
        _has_negative_upper_bound = _has_negative_upper_bound || upper_bound < T();
        _solution_cached = false;
        return id;
    }

    int add_constraint(int a, int b, T upper_bound) {
        return add_upper_bound(a, b, upper_bound);
    }

    int add_lower_bound(int a, int b, T lower_bound) {
        return add_upper_bound(b, a, negate(lower_bound));
    }

    void add_bounds(int a, int b, T lower_bound, T upper_bound) {
        assert(lower_bound <= upper_bound);
        add_lower_bound(a, b, lower_bound);
        add_upper_bound(a, b, upper_bound);
    }

    void add_equality(int a, int b, T difference) {
        add_bounds(a, b, difference, difference);
    }

    CowGameSolution<T> solve() const {
        if (_solution_cached) return _cached_solution;

        _cached_solution.feasible = true;
        _cached_solution.value.assign(_n, T());
        if (_has_negative_upper_bound) {
            auto result = check_feasibility();
            _cached_solution.feasible = !result.has_negative_cycle;
            _cached_solution.value.clear();
            if (_cached_solution.feasible) {
                _cached_solution.value = std::move(result.dist);
            }
        }
        _solution_cached = true;
        return _cached_solution;
    }

    bool is_feasible() const {
        if (!_solution_cached) (void)solve();
        return _cached_solution.feasible;
    }

    CowGameUpperBounds<T> tightest_upper_bounds(int source) const {
        assert_variable(source);
        T inf = std::numeric_limits<T>::max() / T(4);
        CowGameUpperBounds<T> result;
        result.feasible = is_feasible();
        result.inf = inf;
        result.upper_bound.assign(_n, inf);
        if (!result.feasible) return result;

        result.upper_bound = shortest_paths(source, inf);
        return result;
    }

    CowGameDifferenceBounds<T> difference_bounds(int a, int b) const {
        assert_variable(a);
        assert_variable(b);
        T inf = std::numeric_limits<T>::max() / T(4);
        CowGameDifferenceBounds<T> result;
        result.feasible = is_feasible();
        if (!result.feasible) return result;

        auto upper = shortest_paths(b, inf);
        if (upper[a] != inf) result.upper_bound = upper[a];

        auto lower = shortest_paths(a, inf);
        if (lower[b] != inf) result.lower_bound = negate(lower[b]);
        return result;
    }
};

template <class T>
using DifferenceConstraints = CowGame<T>;

}  // namespace graph
}  // namespace m1une


#line 1 "graph/dag_shortest_path.hpp"



#line 9 "graph/dag_shortest_path.hpp"

#line 1 "graph/topological_sort.hpp"



#line 7 "graph/topological_sort.hpp"

#line 9 "graph/topological_sort.hpp"

namespace m1une {
namespace graph {

template <class T>
std::optional<std::vector<int>> topological_sort(const Graph<T>& g) {
    int n = g.size();
    std::vector<int> indeg(n, 0);
    for (int v = 0; v < n; v++) {
        for (const auto& e : g[v]) {
            if (!e.alive) continue;
            indeg[e.to]++;
        }
    }

    std::queue<int> que;
    for (int v = 0; v < n; v++) {
        if (indeg[v] == 0) que.push(v);
    }

    std::vector<int> order;
    order.reserve(n);
    while (!que.empty()) {
        int v = que.front();
        que.pop();
        order.push_back(v);
        for (const auto& e : g[v]) {
            if (!e.alive) continue;
            indeg[e.to]--;
            if (indeg[e.to] == 0) que.push(e.to);
        }
    }

    if (int(order.size()) != n) return std::nullopt;
    return order;
}

template <class T>
bool is_dag(const Graph<T>& g) {
    return topological_sort(g).has_value();
}

}  // namespace graph
}  // namespace m1une


#line 12 "graph/dag_shortest_path.hpp"

namespace m1une {
namespace graph {

template <class T>
struct DagShortestPathResult {
    std::vector<T> dist;
    std::vector<int> parent;
    std::vector<int> parent_edge;
    std::vector<int> topological_order;
    T inf;

    bool reachable(int v) const {
        assert(0 <= v && v < int(dist.size()));
        return dist[v] != inf;
    }

    std::vector<int> path(int t) const {
        assert(reachable(t));
        std::vector<int> result;
        for (int v = t; v != -1; v = parent[v]) result.push_back(v);
        std::reverse(result.begin(), result.end());
        return result;
    }
};

template <class T>
std::optional<DagShortestPathResult<T>> dag_shortest_path(
    const Graph<T>& g, const std::vector<int>& sources, T inf = std::numeric_limits<T>::max() / T(4)) {
    int n = g.size();
    auto order = topological_sort(g);
    if (!order) return std::nullopt;

    DagShortestPathResult<T> result;
    result.dist.assign(n, inf);
    result.parent.assign(n, -1);
    result.parent_edge.assign(n, -1);
    result.topological_order = *order;
    result.inf = inf;

    for (int s : sources) {
        assert(0 <= s && s < n);
        if (result.dist[s] == T(0)) continue;
        result.dist[s] = T(0);
    }

    for (int v : *order) {
        if (result.dist[v] == inf) continue;
        for (const auto& e : g[v]) {
            if (!e.alive) continue;
            T nd = result.dist[v] + e.cost;
            if (result.dist[e.to] <= nd) continue;
            result.dist[e.to] = nd;
            result.parent[e.to] = v;
            result.parent_edge[e.to] = e.id;
        }
    }

    return result;
}

template <class T>
std::optional<DagShortestPathResult<T>> dag_shortest_path(
    const Graph<T>& g, int s, T inf = std::numeric_limits<T>::max() / T(4)) {
    return dag_shortest_path(g, std::vector<int>{s}, inf);
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/k_shortest_walk.hpp"



#line 10 "graph/k_shortest_walk.hpp"

#line 12 "graph/k_shortest_walk.hpp"

namespace m1une {
namespace graph {

namespace internal {

template <class T>
class KShortestWalkHeap {
    struct Node {
        T key;
        int to;
        int left;
        int right;
        int rank;
    };

    std::vector<Node> _nodes;

    int rank(int root) const {
        return root == -1 ? 0 : _nodes[root].rank;
    }

   public:
    int make_node(T key, int to) {
        int result = int(_nodes.size());
        _nodes.push_back(Node{key, to, -1, -1, 1});
        return result;
    }

    int meld_mutable(int first, int second) {
        if (first == -1) return second;
        if (second == -1) return first;
        if (_nodes[second].key < _nodes[first].key) std::swap(first, second);
        _nodes[first].right = meld_mutable(_nodes[first].right, second);
        if (rank(_nodes[first].left) < rank(_nodes[first].right)) {
            std::swap(_nodes[first].left, _nodes[first].right);
        }
        _nodes[first].rank = rank(_nodes[first].right) + 1;
        return first;
    }

    int meld_persistent(int first, int second) {
        if (first == -1) return second;
        if (second == -1) return first;
        if (_nodes[second].key < _nodes[first].key) std::swap(first, second);
        int result = int(_nodes.size());
        _nodes.push_back(_nodes[first]);
        _nodes[result].right = meld_persistent(_nodes[result].right, second);
        if (rank(_nodes[result].left) < rank(_nodes[result].right)) {
            std::swap(_nodes[result].left, _nodes[result].right);
        }
        _nodes[result].rank = rank(_nodes[result].right) + 1;
        return result;
    }

    const Node& operator[](int index) const {
        return _nodes[index];
    }
};

}  // namespace internal

template <class T>
std::vector<T> k_shortest_walk(
    const Graph<T>& g,
    int s,
    int t,
    int k,
    T inf = std::numeric_limits<T>::max() / T(4)
) {
    int n = g.size();
    assert(0 <= s && s < n);
    assert(0 <= t && t < n);
    assert(0 <= k);
    if (k == 0) return {};

    struct ReverseEdge {
        int from;
        int index;
        T cost;
    };
    std::vector<std::vector<ReverseEdge>> reverse_graph(n);
    for (int from = 0; from < n; from++) {
        for (int index = 0; index < int(g[from].size()); index++) {
            const auto& edge = g[from][index];
            if (!edge.alive) continue;
            assert(T(0) <= edge.cost);
            reverse_graph[edge.to].push_back(ReverseEdge{from, index, edge.cost});
        }
    }

    std::vector<T> dist(n, inf);
    std::vector<int> tree_edge(n, -1);
    std::vector<int> order;
    order.reserve(n);
    using QueueEntry = std::pair<T, int>;
    std::priority_queue<QueueEntry, std::vector<QueueEntry>, std::greater<QueueEntry>> queue;
    dist[t] = T(0);
    queue.emplace(T(0), t);
    while (!queue.empty()) {
        auto [current_dist, vertex] = queue.top();
        queue.pop();
        if (dist[vertex] != current_dist) continue;
        order.push_back(vertex);
        for (const auto& edge : reverse_graph[vertex]) {
            T next_dist = current_dist + edge.cost;
            if (dist[edge.from] <= next_dist) continue;
            dist[edge.from] = next_dist;
            tree_edge[edge.from] = edge.index;
            queue.emplace(next_dist, edge.from);
        }
    }
    if (dist[s] == inf) return {};

    internal::KShortestWalkHeap<T> heap_pool;
    std::vector<int> local_heap(n, -1);
    for (int vertex : order) {
        for (int index = 0; index < int(g[vertex].size()); index++) {
            const auto& edge = g[vertex][index];
            if (!edge.alive || dist[edge.to] == inf || index == tree_edge[vertex]) continue;
            T extra = edge.cost + dist[edge.to] - dist[vertex];
            assert(T(0) <= extra);
            int node = heap_pool.make_node(extra, edge.to);
            local_heap[vertex] = heap_pool.meld_mutable(local_heap[vertex], node);
        }
    }

    std::vector<int> path_heap(n, -1);
    for (int vertex : order) {
        int inherited = -1;
        if (tree_edge[vertex] != -1) inherited = path_heap[g[vertex][tree_edge[vertex]].to];
        path_heap[vertex] = heap_pool.meld_persistent(inherited, local_heap[vertex]);
    }

    std::vector<T> result;
    result.reserve(k);
    result.push_back(dist[s]);
    std::priority_queue<QueueEntry, std::vector<QueueEntry>, std::greater<QueueEntry>> candidates;
    if (path_heap[s] != -1) {
        candidates.emplace(dist[s] + heap_pool[path_heap[s]].key, path_heap[s]);
    }
    while (int(result.size()) < k && !candidates.empty()) {
        auto [cost, node_index] = candidates.top();
        candidates.pop();
        result.push_back(cost);
        const auto& node = heap_pool[node_index];
        if (node.left != -1) {
            candidates.emplace(cost - node.key + heap_pool[node.left].key, node.left);
        }
        if (node.right != -1) {
            candidates.emplace(cost - node.key + heap_pool[node.right].key, node.right);
        }
        int next_heap = path_heap[node.to];
        if (next_heap != -1) {
            candidates.emplace(cost + heap_pool[next_heap].key, next_heap);
        }
    }
    return result;
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/warshall_floyd.hpp"



#line 8 "graph/warshall_floyd.hpp"

#line 10 "graph/warshall_floyd.hpp"

namespace m1une {
namespace graph {

template <class T>
std::vector<std::vector<T>> warshall_floyd(std::vector<std::vector<T>> dist,
                                           T inf = std::numeric_limits<T>::max() / T(4)) {
    int n = int(dist.size());
    for (int k = 0; k < n; k++) {
        for (int i = 0; i < n; i++) {
            if (dist[i][k] == inf) continue;
            for (int j = 0; j < n; j++) {
                if (dist[k][j] == inf) continue;
                T nd = dist[i][k] + dist[k][j];
                if (nd < dist[i][j]) dist[i][j] = nd;
            }
        }
    }
    return dist;
}

template <class T>
std::vector<std::vector<T>> warshall_floyd(const Graph<T>& g, T inf = std::numeric_limits<T>::max() / T(4)) {
    int n = g.size();
    std::vector<std::vector<T>> dist(n, std::vector<T>(n, inf));
    for (int i = 0; i < n; i++) dist[i][i] = T(0);
    for (int v = 0; v < n; v++) {
        for (const auto& e : g[v]) {
            if (!e.alive) continue;
            if (e.cost < dist[e.from][e.to]) dist[e.from][e.to] = e.cost;
        }
    }
    return warshall_floyd(std::move(dist), inf);
}

template <class T>
bool warshall_floyd_add_directed_edge(std::vector<std::vector<T>>& dist, int from, int to, T cost,
                                      T inf = std::numeric_limits<T>::max() / T(4)) {
    int n = int(dist.size());
    assert(0 <= from && from < n);
    assert(0 <= to && to < n);

    std::vector<T> to_from(n), from_to(n);
    for (int i = 0; i < n; i++) {
        to_from[i] = dist[i][from];
        from_to[i] = dist[to][i];
    }

    bool updated = false;
    for (int i = 0; i < n; i++) {
        if (to_from[i] == inf) continue;
        for (int j = 0; j < n; j++) {
            if (from_to[j] == inf) continue;
            T nd = to_from[i] + cost + from_to[j];
            if (nd < dist[i][j]) {
                dist[i][j] = nd;
                updated = true;
            }
        }
    }
    return updated;
}

template <class T>
bool warshall_floyd_add_undirected_edge(std::vector<std::vector<T>>& dist, int u, int v, T cost,
                                        T inf = std::numeric_limits<T>::max() / T(4)) {
    int n = int(dist.size());
    assert(0 <= u && u < n);
    assert(0 <= v && v < n);

    std::vector<T> to_u(n), from_u(n), to_v(n), from_v(n);
    for (int i = 0; i < n; i++) {
        to_u[i] = dist[i][u];
        from_u[i] = dist[u][i];
        to_v[i] = dist[i][v];
        from_v[i] = dist[v][i];
    }

    bool updated = false;
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            if (to_u[i] != inf && from_v[j] != inf) {
                T nd = to_u[i] + cost + from_v[j];
                if (nd < dist[i][j]) {
                    dist[i][j] = nd;
                    updated = true;
                }
            }
            if (to_v[i] != inf && from_u[j] != inf) {
                T nd = to_v[i] + cost + from_u[j];
                if (nd < dist[i][j]) {
                    dist[i][j] = nd;
                    updated = true;
                }
            }
        }
    }
    return updated;
}

template <class T>
bool has_negative_cycle(const std::vector<std::vector<T>>& dist) {
    int n = int(dist.size());
    for (int i = 0; i < n; i++) {
        if (dist[i][i] < T(0)) return true;
    }
    return false;
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/zero_one_bfs.hpp"



#line 6 "graph/zero_one_bfs.hpp"
#include <deque>
#line 9 "graph/zero_one_bfs.hpp"

#line 11 "graph/zero_one_bfs.hpp"

namespace m1une {
namespace graph {

struct ZeroOneBfsResult {
    std::vector<int> dist;
    std::vector<int> parent;
    std::vector<int> parent_edge;
    int inf;

    bool reachable(int v) const {
        assert(0 <= v && v < int(dist.size()));
        return dist[v] != inf;
    }

    std::vector<int> path(int t) const {
        assert(reachable(t));
        std::vector<int> result;
        for (int v = t; v != -1; v = parent[v]) result.push_back(v);
        std::reverse(result.begin(), result.end());
        return result;
    }
};

template <class T>
ZeroOneBfsResult zero_one_bfs(const Graph<T>& g, const std::vector<int>& sources,
                              int inf = std::numeric_limits<int>::max() / 2) {
    int n = g.size();
    ZeroOneBfsResult result;
    result.dist.assign(n, inf);
    result.parent.assign(n, -1);
    result.parent_edge.assign(n, -1);
    result.inf = inf;

    std::deque<int> deq;
    for (int s : sources) {
        assert(0 <= s && s < n);
        if (result.dist[s] == 0) continue;
        result.dist[s] = 0;
        deq.push_back(s);
    }

    while (!deq.empty()) {
        int v = deq.front();
        deq.pop_front();
        for (const auto& e : g[v]) {
            if (!e.alive) continue;
            int w;
            if (e.cost == T(0)) {
                w = 0;
            } else {
                assert(e.cost == T(1));
                w = 1;
            }
            int nd = result.dist[v] + w;
            if (result.dist[e.to] <= nd) continue;
            result.dist[e.to] = nd;
            result.parent[e.to] = v;
            result.parent_edge[e.to] = e.id;
            if (w == 0) {
                deq.push_front(e.to);
            } else {
                deq.push_back(e.to);
            }
        }
    }

    return result;
}

template <class T>
ZeroOneBfsResult zero_one_bfs(const Graph<T>& g, int s, int inf = std::numeric_limits<int>::max() / 2) {
    return zero_one_bfs(g, std::vector<int>{s}, inf);
}

}  // namespace graph
}  // namespace m1une


#line 12 "graph/shortest_path.hpp"


#line 1 "graph/st_numbering.hpp"



#line 6 "graph/st_numbering.hpp"

#line 8 "graph/st_numbering.hpp"

namespace m1une {
namespace graph {

// Returns ranks p with p[source] = 0 and p[sink] = n - 1 such that every
// other vertex has neighbors of both smaller and larger rank. Returns an empty
// vector when no such numbering exists.
template <class T>
std::vector<int> st_numbering(
    const Graph<T>& graph,
    int source,
    int sink
) {
    const int n = graph.size();
    assert(0 < n);
    assert(0 <= source && source < n);
    assert(0 <= sink && sink < n);
    assert(source != sink);

#ifndef NDEBUG
    std::vector<int> incidence_count(graph.edge_count(), 0);
    for (int vertex = 0; vertex < n; vertex++) {
        for (const Edge<T>& edge : graph[vertex]) {
            if (!edge.alive) continue;
            assert(0 <= edge.id && edge.id < graph.edge_count());
            incidence_count[edge.id]++;
        }
    }
    for (int edge_id = 0; edge_id < graph.edge_count(); edge_id++) {
        if (graph.is_edge_alive(edge_id)) {
            assert(incidence_count[edge_id] == 2);
        }
    }
#endif

    std::vector<int> parent(n, -1);
    std::vector<int> preorder(n, -1);
    std::vector<int> low_vertex(n, -1);
    std::vector<int> next_edge(n, 0);
    std::vector<int> traversal;
    traversal.reserve(n);

    preorder[source] = 0;
    low_vertex[source] = source;
    traversal.push_back(source);
    preorder[sink] = 1;
    low_vertex[sink] = sink;
    traversal.push_back(sink);

    std::vector<int> stack(1, sink);
    while (!stack.empty()) {
        const int vertex = stack.back();
        if (next_edge[vertex] < int(graph[vertex].size())) {
            const Edge<T>& edge = graph[vertex][next_edge[vertex]++];
            if (!edge.alive || edge.to == vertex) continue;
            const int to = edge.to;
            if (preorder[to] == -1) {
                parent[to] = vertex;
                preorder[to] = int(traversal.size());
                low_vertex[to] = to;
                traversal.push_back(to);
                stack.push_back(to);
            } else if (preorder[to] < preorder[low_vertex[vertex]]) {
                low_vertex[vertex] = to;
            }
            continue;
        }

        stack.pop_back();
        const int parent_vertex = parent[vertex];
        if (parent_vertex != -1 &&
            preorder[low_vertex[vertex]] <
                preorder[low_vertex[parent_vertex]]) {
            low_vertex[parent_vertex] = low_vertex[vertex];
        }
    }
    if (int(traversal.size()) != n) return {};

    std::vector<int> next(n, -1);
    std::vector<int> previous(n, -1);
    std::vector<int> sign(n, 0);
    next[source] = sink;
    previous[sink] = source;
    sign[source] = -1;

    for (int index = 2; index < n; index++) {
        const int vertex = traversal[index];
        const int parent_vertex = parent[vertex];
        assert(parent_vertex != -1);
        if (sign[low_vertex[vertex]] == -1) {
            const int before = previous[parent_vertex];
            if (before == -1) return {};
            next[before] = vertex;
            next[vertex] = parent_vertex;
            previous[vertex] = before;
            previous[parent_vertex] = vertex;
            sign[parent_vertex] = 1;
        } else {
            const int after = next[parent_vertex];
            if (after == -1) return {};
            next[parent_vertex] = vertex;
            next[vertex] = after;
            previous[vertex] = parent_vertex;
            previous[after] = vertex;
            sign[parent_vertex] = -1;
        }
    }

    std::vector<int> order;
    order.reserve(n);
    int vertex = source;
    while (vertex != -1 && int(order.size()) <= n) {
        order.push_back(vertex);
        if (vertex == sink) break;
        vertex = next[vertex];
    }
    if (int(order.size()) != n || order.back() != sink) return {};

    std::vector<int> rank(n, -1);
    for (int index = 0; index < n; index++) rank[order[index]] = index;

    for (int index = 0; index < n; index++) {
        const int current = order[index];
        bool has_smaller = false;
        bool has_larger = false;
        for (const Edge<T>& edge : graph[current]) {
            if (!edge.alive || edge.to == current) continue;
            has_smaller = has_smaller || rank[edge.to] < index;
            has_larger = has_larger || index < rank[edge.to];
        }
        if (index > 0 && !has_smaller) return {};
        if (index + 1 < n && !has_larger) return {};
    }
    return rank;
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/three_edge_connected_components.hpp"



#line 9 "graph/three_edge_connected_components.hpp"

#line 11 "graph/three_edge_connected_components.hpp"

namespace m1une {
namespace graph {

struct ThreeEdgeConnectedComponentsResult {
    std::vector<std::vector<int>> components;
    std::vector<int> component_of_vertex;

    int component_count() const {
        return int(components.size());
    }

    bool same(int first, int second) const {
        assert(0 <= first && first < int(component_of_vertex.size()));
        assert(0 <= second && second < int(component_of_vertex.size()));
        return component_of_vertex[first] == component_of_vertex[second];
    }
};

namespace internal {

// Maintains every component as a circular linked list. Swapping two successors
// concatenates two different lists in O(1) time.
struct ThreeEdgeComponentCycles {
    std::vector<int> next;

    explicit ThreeEdgeComponentCycles(int n) : next(n) {
        std::iota(next.begin(), next.end(), 0);
    }

    void unite(int first, int second) {
        std::swap(next[first], next[second]);
    }

    ThreeEdgeConnectedComponentsResult build_result() const {
        const int n = int(next.size());
        ThreeEdgeConnectedComponentsResult result;
        result.component_of_vertex.assign(n, -1);
        for (int first = 0; first < n; first++) {
            if (result.component_of_vertex[first] != -1) continue;
            const int component = result.component_count();
            result.components.emplace_back();
            int vertex = first;
            do {
                result.component_of_vertex[vertex] = component;
                result.components.back().push_back(vertex);
                vertex = next[vertex];
            } while (vertex != first);
        }
        return result;
    }
};

}  // namespace internal

// Decomposes an undirected multigraph into maximal vertex sets joined by at
// least three edge-disjoint paths. This is an iterative form of Tsin's
// one-pass contraction algorithm.
template <class T>
ThreeEdgeConnectedComponentsResult three_edge_connected_components(
    const Graph<T>& graph
) {
    const int n = graph.size();
    const int edge_count = graph.edge_count();

#ifndef NDEBUG
    std::vector<int> incidence_count(edge_count, 0);
    for (int vertex = 0; vertex < n; vertex++) {
        for (const Edge<T>& edge : graph[vertex]) {
            if (!edge.alive) continue;
            assert(edge.from == vertex);
            assert(0 <= edge.to && edge.to < n);
            assert(0 <= edge.id && edge.id < edge_count);
            incidence_count[edge.id]++;
        }
    }
    for (int edge_id = 0; edge_id < edge_count; edge_id++) {
        if (incidence_count[edge_id] != 0) assert(incidence_count[edge_id] == 2);
    }
#endif

    const int none = n;
    std::vector<int> enter(n, -1);
    std::vector<int> leave(n, 0);
    std::vector<int> low(n, none);
    std::vector<int> degree(n, 0);
    std::vector<int> path(n, none);
    std::vector<int> parent(n, -1);
    std::vector<int> parent_edge(n, -1);
    std::vector<int> next_edge(n, 0);
    std::vector<int> dfs_stack;
    internal::ThreeEdgeComponentCycles component_cycles(n);
    int timer = 0;

    auto absorb = [&](int vertex, int other) {
        component_cycles.unite(vertex, other);
        degree[vertex] += degree[other];
    };

    auto process_visited_edge = [&](int vertex, int to) {
        if (enter[to] < enter[vertex]) {
            degree[vertex]++;
            low[vertex] = std::min(low[vertex], enter[to]);
            return;
        }

        degree[vertex]--;
        int current = path[vertex];
        while (current != none && enter[current] <= enter[to] && enter[to] < leave[current]) {
            absorb(vertex, current);
            current = path[current];
        }
        path[vertex] = current;
    };

    auto process_child = [&](int vertex, int child) {
        if (path[child] == none && degree[child] <= 1) {
            degree[vertex] += degree[child];
            low[vertex] = std::min(low[vertex], low[child]);
            return;
        }

        int current = child;
        if (degree[child] == 0) current = path[child];
        assert(current != none);
        if (low[current] < low[vertex]) {
            low[vertex] = low[current];
            std::swap(current, path[vertex]);
        }
        while (current != none) {
            absorb(vertex, current);
            current = path[current];
        }
    };

    for (int root = 0; root < n; root++) {
        if (enter[root] != -1) continue;
        enter[root] = timer++;
        dfs_stack.push_back(root);

        while (!dfs_stack.empty()) {
            const int vertex = dfs_stack.back();
            if (next_edge[vertex] < int(graph[vertex].size())) {
                const Edge<T>& edge = graph[vertex][next_edge[vertex]++];
                if (!edge.alive || edge.from == edge.to || edge.id == parent_edge[vertex]) continue;
                const int to = edge.to;
                if (enter[to] == -1) {
                    parent[to] = vertex;
                    parent_edge[to] = edge.id;
                    enter[to] = timer++;
                    dfs_stack.push_back(to);
                } else {
                    process_visited_edge(vertex, to);
                }
                continue;
            }

            leave[vertex] = timer;
            dfs_stack.pop_back();
            if (parent[vertex] != -1) process_child(parent[vertex], vertex);
        }
    }

    return component_cycles.build_result();
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/two_edge_connected_components.hpp"



#line 6 "graph/two_edge_connected_components.hpp"

#line 8 "graph/two_edge_connected_components.hpp"

namespace m1une {
namespace graph {

struct TwoEdgeConnectedBridge {
    int from;
    int to;
    int edge_id;
};

struct TwoEdgeConnectedComponentsResult {
    std::vector<std::vector<int>> components;
    std::vector<int> component_of_vertex;
    std::vector<int> bridge_ids;
    std::vector<char> bridge;
    std::vector<TwoEdgeConnectedBridge> bridge_forest_edges;
    std::vector<int> ord;
    std::vector<int> low;

    int component_count() const {
        return int(components.size());
    }

    bool same(int first, int second) const {
        assert(0 <= first && first < int(component_of_vertex.size()));
        assert(0 <= second && second < int(component_of_vertex.size()));
        return component_of_vertex[first] == component_of_vertex[second];
    }

    bool is_bridge(int edge_id) const {
        assert(0 <= edge_id && edge_id < int(bridge.size()));
        return bridge[edge_id];
    }
};

// Removes every active bridge and returns the remaining connected components.
// The first lowlink traversal and the component traversal are both iterative.
template <class T>
TwoEdgeConnectedComponentsResult two_edge_connected_components(
    const Graph<T>& graph
) {
    const int n = graph.size();
    const int edge_count = graph.edge_count();

    TwoEdgeConnectedComponentsResult result;
    result.component_of_vertex.assign(n, -1);
    result.bridge.assign(edge_count, false);
    result.ord.assign(n, -1);
    result.low.assign(n, -1);

    std::vector<int> edge_from(edge_count, -1);
    std::vector<int> edge_to(edge_count, -1);
    std::vector<int> incidence_count(edge_count, 0);
    for (int vertex = 0; vertex < n; vertex++) {
        for (const Edge<T>& edge : graph[vertex]) {
            if (!edge.alive) continue;
            assert(0 <= edge.id && edge.id < edge_count);
            if (incidence_count[edge.id] == 0) {
                edge_from[edge.id] = edge.from;
                edge_to[edge.id] = edge.to;
            }
            incidence_count[edge.id]++;
        }
    }
#ifndef NDEBUG
    for (int edge_id = 0; edge_id < edge_count; edge_id++) {
        if (incidence_count[edge_id] != 0) assert(incidence_count[edge_id] == 2);
    }
#endif

    std::vector<int> parent(n, -1);
    std::vector<int> parent_edge(n, -1);
    std::vector<int> next_edge(n, 0);
    std::vector<int> stack;
    int timer = 0;

    for (int root = 0; root < n; root++) {
        if (result.ord[root] != -1) continue;
        result.ord[root] = result.low[root] = timer++;
        stack.push_back(root);
        while (!stack.empty()) {
            const int vertex = stack.back();
            if (next_edge[vertex] < int(graph[vertex].size())) {
                const Edge<T>& edge = graph[vertex][next_edge[vertex]++];
                if (!edge.alive || edge.id == parent_edge[vertex]) continue;
                const int to = edge.to;
                if (result.ord[to] == -1) {
                    parent[to] = vertex;
                    parent_edge[to] = edge.id;
                    result.ord[to] = result.low[to] = timer++;
                    stack.push_back(to);
                } else if (result.ord[to] < result.low[vertex]) {
                    result.low[vertex] = result.ord[to];
                }
                continue;
            }

            stack.pop_back();
            const int parent_vertex = parent[vertex];
            if (parent_vertex == -1) continue;
            if (result.low[vertex] < result.low[parent_vertex]) {
                result.low[parent_vertex] = result.low[vertex];
            }
            if (result.ord[parent_vertex] < result.low[vertex]) {
                result.bridge[parent_edge[vertex]] = true;
            }
        }
    }

    for (int root = 0; root < n; root++) {
        if (result.component_of_vertex[root] != -1) continue;
        const int component = result.component_count();
        result.components.emplace_back();
        result.component_of_vertex[root] = component;
        stack.push_back(root);
        while (!stack.empty()) {
            const int vertex = stack.back();
            stack.pop_back();
            result.components.back().push_back(vertex);
            for (const Edge<T>& edge : graph[vertex]) {
                if (!edge.alive || result.bridge[edge.id]) continue;
                if (result.component_of_vertex[edge.to] != -1) continue;
                result.component_of_vertex[edge.to] = component;
                stack.push_back(edge.to);
            }
        }
    }

    for (int edge_id = 0; edge_id < edge_count; edge_id++) {
        if (!result.bridge[edge_id]) continue;
        result.bridge_ids.push_back(edge_id);
        const int first_component = result.component_of_vertex[edge_from[edge_id]];
        const int second_component = result.component_of_vertex[edge_to[edge_id]];
        assert(first_component != second_component);
        result.bridge_forest_edges.push_back(
            TwoEdgeConnectedBridge{first_component, second_component, edge_id});
    }
    return result;
}

}  // namespace graph
}  // namespace m1une


#line 32 "graph/undirected.hpp"

Back to top page