Directed Graph Algorithms
(graph/directed.hpp)
- View this file on GitHub
- Last update: 2026-08-24 02:34:24+09:00
- Include:
#include "graph/directed.hpp"
Overview
graph/directed.hpp includes algorithms whose main interpretation is directed,
plus direction-respecting shortest paths.
Use this header when the input edges are one-way, or when reachability/order depends on edge direction.
Included Headers
| Header | Graph orientation | Contents |
|---|---|---|
graph/dag.hpp |
Directed DAG only | Ordering, shortest and longest paths, path counts, reachability, transitive reduction, and minimum path cover. |
graph/shortest_path.hpp |
Direction-respecting / DAG-specific | BFS, 0-1 BFS, DAG shortest path, Dijkstra, Bellman-Ford, and Warshall-Floyd. |
graph/dfs.hpp |
Direction-respecting | Iterative DFS forests with parent paths, timestamps, and traversal orders. |
graph/directed_mst.hpp |
Directed rooted graph | Minimum-cost spanning arborescence with edge reconstruction. |
graph/matrix_tree_theorem.hpp |
Directed rooted graph | Counts weighted inward and outward spanning arborescences. |
graph/scc.hpp |
Directed only | Strongly connected components and condensation DAG. |
graph/incremental_scc.hpp |
Directed only | Offline SCC merge times under edge insertions. |
graph/functional_graph.hpp |
One successor per vertex | Cycle decomposition, large jumps, paths and orbits, visit counts, reachability distances, and synchronized meetings. |
graph/two_sat.hpp |
Implication graph | 2-SAT clauses, satisfiability, and one assignment. |
graph/cycle_detection.hpp |
Directed and undirected variants | Use find_directed_cycle(g) for directed graphs. |
graph/eulerian_trail.hpp |
Directed and undirected variants | Use directed_eulerian_trail(g) for directed 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
- Bellman-Ford
(graph/bellman_ford.hpp)
- BFS
(graph/bfs.hpp)
- Bipartite Graph
(graph/bipartite.hpp)
- Cow Game (Difference Constraints)
(graph/cow_game.hpp)
- Cycle Detection
(graph/cycle_detection.hpp)
- DAG Algorithms
(graph/dag.hpp)
- DAG Longest Path
(graph/dag_longest_path.hpp)
- DAG Path Count
(graph/dag_path_count.hpp)
- Minimum DAG Path Cover
(graph/dag_path_cover.hpp)
- DAG Reachability and Transitive Reduction
(graph/dag_reachability.hpp)
- DAG Shortest Path
(graph/dag_shortest_path.hpp)
- DFS
(graph/dfs.hpp)
- Dijkstra
(graph/dijkstra.hpp)
- Directed Minimum Spanning Tree
(graph/directed_mst.hpp)
- Eulerian Trail
(graph/eulerian_trail.hpp)
- Functional Graph
(graph/functional_graph.hpp)
- Graph
(graph/graph.hpp)
- Incremental Strongly Connected Components
(graph/incremental_scc.hpp)
- K-Shortest Walk
(graph/k_shortest_walk.hpp)
- Matrix-Tree Theorem
(graph/matrix_tree_theorem.hpp)
- Strongly Connected Components
(graph/scc.hpp)
- Shortest Path
(graph/shortest_path.hpp)
- Topological Sort
(graph/topological_sort.hpp)
- Two-Satisfiability
(graph/two_sat.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)
- Dynamic Bitset
(utilities/dynamic_bitset.hpp)
Required by
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_DIRECTED_HPP
#define M1UNE_GRAPH_DIRECTED_HPP 1
#include "cycle_detection.hpp"
#include "dag.hpp"
#include "dfs.hpp"
#include "directed_mst.hpp"
#include "eulerian_trail.hpp"
#include "functional_graph.hpp"
#include "graph.hpp"
#include "incremental_scc.hpp"
#include "matrix_tree_theorem.hpp"
#include "scc.hpp"
#include "shortest_path.hpp"
#include "two_sat.hpp"
#endif // M1UNE_GRAPH_DIRECTED_HPP#line 1 "graph/directed.hpp"
#line 1 "graph/cycle_detection.hpp"
#include <algorithm>
#include <cstddef>
#include <vector>
#line 1 "graph/graph.hpp"
#include <array>
#include <cassert>
#include <utility>
#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 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/dag.hpp"
#line 1 "graph/dag_longest_path.hpp"
#line 6 "graph/dag_longest_path.hpp"
#include <limits>
#include <optional>
#line 9 "graph/dag_longest_path.hpp"
#line 1 "graph/topological_sort.hpp"
#line 5 "graph/topological_sort.hpp"
#include <queue>
#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_longest_path.hpp"
namespace m1une {
namespace graph {
template <class T>
struct DagLongestPathResult {
std::vector<T> dist;
std::vector<int> parent;
std::vector<int> parent_edge;
std::vector<int> topological_order;
T neg_inf;
bool reachable(int v) const {
assert(0 <= v && v < int(dist.size()));
return dist[v] != neg_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<DagLongestPathResult<T>> dag_longest_path(
const Graph<T>& g,
const std::vector<int>& sources,
T neg_inf = std::numeric_limits<T>::lowest() / T(4)
) {
const int n = g.size();
auto order = topological_sort(g);
if (!order) return std::nullopt;
DagLongestPathResult<T> result;
result.dist.assign(n, neg_inf);
result.parent.assign(n, -1);
result.parent_edge.assign(n, -1);
result.topological_order = *order;
result.neg_inf = neg_inf;
for (int s : sources) {
assert(0 <= s && s < n);
result.dist[s] = T(0);
}
for (int v : *order) {
if (result.dist[v] == neg_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<DagLongestPathResult<T>> dag_longest_path(
const Graph<T>& g,
int s,
T neg_inf = std::numeric_limits<T>::lowest() / T(4)
) {
return dag_longest_path(g, std::vector<int>{s}, neg_inf);
}
} // namespace graph
} // namespace m1une
#line 1 "graph/dag_path_count.hpp"
#line 7 "graph/dag_path_count.hpp"
#line 10 "graph/dag_path_count.hpp"
namespace m1une {
namespace graph {
template <class Count = long long, class T>
std::optional<std::vector<Count>> dag_path_count(
const Graph<T>& g,
const std::vector<int>& sources
) {
const int n = g.size();
auto order = topological_sort(g);
if (!order) return std::nullopt;
std::vector<Count> ways(n, Count(0));
std::vector<char> used_source(n, false);
for (int s : sources) {
assert(0 <= s && s < n);
if (used_source[s]) continue;
used_source[s] = true;
ways[s] += Count(1);
}
for (int v : *order) {
for (const auto& e : g[v]) {
if (e.alive) ways[e.to] += ways[v];
}
}
return ways;
}
template <class Count = long long, class T>
std::optional<std::vector<Count>> dag_path_count(const Graph<T>& g, int s) {
return dag_path_count<Count>(g, std::vector<int>{s});
}
} // namespace graph
} // namespace m1une
#line 1 "graph/dag_path_cover.hpp"
#line 7 "graph/dag_path_cover.hpp"
#line 1 "graph/bipartite.hpp"
#line 7 "graph/bipartite.hpp"
#include <cstdint>
#line 13 "graph/bipartite.hpp"
#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 11 "graph/dag_path_cover.hpp"
namespace m1une {
namespace graph {
struct DagPathCoverResult {
std::vector<std::vector<int>> paths;
std::vector<std::vector<int>> path_edge_ids;
std::vector<int> predecessor;
std::vector<int> successor;
std::vector<int> predecessor_edge;
std::vector<int> successor_edge;
int size() const {
return int(paths.size());
}
};
template <class T>
std::optional<DagPathCoverResult> minimum_dag_path_cover(const Graph<T>& g) {
const int n = g.size();
if (!topological_sort(g)) return std::nullopt;
BipartiteMatching matching(n, n);
std::vector<int> original_edge_id;
for (int v = 0; v < n; v++) {
for (const auto& e : g[v]) {
if (!e.alive) continue;
matching.add_edge(v, e.to);
original_edge_id.push_back(e.id);
}
}
DagPathCoverResult result;
result.predecessor.assign(n, -1);
result.successor.assign(n, -1);
result.predecessor_edge.assign(n, -1);
result.successor_edge.assign(n, -1);
for (const auto& pair : matching.matching()) {
const int edge_id = original_edge_id[pair.edge_id];
result.successor[pair.left] = pair.right;
result.successor_edge[pair.left] = edge_id;
result.predecessor[pair.right] = pair.left;
result.predecessor_edge[pair.right] = edge_id;
}
int covered = 0;
for (int s = 0; s < n; s++) {
if (result.predecessor[s] != -1) continue;
result.paths.emplace_back();
result.path_edge_ids.emplace_back();
for (int v = s; v != -1; v = result.successor[v]) {
result.paths.back().push_back(v);
covered++;
if (result.successor_edge[v] != -1) {
result.path_edge_ids.back().push_back(result.successor_edge[v]);
}
}
}
assert(covered == n);
return result;
}
} // namespace graph
} // namespace m1une
#line 1 "graph/dag_reachability.hpp"
#line 9 "graph/dag_reachability.hpp"
#line 1 "utilities/dynamic_bitset.hpp"
#line 9 "utilities/dynamic_bitset.hpp"
namespace m1une {
namespace utilities {
struct DynamicBitset {
private:
static constexpr int BITS_PER_BLOCK = 64;
static constexpr uint64_t FULL_BLOCK = ~uint64_t{0};
int _n;
std::vector<uint64_t> blocks;
static int block_count(int n) {
assert(n >= 0);
return (n + BITS_PER_BLOCK - 1) >> 6;
}
uint64_t tail_mask() const {
const int rem = _n & (BITS_PER_BLOCK - 1);
return rem == 0 ? FULL_BLOCK : ((uint64_t{1} << rem) - 1);
}
// Keep unused bits in the last block equal to zero.
void clean() {
if (!blocks.empty()) blocks.back() &= tail_mask();
}
public:
DynamicBitset() : _n(0), blocks() {}
explicit DynamicBitset(int n, bool val = false) : _n(n), blocks(block_count(n), val ? FULL_BLOCK : 0) {
if (val) clean();
}
// Returns the logical number of bits.
int size() const {
return _n;
}
// Returns whether the bit at index i is set.
bool test(int i) const {
assert(0 <= i && i < _n);
return (blocks[i >> 6] >> (i & (BITS_PER_BLOCK - 1))) & 1;
}
// Sets the bit at index i to true.
void set(int i) {
assert(0 <= i && i < _n);
blocks[i >> 6] |= uint64_t{1} << (i & (BITS_PER_BLOCK - 1));
}
// Sets all bits to true.
void set() {
std::fill(blocks.begin(), blocks.end(), FULL_BLOCK);
clean();
}
// Sets the bit at index i to false.
void reset(int i) {
assert(0 <= i && i < _n);
blocks[i >> 6] &= ~(uint64_t{1} << (i & (BITS_PER_BLOCK - 1)));
}
// Sets all bits to false.
void reset() {
std::fill(blocks.begin(), blocks.end(), uint64_t{0});
}
// Flips the bit at index i.
void flip(int i) {
assert(0 <= i && i < _n);
blocks[i >> 6] ^= uint64_t{1} << (i & (BITS_PER_BLOCK - 1));
}
// Flips all bits.
void flip() {
for (uint64_t& block : blocks) block = ~block;
clean();
}
// Returns the number of set bits.
int popcount() const {
int res = 0;
for (uint64_t block : blocks) res += __builtin_popcountll(block);
return res;
}
// Returns the index of the least significant set bit, or -1 if no bit is set.
int lowbit() const {
const int m = static_cast<int>(blocks.size());
for (int i = 0; i < m; ++i) {
if (blocks[i] != 0) return (i << 6) + __builtin_ctzll(blocks[i]);
}
return -1;
}
// Returns the index of the most significant set bit, or -1 if no bit is set.
int topbit() const {
for (int i = static_cast<int>(blocks.size()) - 1; i >= 0; --i) {
if (blocks[i] != 0) return (i << 6) + (BITS_PER_BLOCK - 1 - __builtin_clzll(blocks[i]));
}
return -1;
}
// Returns whether at least one bit is set.
bool any() const {
for (uint64_t block : blocks) {
if (block != 0) return true;
}
return false;
}
// Returns whether every logical bit is set.
bool all() const {
if (_n == 0) return true;
const int m = static_cast<int>(blocks.size());
for (int i = 0; i + 1 < m; ++i) {
if (blocks[i] != FULL_BLOCK) return false;
}
return blocks.back() == tail_mask();
}
// Returns whether no bit is set.
bool none() const {
return !any();
}
DynamicBitset& operator&=(const DynamicBitset& other) {
assert(_n == other._n);
const std::size_t m = blocks.size();
for (std::size_t i = 0; i < m; ++i) blocks[i] &= other.blocks[i];
return *this;
}
DynamicBitset& operator|=(const DynamicBitset& other) {
assert(_n == other._n);
const std::size_t m = blocks.size();
for (std::size_t i = 0; i < m; ++i) blocks[i] |= other.blocks[i];
return *this;
}
DynamicBitset& operator^=(const DynamicBitset& other) {
assert(_n == other._n);
const std::size_t m = blocks.size();
for (std::size_t i = 0; i < m; ++i) blocks[i] ^= other.blocks[i];
return *this;
}
DynamicBitset operator~() const {
DynamicBitset res = *this;
res.flip();
return res;
}
friend DynamicBitset operator&(DynamicBitset lhs, const DynamicBitset& rhs) {
lhs &= rhs;
return lhs;
}
friend DynamicBitset operator|(DynamicBitset lhs, const DynamicBitset& rhs) {
lhs |= rhs;
return lhs;
}
friend DynamicBitset operator^(DynamicBitset lhs, const DynamicBitset& rhs) {
lhs ^= rhs;
return lhs;
}
};
} // namespace utilities
} // namespace m1une
#line 13 "graph/dag_reachability.hpp"
namespace m1une {
namespace graph {
struct DagReachability {
std::vector<utilities::DynamicBitset> reachable_vertices;
std::vector<int> topological_order;
int size() const {
return int(reachable_vertices.size());
}
bool reachable(int from, int to) const {
assert(0 <= from && from < size());
assert(0 <= to && to < size());
return reachable_vertices[from].test(to);
}
};
template <class T>
std::optional<DagReachability> dag_reachability(const Graph<T>& g) {
const int n = g.size();
auto order = topological_sort(g);
if (!order) return std::nullopt;
DagReachability result;
result.reachable_vertices.assign(n, utilities::DynamicBitset(n));
result.topological_order = *order;
for (int i = n - 1; i >= 0; i--) {
int v = (*order)[i];
result.reachable_vertices[v].set(v);
for (const auto& e : g[v]) {
if (e.alive) result.reachable_vertices[v] |= result.reachable_vertices[e.to];
}
}
return result;
}
template <class T>
struct DagTransitiveReductionResult {
Graph<T> graph;
std::vector<int> original_edge_ids;
};
template <class T>
std::optional<DagTransitiveReductionResult<T>> dag_transitive_reduction(const Graph<T>& g) {
auto reachability = dag_reachability(g);
if (!reachability) return std::nullopt;
const int n = g.size();
std::vector<int> position(n);
for (int i = 0; i < n; i++) position[reachability->topological_order[i]] = i;
std::vector<char> kept(g.edge_count(), false);
for (int v = 0; v < n; v++) {
std::vector<const Edge<T>*> outgoing;
outgoing.reserve(g[v].size());
for (const auto& e : g[v]) {
if (e.alive) outgoing.push_back(&e);
}
std::stable_sort(outgoing.begin(), outgoing.end(), [&](const auto* lhs, const auto* rhs) {
return position[lhs->to] < position[rhs->to];
});
utilities::DynamicBitset covered(n);
for (const auto* e : outgoing) {
if (covered.test(e->to)) continue;
kept[e->id] = true;
covered |= reachability->reachable_vertices[e->to];
}
}
DagTransitiveReductionResult<T> result;
result.graph = Graph<T>(n);
for (int v = 0; v < n; v++) {
for (const auto& e : g[v]) {
if (!e.alive || !kept[e.id]) continue;
result.graph.add_directed_edge(e.from, e.to, e.cost);
result.original_edge_ids.push_back(e.id);
}
}
return result;
}
} // namespace graph
} // namespace m1une
#line 1 "graph/dag_shortest_path.hpp"
#line 9 "graph/dag_shortest_path.hpp"
#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 10 "graph/dag.hpp"
#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/directed_mst.hpp"
#line 8 "graph/directed_mst.hpp"
#line 10 "graph/directed_mst.hpp"
namespace m1une {
namespace graph {
template <class T>
struct DirectedMinimumSpanningTree {
T cost;
std::vector<int> parent;
std::vector<int> parent_edge;
std::vector<Edge<T>> edges;
int root;
};
namespace internal {
template <class T>
struct DirectedMstEdge {
int from = -1;
int to = -1;
T cost = T(0);
int id = -1;
};
template <class T>
struct DirectedMstHeapPool {
using StoredEdge = DirectedMstEdge<T>;
struct Node {
StoredEdge edge;
T offset = T(0);
int child = -1;
int sibling = -1;
};
struct Heap {
int root = -1;
int size = 0;
};
std::vector<Node> nodes;
explicit DirectedMstHeapPool(int capacity = 0) {
nodes.reserve(capacity);
}
T key(int node) const {
return nodes[node].edge.cost + nodes[node].offset;
}
int meld_roots(int first, int second) {
if (first == -1) return second;
if (second == -1) return first;
if (key(second) < key(first)) std::swap(first, second);
nodes[second].offset -= nodes[first].offset;
nodes[second].sibling = nodes[first].child;
nodes[first].child = second;
return first;
}
void push(Heap& heap, const StoredEdge& edge) {
const int node = int(nodes.size());
nodes.push_back(Node{edge, T(0), -1, -1});
heap.root = meld_roots(heap.root, node);
heap.size++;
}
void meld(Heap& destination, Heap& source) {
destination.root = meld_roots(destination.root, source.root);
destination.size += source.size;
source.root = -1;
source.size = 0;
}
const StoredEdge& top(const Heap& heap) const {
assert(heap.root != -1);
return nodes[heap.root].edge;
}
T top_key(const Heap& heap) const {
assert(heap.root != -1);
return key(heap.root);
}
void add_all(Heap& heap, const T& delta) {
assert(heap.root != -1);
nodes[heap.root].offset += delta;
}
void pop(Heap& heap) {
assert(heap.root != -1 && heap.size > 0);
const int old_root = heap.root;
int child = nodes[old_root].child;
std::vector<int> pairs;
while (child != -1) {
int first = child;
child = nodes[first].sibling;
nodes[first].sibling = -1;
nodes[first].offset += nodes[old_root].offset;
if (child != -1) {
int second = child;
child = nodes[second].sibling;
nodes[second].sibling = -1;
nodes[second].offset += nodes[old_root].offset;
first = meld_roots(first, second);
}
pairs.push_back(first);
}
heap.root = -1;
for (auto it = pairs.rbegin(); it != pairs.rend(); ++it) {
heap.root = meld_roots(*it, heap.root);
}
heap.size--;
}
};
struct DirectedMstDsu {
std::vector<int> parent;
explicit DirectedMstDsu(int n) : parent(n, -1) {}
int leader(int vertex) {
int root = vertex;
while (parent[root] != -1) root = parent[root];
while (vertex != root) {
int next = parent[vertex];
parent[vertex] = root;
vertex = next;
}
return root;
}
};
template <class T>
struct DirectedMstRootlessCost {
int artificial_edges;
T original_cost;
DirectedMstRootlessCost() : artificial_edges(0), original_cost(T(0)) {}
explicit DirectedMstRootlessCost(int zero)
: artificial_edges(zero), original_cost(T(0)) {
assert(zero == 0);
}
DirectedMstRootlessCost(int artificial_edges_, const T& original_cost_)
: artificial_edges(artificial_edges_), original_cost(original_cost_) {}
DirectedMstRootlessCost& operator+=(const DirectedMstRootlessCost& other) {
artificial_edges += other.artificial_edges;
original_cost += other.original_cost;
return *this;
}
DirectedMstRootlessCost& operator-=(const DirectedMstRootlessCost& other) {
artificial_edges -= other.artificial_edges;
original_cost -= other.original_cost;
return *this;
}
friend DirectedMstRootlessCost operator+(
DirectedMstRootlessCost first,
const DirectedMstRootlessCost& second
) {
return first += second;
}
friend DirectedMstRootlessCost operator-(
DirectedMstRootlessCost first,
const DirectedMstRootlessCost& second
) {
return first -= second;
}
friend bool operator<(
const DirectedMstRootlessCost& first,
const DirectedMstRootlessCost& second
) {
if (first.artificial_edges != second.artificial_edges) {
return first.artificial_edges < second.artificial_edges;
}
return first.original_cost < second.original_cost;
}
};
} // namespace internal
// Returns a minimum-cost spanning arborescence rooted at root, or nullopt when
// some vertex is unreachable from the root using active directed edges.
template <class T>
std::optional<DirectedMinimumSpanningTree<T>> directed_mst(
const Graph<T>& graph,
int root
) {
const int n = graph.size();
assert(0 <= root && root < n);
const int maximum_node_count = 2 * n;
int active_edge_count = 0;
#ifndef NDEBUG
std::vector<int> incidence(graph.edge_count(), 0);
#endif
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());
#ifndef NDEBUG
incidence[edge.id]++;
#endif
active_edge_count++;
}
}
#ifndef NDEBUG
for (int count : incidence) {
if (count != 0) assert(count == 1);
}
#endif
using StoredEdge = internal::DirectedMstEdge<T>;
using HeapPool = internal::DirectedMstHeapPool<T>;
HeapPool pool(active_edge_count);
std::vector<typename HeapPool::Heap> heaps(maximum_node_count);
for (int vertex = 0; vertex < n; vertex++) {
for (const Edge<T>& edge : graph[vertex]) {
if (!edge.alive) continue;
pool.push(heaps[edge.to], StoredEdge{edge.from, edge.to, edge.cost, edge.id});
}
}
internal::DirectedMstDsu dsu(maximum_node_count);
std::vector<int> contraction_parent(maximum_node_count, -1);
std::vector<int> visited(maximum_node_count, 0);
std::vector<StoredEdge> selected(maximum_node_count);
int node_count = n;
int visit_token = 1;
visited[root] = 1;
for (int start = 0; start < n; start++) {
if (visited[start] != 0) continue;
visit_token++;
int component = start;
while (visited[component] == 0 || visited[component] == visit_token) {
if (visited[component] == visit_token) {
if (node_count == maximum_node_count) return std::nullopt;
const int contracted = node_count++;
int current = component;
do {
const T reduction = T(0) - pool.top_key(heaps[current]);
pool.add_all(heaps[current], reduction);
pool.meld(heaps[contracted], heaps[current]);
contraction_parent[current] = contracted;
dsu.parent[current] = contracted;
current = dsu.leader(selected[current].from);
} while (current != contracted);
component = contracted;
}
assert(visited[component] == 0);
visited[component] = visit_token;
while (heaps[component].size > 0 &&
dsu.leader(pool.top(heaps[component]).from) == component) {
pool.pop(heaps[component]);
}
if (heaps[component].size == 0) return std::nullopt;
selected[component] = pool.top(heaps[component]);
component = dsu.leader(selected[component].from);
}
}
DirectedMinimumSpanningTree<T> result;
result.cost = T(0);
result.parent.assign(n, -1);
result.parent_edge.assign(n, -1);
result.root = root;
result.parent[root] = root;
std::vector<char> expanded(node_count, false);
std::vector<StoredEdge> chosen(n);
for (int component = node_count - 1; component >= 0; component--) {
if (component == root || expanded[component]) continue;
const StoredEdge& edge = selected[component];
if (edge.id == -1) return std::nullopt;
int vertex = edge.to;
while (vertex != -1 && !expanded[vertex]) {
expanded[vertex] = true;
vertex = contraction_parent[vertex];
}
result.cost += edge.cost;
result.parent[edge.to] = edge.from;
result.parent_edge[edge.to] = edge.id;
chosen[edge.to] = edge;
}
result.edges.reserve(n - 1);
for (int vertex = 0; vertex < n; vertex++) {
if (vertex == root) continue;
if (result.parent[vertex] == -1) return std::nullopt;
const StoredEdge& edge = chosen[vertex];
result.edges.emplace_back(edge.from, edge.to, edge.cost, edge.id, true);
}
return result;
}
// Chooses the root that gives a minimum-cost spanning arborescence.
template <class T>
std::optional<DirectedMinimumSpanningTree<T>> directed_mst(
const Graph<T>& graph
) {
const int n = graph.size();
if (n == 0) return std::nullopt;
using Cost = internal::DirectedMstRootlessCost<T>;
Graph<Cost> augmented(n + 1);
std::vector<int> original_edge_id;
original_edge_id.reserve(graph.edge_count() + n);
#ifndef NDEBUG
std::vector<int> incidence(graph.edge_count(), 0);
#endif
for (int vertex = 0; vertex < n; vertex++) {
for (const Edge<T>& edge : graph[vertex]) {
if (!edge.alive) continue;
#ifndef NDEBUG
assert(0 <= edge.id && edge.id < graph.edge_count());
incidence[edge.id]++;
#endif
augmented.add_directed_edge(
edge.from,
edge.to,
Cost(0, edge.cost)
);
original_edge_id.push_back(edge.id);
}
}
#ifndef NDEBUG
for (int count : incidence) {
if (count != 0) assert(count == 1);
}
#endif
const int artificial_root = n;
for (int vertex = 0; vertex < n; vertex++) {
augmented.add_directed_edge(
artificial_root,
vertex,
Cost(1, T(0))
);
original_edge_id.push_back(-1);
}
auto augmented_result = directed_mst(augmented, artificial_root);
if (!augmented_result || augmented_result->cost.artificial_edges != 1) {
return std::nullopt;
}
DirectedMinimumSpanningTree<T> result;
result.cost = augmented_result->cost.original_cost;
result.parent.assign(n, -1);
result.parent_edge.assign(n, -1);
result.root = -1;
result.edges.reserve(n - 1);
for (int vertex = 0; vertex < n; vertex++) {
int augmented_edge_id = augmented_result->parent_edge[vertex];
assert(0 <= augmented_edge_id &&
augmented_edge_id < int(original_edge_id.size()));
int edge_id = original_edge_id[augmented_edge_id];
if (edge_id == -1) {
assert(result.root == -1);
result.root = vertex;
result.parent[vertex] = vertex;
continue;
}
result.parent[vertex] = augmented_result->parent[vertex];
result.parent_edge[vertex] = edge_id;
result.edges.emplace_back(
result.parent[vertex],
vertex,
augmented_result->edges[vertex].cost.original_cost,
edge_id,
true
);
}
assert(result.root != -1);
return result;
}
} // 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/functional_graph.hpp"
#line 10 "graph/functional_graph.hpp"
namespace m1une {
namespace graph {
struct FunctionalGraph {
int component_count;
std::vector<int> successor;
std::vector<std::vector<int>> predecessors;
std::vector<std::vector<int>> cycles;
std::vector<int> component;
std::vector<int> component_size;
std::vector<int> cycle_entry;
std::vector<int> cycle_position;
std::vector<int> distance_to_cycle;
private:
std::vector<std::vector<int>> _up;
void check_vertex(int vertex) const {
assert(0 <= vertex && vertex < size());
}
int advance_before_cycle(int vertex, int steps) const {
assert(0 <= steps && steps <= distance_to_cycle[vertex]);
int bit = 0;
while (steps > 0) {
if (steps & 1) vertex = _up[bit][vertex];
steps >>= 1;
bit++;
}
return vertex;
}
public:
FunctionalGraph() : component_count(0) {}
explicit FunctionalGraph(const std::vector<int>& successor_) {
build(successor_);
}
void build(const std::vector<int>& successor_) {
successor = successor_;
const int n = size();
for (int to : successor) assert(0 <= to && to < n);
component_count = 0;
predecessors.assign(n, {});
cycles.clear();
component.assign(n, -1);
cycle_entry.assign(n, -1);
cycle_position.assign(n, -1);
distance_to_cycle.assign(n, -1);
std::vector<int> indegree(n, 0);
for (int vertex = 0; vertex < n; vertex++) {
predecessors[successor[vertex]].push_back(vertex);
indegree[successor[vertex]]++;
}
std::queue<int> queue;
std::vector<char> removed(n, false);
for (int vertex = 0; vertex < n; vertex++) {
if (indegree[vertex] == 0) queue.push(vertex);
}
while (!queue.empty()) {
const int vertex = queue.front();
queue.pop();
removed[vertex] = true;
const int to = successor[vertex];
indegree[to]--;
if (indegree[to] == 0) queue.push(to);
}
for (int start = 0; start < n; start++) {
if (removed[start] || component[start] != -1) continue;
const int component_id = int(cycles.size());
std::vector<int> cycle;
int vertex = start;
do {
const int position = int(cycle.size());
cycle.push_back(vertex);
component[vertex] = component_id;
cycle_entry[vertex] = vertex;
cycle_position[vertex] = position;
distance_to_cycle[vertex] = 0;
vertex = successor[vertex];
} while (vertex != start);
cycles.push_back(std::move(cycle));
}
component_count = int(cycles.size());
for (const std::vector<int>& cycle : cycles) {
for (int vertex : cycle) queue.push(vertex);
}
while (!queue.empty()) {
const int vertex = queue.front();
queue.pop();
for (int from : predecessors[vertex]) {
if (component[from] != -1) continue;
component[from] = component[vertex];
cycle_entry[from] = cycle_entry[vertex];
cycle_position[from] = cycle_position[vertex];
distance_to_cycle[from] = distance_to_cycle[vertex] + 1;
queue.push(from);
}
}
component_size.assign(component_count, 0);
for (int component_id : component) component_size[component_id]++;
int log = 1;
while ((std::uint64_t(1) << log) <= std::uint64_t(n)) log++;
_up.assign(log, successor);
for (int bit = 1; bit < log; bit++) {
for (int vertex = 0; vertex < n; vertex++) {
_up[bit][vertex] = _up[bit - 1][_up[bit - 1][vertex]];
}
}
}
int size() const {
return int(successor.size());
}
bool empty() const {
return successor.empty();
}
bool same_component(int first, int second) const {
check_vertex(first);
check_vertex(second);
return component[first] == component[second];
}
bool on_cycle(int vertex) const {
check_vertex(vertex);
return distance_to_cycle[vertex] == 0;
}
int cycle_size(int vertex) const {
check_vertex(vertex);
return int(cycles[component[vertex]].size());
}
int orbit_size(int vertex) const {
check_vertex(vertex);
return distance_to_cycle[vertex] + cycle_size(vertex);
}
int jump(int vertex, std::uint64_t steps) const {
check_vertex(vertex);
const int tail_length = distance_to_cycle[vertex];
if (steps < std::uint64_t(tail_length)) {
return advance_before_cycle(vertex, int(steps));
}
steps -= std::uint64_t(tail_length);
const int entry = cycle_entry[vertex];
const int length = cycle_size(entry);
const int offset = int(steps % std::uint64_t(length));
const int position = (cycle_position[entry] + offset) % length;
return cycles[component[vertex]][position];
}
long long distance(int from, int to) const {
check_vertex(from);
check_vertex(to);
if (!same_component(from, to)) return -1;
if (!on_cycle(to)) {
if (distance_to_cycle[from] < distance_to_cycle[to]) return -1;
const int difference = distance_to_cycle[from] - distance_to_cycle[to];
return advance_before_cycle(from, difference) == to ? difference : -1;
}
const int entry = cycle_entry[from];
const int length = cycle_size(from);
int cycle_distance = cycle_position[to] - cycle_position[entry];
if (cycle_distance < 0) cycle_distance += length;
return static_cast<long long>(distance_to_cycle[from]) + cycle_distance;
}
bool reachable(int from, int to) const {
return distance(from, to) != -1;
}
std::vector<int> path(int from, int to) const {
const long long path_length = distance(from, to);
if (path_length == -1) return {};
std::vector<int> result;
result.reserve(path_length + 1);
for (long long step = 0; step <= path_length; step++) {
result.push_back(from);
from = successor[from];
}
return result;
}
std::vector<int> orbit(int vertex) const {
check_vertex(vertex);
const int length = orbit_size(vertex);
std::vector<int> result;
result.reserve(length);
for (int step = 0; step < length; step++) {
result.push_back(vertex);
vertex = successor[vertex];
}
return result;
}
std::uint64_t visit_count(
int from,
int to,
std::uint64_t step_count
) const {
const long long first_visit = distance(from, to);
if (first_visit == -1 ||
std::uint64_t(first_visit) >= step_count) {
return 0;
}
if (!on_cycle(to)) return 1;
const std::uint64_t remaining =
step_count - 1 - std::uint64_t(first_visit);
return 1 + remaining / std::uint64_t(cycle_size(to));
}
long long first_meeting_time(int first, int second) const {
check_vertex(first);
check_vertex(second);
if (!same_component(first, second)) return -1;
if (first == second) return 0;
const int first_depth = distance_to_cycle[first];
const int second_depth = distance_to_cycle[second];
if (first_depth == second_depth &&
cycle_entry[first] == cycle_entry[second]) {
int elapsed = 0;
for (int bit = int(_up.size()) - 1; bit >= 0; bit--) {
const int steps = 1 << bit;
if (first_depth - elapsed < steps) continue;
const int next_first = _up[bit][first];
const int next_second = _up[bit][second];
if (next_first == next_second) continue;
first = next_first;
second = next_second;
elapsed += steps;
}
return elapsed + 1;
}
const int length = cycle_size(first);
int first_phase =
cycle_position[first] - first_depth % length;
int second_phase =
cycle_position[second] - second_depth % length;
if (first_phase < 0) first_phase += length;
if (second_phase < 0) second_phase += length;
if (first_phase != second_phase) return -1;
return std::max(first_depth, second_depth);
}
int first_meeting_vertex(int first, int second) const {
const long long time = first_meeting_time(first, second);
if (time == -1) return -1;
return jump(first, std::uint64_t(time));
}
};
} // namespace graph
} // namespace m1une
#line 1 "graph/incremental_scc.hpp"
#line 9 "graph/incremental_scc.hpp"
#line 11 "graph/incremental_scc.hpp"
namespace m1une {
namespace graph {
namespace incremental_scc_detail {
struct EdgeEvent {
int id;
int from;
int to;
};
inline std::vector<int> component_ids(
int vertex_count,
const std::vector<EdgeEvent>& edges,
int time
) {
std::vector<int> begin(vertex_count + 1, 0);
std::vector<int> reverse_begin(vertex_count + 1, 0);
int edge_count = 0;
for (const EdgeEvent& edge : edges) {
if (edge.id >= time) continue;
begin[edge.from + 1]++;
reverse_begin[edge.to + 1]++;
edge_count++;
}
for (int vertex = 0; vertex < vertex_count; vertex++) {
begin[vertex + 1] += begin[vertex];
reverse_begin[vertex + 1] += reverse_begin[vertex];
}
std::vector<int> adjacency(edge_count);
std::vector<int> reverse_adjacency(edge_count);
std::vector<int> cursor = begin;
std::vector<int> reverse_cursor = reverse_begin;
for (const EdgeEvent& edge : edges) {
if (edge.id >= time) continue;
adjacency[cursor[edge.from]++] = edge.to;
reverse_adjacency[reverse_cursor[edge.to]++] = edge.from;
}
std::vector<int>().swap(cursor);
std::vector<int>().swap(reverse_cursor);
std::vector<char> visited(vertex_count, false);
std::vector<int> next_position(begin.begin(), begin.end() - 1);
std::vector<int> order;
order.reserve(vertex_count);
std::vector<int> stack;
for (int start = 0; start < vertex_count; start++) {
if (visited[start]) continue;
visited[start] = true;
stack.push_back(start);
while (!stack.empty()) {
const int vertex = stack.back();
int& position = next_position[vertex];
if (position < begin[vertex + 1]) {
const int to = adjacency[position++];
if (!visited[to]) {
visited[to] = true;
stack.push_back(to);
}
} else {
order.push_back(vertex);
stack.pop_back();
}
}
}
std::vector<int> component(vertex_count, -1);
int component_count = 0;
for (auto iterator = order.rbegin(); iterator != order.rend(); ++iterator) {
const int start = *iterator;
if (component[start] != -1) continue;
component[start] = component_count;
stack.push_back(start);
while (!stack.empty()) {
const int vertex = stack.back();
stack.pop_back();
for (int position = reverse_begin[vertex];
position < reverse_begin[vertex + 1]; position++) {
const int to = reverse_adjacency[position];
if (component[to] != -1) continue;
component[to] = component_count;
stack.push_back(to);
}
}
component_count++;
}
return component;
}
} // namespace incremental_scc_detail
// For every directed edge e, returns the first time t after e is inserted such
// that its endpoints are in the same SCC. At time t, edges with IDs less than
// t have been inserted. edge_count() + 1 means this never happens.
template <class T>
std::vector<int> incremental_scc(const Graph<T>& graph) {
using incremental_scc_detail::EdgeEvent;
using incremental_scc_detail::component_ids;
const int vertex_count = graph.size();
const int edge_count = graph.edge_count();
const int never = edge_count + 1;
std::vector<int> merge_time(edge_count, never);
if (edge_count == 0) return merge_time;
std::vector<EdgeEvent> edges_by_id(edge_count);
std::vector<char> initialized(edge_count, false);
for (int vertex = 0; vertex < vertex_count; vertex++) {
for (const Edge<T>& edge : graph[vertex]) {
assert(0 <= edge.id && edge.id < edge_count);
assert(!initialized[edge.id]);
if (initialized[edge.id]) continue;
initialized[edge.id] = true;
edges_by_id[edge.id] = EdgeEvent{edge.id, edge.from, edge.to};
}
}
std::vector<EdgeEvent> events;
events.reserve(edge_count);
for (int edge_id = 0; edge_id < edge_count; edge_id++) {
assert(initialized[edge_id]);
if (graph.is_edge_alive(edge_id)) {
events.push_back(edges_by_id[edge_id]);
}
}
std::vector<EdgeEvent>().swap(edges_by_id);
std::vector<char>().swap(initialized);
std::vector<int> new_index(vertex_count, -1);
auto divide = [&](
auto&& self,
std::vector<EdgeEvent> current,
int left,
int right
) -> void {
if (current.empty() || right == left + 1) return;
const int middle = left + (right - left) / 2;
std::vector<int> touched;
touched.reserve(std::min(
std::size_t(vertex_count),
current.size() * 2
));
int compressed_count = 0;
for (const EdgeEvent& edge : current) {
if (new_index[edge.from] == -1) {
new_index[edge.from] = compressed_count++;
touched.push_back(edge.from);
}
if (new_index[edge.to] == -1) {
new_index[edge.to] = compressed_count++;
touched.push_back(edge.to);
}
}
for (EdgeEvent& edge : current) {
edge.from = new_index[edge.from];
edge.to = new_index[edge.to];
}
for (int vertex : touched) new_index[vertex] = -1;
std::vector<EdgeEvent> earlier;
std::vector<EdgeEvent> later;
earlier.reserve(current.size() / 2);
later.reserve(current.size() / 2);
{
std::vector<int> component =
component_ids(compressed_count, current, middle);
for (const EdgeEvent& edge : current) {
const int from_component = component[edge.from];
const int to_component = component[edge.to];
if (edge.id < middle &&
from_component == to_component) {
merge_time[edge.id] =
std::min(merge_time[edge.id], middle);
earlier.push_back(edge);
} else {
later.push_back(EdgeEvent{
edge.id,
from_component,
to_component
});
}
}
}
std::vector<EdgeEvent>().swap(current);
self(self, std::move(earlier), left, middle);
self(self, std::move(later), middle, right);
};
divide(divide, std::move(events), 0, edge_count + 1);
return merge_time;
}
} // 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 5 "math/matrix/linear_algebra.hpp"
#include <type_traits>
#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/scc.hpp"
#line 9 "graph/scc.hpp"
#line 11 "graph/scc.hpp"
namespace m1une {
namespace graph {
struct SccResult {
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>
Graph<int> dag(const Graph<T>& g) const {
std::vector<std::pair<int, int>> edges;
for (int v = 0; v < g.size(); v++) {
for (const auto& e : g[v]) {
if (!e.alive) continue;
int a = comp[e.from], b = comp[e.to];
if (a != b) edges.emplace_back(a, b);
}
}
std::sort(edges.begin(), edges.end());
edges.erase(std::unique(edges.begin(), edges.end()), edges.end());
Graph<int> result(count);
for (auto [a, b] : edges) result.add_directed_edge(a, b);
return result;
}
};
template <class T>
SccResult strongly_connected_components(const Graph<T>& g) {
const int n = g.size();
std::vector<std::vector<int>> reverse_graph(n);
for (int vertex = 0; vertex < n; vertex++) {
for (const auto& edge : g[vertex]) {
if (edge.alive) reverse_graph[edge.to].push_back(vertex);
}
}
std::vector<char> seen(n, false);
std::vector<int> order;
order.reserve(n);
std::vector<std::pair<int, std::size_t>> dfs_stack;
for (int start = 0; start < n; start++) {
if (seen[start]) continue;
seen[start] = true;
dfs_stack.emplace_back(start, 0);
while (!dfs_stack.empty()) {
int vertex = dfs_stack.back().first;
std::size_t& edge_index = dfs_stack.back().second;
while (edge_index < g[vertex].size() &&
!g[vertex][edge_index].alive) {
edge_index++;
}
if (edge_index == g[vertex].size()) {
order.push_back(vertex);
dfs_stack.pop_back();
continue;
}
const int to = g[vertex][edge_index++].to;
if (!seen[to]) {
seen[to] = true;
dfs_stack.emplace_back(to, 0);
}
}
}
std::vector<int> comp(n, -1);
std::vector<std::vector<int>> groups;
std::vector<int> stack;
for (auto iterator = order.rbegin(); iterator != order.rend(); ++iterator) {
const int start = *iterator;
if (comp[start] != -1) continue;
const int component = int(groups.size());
groups.emplace_back();
comp[start] = component;
stack.push_back(start);
while (!stack.empty()) {
const int vertex = stack.back();
stack.pop_back();
groups.back().push_back(vertex);
for (int to : reverse_graph[vertex]) {
if (comp[to] != -1) continue;
comp[to] = component;
stack.push_back(to);
}
}
}
return SccResult{int(groups.size()), std::move(comp), std::move(groups)};
}
} // 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/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 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/two_sat.hpp"
#line 9 "graph/two_sat.hpp"
namespace m1une {
namespace graph {
// A 2-SAT solver using iterative strongly connected components.
struct TwoSat {
private:
struct Csr {
std::vector<int> start;
std::vector<int> to;
};
int _n;
std::vector<std::pair<int, int>> _edges;
bool _solved;
bool _satisfiable;
std::vector<bool> _answer;
int node(int variable, bool value) const {
assert(0 <= variable && variable < _n);
return 2 * variable + int(value);
}
void add_edge(int from, int to) {
_edges.emplace_back(from, to);
_solved = false;
_answer.clear();
}
Csr build_csr(bool reverse) const {
int vertices = 2 * _n;
Csr graph;
graph.start.assign(vertices + 1, 0);
graph.to.resize(_edges.size());
for (auto [from, to] : _edges) {
int source = reverse ? to : from;
graph.start[source + 1]++;
}
for (int v = 0; v < vertices; v++) {
graph.start[v + 1] += graph.start[v];
}
std::vector<int> cursor = graph.start;
for (auto [from, to] : _edges) {
int source = reverse ? to : from;
int target = reverse ? from : to;
graph.to[cursor[source]++] = target;
}
return graph;
}
public:
TwoSat() : TwoSat(0) {}
explicit TwoSat(int n)
: _n(n), _solved(false), _satisfiable(false) {
assert(0 <= n);
assert(n <= std::numeric_limits<int>::max() / 2);
}
int size() const {
return _n;
}
bool empty() const {
return _n == 0;
}
// Reserves space for approximately `clause_count` two-literal clauses.
void reserve(std::size_t clause_count) {
assert(clause_count <= std::size_t(std::numeric_limits<int>::max()) / 2);
_edges.reserve(2 * clause_count);
}
// Adds (variable i == f) OR (variable j == g).
void add_clause(int i, bool f, int j, bool g) {
int a = node(i, f);
int b = node(j, g);
add_edge(a ^ 1, b);
add_edge(b ^ 1, a);
}
// Adds (variable i == f) => (variable j == g).
void add_implication(int i, bool f, int j, bool g) {
add_clause(i, !f, j, g);
}
// Forces variable i to equal value.
void set_value(int i, bool value) {
add_clause(i, value, i, value);
}
// Forces variables i and j to have equal values.
void add_equal(int i, int j) {
add_clause(i, false, j, true);
add_clause(i, true, j, false);
}
// Forces variables i and j to have different values.
void add_not_equal(int i, int j) {
add_clause(i, true, j, true);
add_clause(i, false, j, false);
}
bool satisfiable() {
if (_solved) return _satisfiable;
assert(_edges.size() <= std::size_t(std::numeric_limits<int>::max()));
int vertices = 2 * _n;
Csr graph = build_csr(false);
Csr reverse_graph = build_csr(true);
std::vector<char> seen(vertices, false);
std::vector<int> order;
order.reserve(vertices);
std::vector<std::pair<int, int>> stack;
stack.reserve(vertices);
for (int start = 0; start < vertices; start++) {
if (seen[start]) continue;
seen[start] = true;
stack.emplace_back(start, graph.start[start]);
while (!stack.empty()) {
int v = stack.back().first;
int& edge = stack.back().second;
if (edge == graph.start[v + 1]) {
order.push_back(v);
stack.pop_back();
continue;
}
int to = graph.to[edge++];
if (!seen[to]) {
seen[to] = true;
stack.emplace_back(to, graph.start[to]);
}
}
}
std::vector<int> component(vertices, -1);
std::vector<int> vertices_stack;
vertices_stack.reserve(vertices);
int component_count = 0;
for (int index = vertices - 1; index >= 0; index--) {
int start = order[index];
if (component[start] != -1) continue;
component[start] = component_count;
vertices_stack.push_back(start);
while (!vertices_stack.empty()) {
int v = vertices_stack.back();
vertices_stack.pop_back();
for (int edge = reverse_graph.start[v];
edge < reverse_graph.start[v + 1];
edge++) {
int to = reverse_graph.to[edge];
if (component[to] == -1) {
component[to] = component_count;
vertices_stack.push_back(to);
}
}
}
component_count++;
}
_answer.assign(_n, false);
_satisfiable = true;
for (int i = 0; i < _n; i++) {
if (component[2 * i] == component[2 * i + 1]) {
_satisfiable = false;
_answer.clear();
break;
}
_answer[i] = component[2 * i] < component[2 * i + 1];
}
_solved = true;
return _satisfiable;
}
const std::vector<bool>& answer() const {
assert(_solved && _satisfiable);
return _answer;
}
bool value(int variable) const {
assert(_solved && _satisfiable);
assert(0 <= variable && variable < _n);
return _answer[variable];
}
};
} // namespace graph
} // namespace m1une
#line 16 "graph/directed.hpp"