bipartite_decomposition.hpp¶
二分图最大匹配、最小点覆盖、最大独立集与最小边覆盖；另提供 DAG 原边路径覆盖，以及位集闭包上的最大反链和最小链覆盖。
Complexity: Time: \(O(V + E \sqrt{V})\) for bipartite decomposition / DAG paths. DAG order: \(O(V + (V + E) \lceil V / 64\rceil + R \sqrt{V})\), \(R\) reachable ordered pairs. Space: \(O(V + E)\); DAG order adds \(O(V \lceil V / 64\rceil + R)\).
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: O(V + E sqrt(V)) for bipartite decomposition / DAG paths.
/// DAG order: O(V + (V+E) ceil(V/64) + R sqrt(V)), R reachable ordered pairs.
/// Space: O(V+E); DAG order adds O(V ceil(V/64) + R).
#include "noya/ds/dynamic_bitset.hpp"
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <span>
#include <utility>
#include <vector>
namespace noya {
struct bipartite_decomposition_result {
std::vector<int> ml;
std::vector<int> mr;
std::vector<int> cvl;
std::vector<int> cvr;
std::vector<int> isl;
std::vector<int> isr;
bool ok = false;
std::vector<std::pair<int, int>> ec;
int matching_size() const {
return int(
std::count_if(ml.begin(), ml.end(), [](int val) { return val != -1; }));
}
};
namespace bipartite_decomposition_internal {
struct adjacency {
std::vector<std::size_t> off;
std::vector<int> to;
adjacency(int nl, int nr, const std::vector<std::pair<int, int>> &es)
: off(std::size_t(nl) + 1), to(es.size()) {
for (auto [l, r] : es) {
assert(0 <= l && l < nl && 0 <= r && r < nr);
++off[l + 1];
}
for (int l = 0; l < nl; ++l) off[l + 1] += off[l];
auto cur = off;
for (auto [l, r] : es) to[cur[l]++] = r;
}
explicit adjacency(const std::vector<dynamic_bitset> &rows)
: off(rows.size() + 1) {
for (std::size_t u = 0; u < rows.size(); ++u)
off[u + 1] = off[u] + rows[u].count();
to.reserve(off.back());
for (const auto &row : rows)
for (std::size_t v = row.find_first(); v < row.size(); v = row.find_next(v + 1))
to.push_back(int(v));
}
int size() const { return int(off.size() - 1); }
std::span<const int> operator[](int u) const {
return std::span(to).subspan(off[u], off[u + 1] - off[u]);
}
};
// Return final alternating-reachability distances from unmatched left vertices.
// Per-phase cursors scan each edge at most once during layered augmentation.
inline std::vector<int> match(const adjacency &g, int nr, std::vector<int> &ml,
std::vector<int> &mr) {
int nl = g.size();
ml.assign(nl, -1);
mr.assign(nr, -1);
std::vector<int> active;
for (int l = 0; l < nl; ++l) {
if (g[l].empty()) continue;
active.push_back(l);
for (int r : g[l]) {
if (mr[r] == -1) {
ml[l] = r;
mr[r] = l;
break;
}
}
}
std::vector<int> dis(nl), q;
std::vector<std::size_t> cur;
q.reserve(active.size());
while (true) {
// Isolates keep distance zero: they are unmatched and trivially reachable.
// Only O(E) active vertices need resetting/scanning in each matching phase.
for (int l : active) dis[l] = -1;
q.clear();
for (int l : active) {
if (ml[l] == -1) {
dis[l] = 0;
q.push_back(l);
}
}
int limit = -1;
for (std::size_t head = 0; head < q.size(); ++head) {
int l = q[head];
if (limit != -1 && dis[l] >= limit) continue;
for (int r : g[l]) {
int next = mr[r];
if (next == -1) {
limit = dis[l] + 1;
} else if (dis[next] == -1) {
dis[next] = dis[l] + 1;
q.push_back(next);
}
}
}
if (limit == -1) return dis;
cur.resize(nl);
for (int l : active) cur[l] = 0;
auto dfs = [&](auto &self, int l) -> bool {
auto adj = g[l];
while (cur[l] < adj.size()) {
int r = adj[cur[l]++], next = mr[r];
if (next == -1 ? dis[l] + 1 == limit
: dis[next] == dis[l] + 1 && dis[next] < limit && self(self, next)) {
ml[l] = r;
mr[r] = l;
return true;
}
}
dis[l] = -1;
return false;
};
for (int l : active) if (ml[l] == -1) dfs(dfs, l);
}
}
inline std::vector<int> topology(const adjacency &g) {
int n = g.size();
std::vector<std::size_t> deg(n);
std::vector<int> top;
top.reserve(n);
for (int u = 0; u < n; ++u)
for (int v : g[u]) ++deg[v];
for (int u = 0; u < n; ++u)
if (deg[u] == 0) top.push_back(u);
for (std::size_t at = 0; at < top.size(); ++at)
for (int v : g[top[at]])
if (--deg[v] == 0) top.push_back(v);
assert(int(top.size()) == n);
return top;
}
inline std::vector<std::vector<int>> chains(const std::vector<int> &ml,
const std::vector<int> &mr) {
std::vector<std::vector<int>> res;
for (int s = 0; s < int(mr.size()); ++s) {
if (mr[s] != -1) continue;
res.emplace_back();
for (int u = s; u != -1; u = ml[u]) res.back().push_back(u);
}
return res;
}
} // namespace bipartite_decomposition_internal
/// @brief Maximum matching, Konig minimum vertex cover and its complementary
/// maximum independent set. Also returns a minimum edge cover iff no vertex is
/// isolated (ok=true); otherwise ec is empty. Duplicate edges are allowed.
/// Vertex sets are sorted; matching/edge-cover witnesses need not be canonical.
inline bipartite_decomposition_result
bipartite_decomposition(int nl, int nr,
const std::vector<std::pair<int, int>> &es) {
assert(nl >= 0 && nr >= 0);
using namespace bipartite_decomposition_internal;
adjacency g(nl, nr, es);
bipartite_decomposition_result res;
auto dis = match(g, nr, res.ml, res.mr);
for (int l = 0; l < nl; ++l)
(dis[l] >= 0 ? res.isl : res.cvl).push_back(l);
// At maximum matching every reachable right vertex is matched, and its
// matched left partner is reachable too. No second graph traversal is needed.
for (int r = 0; r < nr; ++r)
(res.mr[r] >= 0 && dis[res.mr[r]] >= 0 ? res.cvr : res.isr).push_back(r);
std::vector<int> first(nr, -1);
for (auto [l, r] : es)
if (first[r] == -1) first[r] = l;
res.ok = true;
for (int l = 0; l < nl; ++l) res.ok &= !g[l].empty();
for (int r = 0; r < nr; ++r) res.ok &= first[r] != -1;
if (res.ok) {
res.ec.reserve(std::size_t(nl) + nr - res.cvl.size() - res.cvr.size());
for (int l = 0; l < nl; ++l)
if (res.ml[l] != -1) res.ec.emplace_back(l, res.ml[l]);
for (int l = 0; l < nl; ++l)
if (res.ml[l] == -1) res.ec.emplace_back(l, g[l][0]);
for (int r = 0; r < nr; ++r)
if (res.mr[r] == -1) res.ec.emplace_back(first[r], r);
}
return res;
}
/// @brief Minimum vertex-disjoint directed path cover of a DAG. Consecutive
/// vertices in a returned path are connected by an original edge. Duplicate
/// edges are allowed; the input must be acyclic.
inline std::vector<std::vector<int>>
minimum_dag_path_cover(int n, const std::vector<std::pair<int, int>> &es) {
assert(n >= 0);
using namespace bipartite_decomposition_internal;
adjacency g(n, n, es);
topology(g);
std::vector<int> ml, mr;
match(g, n, ml, mr);
return chains(ml, mr);
}
struct dag_order_decomposition_result {
std::vector<int> ac;
std::vector<std::vector<int>> cc;
};
/// @brief Maximum reachability antichain and minimum chain cover of a DAG.
/// Chain neighbors need only be reachable, unlike minimum_dag_path_cover.
/// Bitset closure takes O((n+E) ceil(n/64)); matching then takes O(n+R sqrt(n))
/// for R reachable ordered pairs. The input must be acyclic.
inline dag_order_decomposition_result
dag_order_decomposition(int n, const std::vector<std::pair<int, int>> &es) {
assert(n >= 0);
using namespace bipartite_decomposition_internal;
adjacency g(n, n, es);
auto top = topology(g);
std::vector<dynamic_bitset> reach(n, dynamic_bitset(n));
for (auto it = top.rbegin(); it != top.rend(); ++it) {
for (int v : g[*it]) {
if (reach[*it][v]) continue;
reach[*it].set(v);
reach[*it] |= reach[v];
}
}
adjacency comparable(reach);
std::vector<int> ml, mr;
auto dis = match(comparable, n, ml, mr);
dag_order_decomposition_result res;
for (int u = 0; u < n; ++u)
if (dis[u] >= 0 && (mr[u] == -1 || dis[mr[u]] < 0)) res.ac.push_back(u);
res.cc = chains(ml, mr);
return res;
}
} // namespace noya
// NOYA_BIPARTITE_DECOMPOSITION_HPP
#ifndef NOYA_BIPARTITE_DECOMPOSITION_HPP
#define NOYA_BIPARTITE_DECOMPOSITION_HPP 1
/// @complexity Time: O(V + E sqrt(V)) for bipartite decomposition / DAG paths.
/// DAG order: O(V + (V+E) ceil(V/64) + R sqrt(V)), R reachable ordered pairs.
/// Space: O(V+E); DAG order adds O(V ceil(V/64) + R).
#include "noya/ds/dynamic_bitset.hpp"
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <span>
#include <utility>
#include <vector>
namespace noya {
struct bipartite_decomposition_result {
std::vector<int> ml;
std::vector<int> mr;
std::vector<int> cvl;
std::vector<int> cvr;
std::vector<int> isl;
std::vector<int> isr;
bool ok = false;
std::vector<std::pair<int, int>> ec;
int matching_size() const {
return int(
std::count_if(ml.begin(), ml.end(), [](int val) { return val != -1; }));
}
};
namespace bipartite_decomposition_internal {
struct adjacency {
std::vector<std::size_t> off;
std::vector<int> to;
adjacency(int nl, int nr, const std::vector<std::pair<int, int>> &es)
: off(std::size_t(nl) + 1), to(es.size()) {
for (auto [l, r] : es) {
assert(0 <= l && l < nl && 0 <= r && r < nr);
++off[l + 1];
}
for (int l = 0; l < nl; ++l) off[l + 1] += off[l];
auto cur = off;
for (auto [l, r] : es) to[cur[l]++] = r;
}
explicit adjacency(const std::vector<dynamic_bitset> &rows)
: off(rows.size() + 1) {
for (std::size_t u = 0; u < rows.size(); ++u)
off[u + 1] = off[u] + rows[u].count();
to.reserve(off.back());
for (const auto &row : rows)
for (std::size_t v = row.find_first(); v < row.size(); v = row.find_next(v + 1))
to.push_back(int(v));
}
int size() const { return int(off.size() - 1); }
std::span<const int> operator[](int u) const {
return std::span(to).subspan(off[u], off[u + 1] - off[u]);
}
};
// Return final alternating-reachability distances from unmatched left vertices.
// Per-phase cursors scan each edge at most once during layered augmentation.
inline std::vector<int> match(const adjacency &g, int nr, std::vector<int> &ml,
std::vector<int> &mr) {
int nl = g.size();
ml.assign(nl, -1);
mr.assign(nr, -1);
std::vector<int> active;
for (int l = 0; l < nl; ++l) {
if (g[l].empty()) continue;
active.push_back(l);
for (int r : g[l]) {
if (mr[r] == -1) {
ml[l] = r;
mr[r] = l;
break;
}
}
}
std::vector<int> dis(nl), q;
std::vector<std::size_t> cur;
q.reserve(active.size());
while (true) {
// Isolates keep distance zero: they are unmatched and trivially reachable.
// Only O(E) active vertices need resetting/scanning in each matching phase.
for (int l : active) dis[l] = -1;
q.clear();
for (int l : active) {
if (ml[l] == -1) {
dis[l] = 0;
q.push_back(l);
}
}
int limit = -1;
for (std::size_t head = 0; head < q.size(); ++head) {
int l = q[head];
if (limit != -1 && dis[l] >= limit) continue;
for (int r : g[l]) {
int next = mr[r];
if (next == -1) {
limit = dis[l] + 1;
} else if (dis[next] == -1) {
dis[next] = dis[l] + 1;
q.push_back(next);
}
}
}
if (limit == -1) return dis;
cur.resize(nl);
for (int l : active) cur[l] = 0;
auto dfs = [&](auto &self, int l) -> bool {
auto adj = g[l];
while (cur[l] < adj.size()) {
int r = adj[cur[l]++], next = mr[r];
if (next == -1 ? dis[l] + 1 == limit
: dis[next] == dis[l] + 1 && dis[next] < limit && self(self, next)) {
ml[l] = r;
mr[r] = l;
return true;
}
}
dis[l] = -1;
return false;
};
for (int l : active) if (ml[l] == -1) dfs(dfs, l);
}
}
inline std::vector<int> topology(const adjacency &g) {
int n = g.size();
std::vector<std::size_t> deg(n);
std::vector<int> top;
top.reserve(n);
for (int u = 0; u < n; ++u)
for (int v : g[u]) ++deg[v];
for (int u = 0; u < n; ++u)
if (deg[u] == 0) top.push_back(u);
for (std::size_t at = 0; at < top.size(); ++at)
for (int v : g[top[at]])
if (--deg[v] == 0) top.push_back(v);
assert(int(top.size()) == n);
return top;
}
inline std::vector<std::vector<int>> chains(const std::vector<int> &ml,
const std::vector<int> &mr) {
std::vector<std::vector<int>> res;
for (int s = 0; s < int(mr.size()); ++s) {
if (mr[s] != -1) continue;
res.emplace_back();
for (int u = s; u != -1; u = ml[u]) res.back().push_back(u);
}
return res;
}
} // namespace bipartite_decomposition_internal
/// @brief Maximum matching, Konig minimum vertex cover and its complementary
/// maximum independent set. Also returns a minimum edge cover iff no vertex is
/// isolated (ok=true); otherwise ec is empty. Duplicate edges are allowed.
/// Vertex sets are sorted; matching/edge-cover witnesses need not be canonical.
inline bipartite_decomposition_result
bipartite_decomposition(int nl, int nr,
const std::vector<std::pair<int, int>> &es) {
assert(nl >= 0 && nr >= 0);
using namespace bipartite_decomposition_internal;
adjacency g(nl, nr, es);
bipartite_decomposition_result res;
auto dis = match(g, nr, res.ml, res.mr);
for (int l = 0; l < nl; ++l)
(dis[l] >= 0 ? res.isl : res.cvl).push_back(l);
// At maximum matching every reachable right vertex is matched, and its
// matched left partner is reachable too. No second graph traversal is needed.
for (int r = 0; r < nr; ++r)
(res.mr[r] >= 0 && dis[res.mr[r]] >= 0 ? res.cvr : res.isr).push_back(r);
std::vector<int> first(nr, -1);
for (auto [l, r] : es)
if (first[r] == -1) first[r] = l;
res.ok = true;
for (int l = 0; l < nl; ++l) res.ok &= !g[l].empty();
for (int r = 0; r < nr; ++r) res.ok &= first[r] != -1;
if (res.ok) {
res.ec.reserve(std::size_t(nl) + nr - res.cvl.size() - res.cvr.size());
for (int l = 0; l < nl; ++l)
if (res.ml[l] != -1) res.ec.emplace_back(l, res.ml[l]);
for (int l = 0; l < nl; ++l)
if (res.ml[l] == -1) res.ec.emplace_back(l, g[l][0]);
for (int r = 0; r < nr; ++r)
if (res.mr[r] == -1) res.ec.emplace_back(first[r], r);
}
return res;
}
/// @brief Minimum vertex-disjoint directed path cover of a DAG. Consecutive
/// vertices in a returned path are connected by an original edge. Duplicate
/// edges are allowed; the input must be acyclic.
inline std::vector<std::vector<int>>
minimum_dag_path_cover(int n, const std::vector<std::pair<int, int>> &es) {
assert(n >= 0);
using namespace bipartite_decomposition_internal;
adjacency g(n, n, es);
topology(g);
std::vector<int> ml, mr;
match(g, n, ml, mr);
return chains(ml, mr);
}
struct dag_order_decomposition_result {
std::vector<int> ac;
std::vector<std::vector<int>> cc;
};
/// @brief Maximum reachability antichain and minimum chain cover of a DAG.
/// Chain neighbors need only be reachable, unlike minimum_dag_path_cover.
/// Bitset closure takes O((n+E) ceil(n/64)); matching then takes O(n+R sqrt(n))
/// for R reachable ordered pairs. The input must be acyclic.
inline dag_order_decomposition_result
dag_order_decomposition(int n, const std::vector<std::pair<int, int>> &es) {
assert(n >= 0);
using namespace bipartite_decomposition_internal;
adjacency g(n, n, es);
auto top = topology(g);
std::vector<dynamic_bitset> reach(n, dynamic_bitset(n));
for (auto it = top.rbegin(); it != top.rend(); ++it) {
for (int v : g[*it]) {
if (reach[*it][v]) continue;
reach[*it].set(v);
reach[*it] |= reach[v];
}
}
adjacency comparable(reach);
std::vector<int> ml, mr;
auto dis = match(comparable, n, ml, mr);
dag_order_decomposition_result res;
for (int u = 0; u < n; ++u)
if (dis[u] >= 0 && (mr[u] == -1 || dis[mr[u]] < 0)) res.ac.push_back(u);
res.cc = chains(ml, mr);
return res;
}
} // namespace noya
#endif // NOYA_BIPARTITE_DECOMPOSITION_HPP
/// @complexity Time: O(V + E sqrt(V)) for bipartite decomposition / DAG paths.
/// DAG order: O(V + (V+E) ceil(V/64) + R sqrt(V)), R reachable ordered pairs.
/// Space: O(V+E); DAG order adds O(V ceil(V/64) + R).
/// @complexity Time: O(n / 64) for whole-bitset operations; O(1) bit access.
/// Space: O(n / 64).
#include <algorithm>
#include <bit>
#include <cassert>
#include <cstddef>
#include <cstdint>
#include <vector>
namespace noya {
/// @brief Resizable packed bitset with bitwise operations, shifts, population
/// count, and efficient iteration over set bits.
struct dynamic_bitset {
using word_type = std::uint64_t;
static constexpr std::size_t W = 64;
std::size_t sz_ = 0;
std::vector<word_type> a;
dynamic_bitset() = default;
explicit dynamic_bitset(std::size_t siz, bool val = false)
: sz_(siz), a(word_count(siz), val ? ~word_type{} : 0) {
trim();
}
std::size_t size() const { return sz_; }
bool empty() const { return sz_ == 0; }
bool test(std::size_t pos) const {
assert(pos < sz_);
return (a[pos / W] >> (pos % W)) & 1;
}
bool operator[](std::size_t pos) const { return test(pos); }
dynamic_bitset &set(std::size_t pos, bool val = true) {
assert(pos < sz_);
word_type msk = word_type(1) << (pos % W);
if (val) {
a[pos / W] |= msk;
} else {
a[pos / W] &= ~msk;
}
return *this;
}
dynamic_bitset &reset(std::size_t pos) { return set(pos, false); }
dynamic_bitset &flip(std::size_t pos) {
assert(pos < sz_);
a[pos / W] ^= word_type(1) << (pos % W);
return *this;
}
dynamic_bitset &set() {
std::fill(a.begin(), a.end(), ~word_type{});
trim();
return *this;
}
dynamic_bitset &reset() {
std::fill(a.begin(), a.end(), word_type{});
return *this;
}
dynamic_bitset &flip() {
for (word_type &w : a) {
w = ~w;
}
trim();
return *this;
}
std::size_t count() const {
std::size_t res = 0;
for (word_type w : a) {
res += std::popcount(w);
}
return res;
}
bool any() const {
return std::any_of(a.begin(), a.end(), [](word_type w) { return w != 0; });
}
bool none() const { return !any(); }
/// @brief Return the first set position at least position, or size() if no
/// such position exists.
std::size_t find_next(std::size_t pos) const {
if (pos >= sz_) {
return sz_;
}
std::size_t idx = pos / W;
word_type w = a[idx] & (~word_type{} << (pos % W));
if (w != 0) {
return std::min(sz_, idx * W + std::size_t(std::countr_zero(w)));
}
for (idx++; idx < a.size(); idx++) {
if (a[idx] != 0) {
return std::min(sz_, idx * W + std::size_t(std::countr_zero(a[idx])));
}
}
return sz_;
}
std::size_t find_first() const { return find_next(0); }
dynamic_bitset &operator&=(const dynamic_bitset &rhs) {
check_same_size(rhs);
for (std::size_t i = 0; i < a.size(); i++) {
a[i] &= rhs.a[i];
}
return *this;
}
dynamic_bitset &operator|=(const dynamic_bitset &rhs) {
check_same_size(rhs);
for (std::size_t i = 0; i < a.size(); i++) {
a[i] |= rhs.a[i];
}
return *this;
}
dynamic_bitset &operator^=(const dynamic_bitset &rhs) {
check_same_size(rhs);
for (std::size_t i = 0; i < a.size(); i++) {
a[i] ^= rhs.a[i];
}
return *this;
}
dynamic_bitset &operator<<=(std::size_t k) {
if (k == 0)
return *this;
if (k >= sz_)
return reset();
std::size_t blk = k / W;
int rem = int(k % W);
if (rem == 0) {
std::move_backward(a.begin(), a.end() - blk, a.end());
} else {
for (std::size_t i = a.size() - 1; i > blk; --i)
a[i] = (a[i - blk] << rem) | (a[i - blk - 1] >> (W - rem));
a[blk] = a[0] << rem;
}
std::fill(a.begin(), a.begin() + blk, word_type{});
trim();
return *this;
}
dynamic_bitset &operator>>=(std::size_t k) {
if (k == 0)
return *this;
if (k >= sz_)
return reset();
std::size_t blk = k / W;
std::size_t count = a.size() - blk;
int rem = int(k % W);
if (rem == 0) {
std::move(a.begin() + blk, a.end(), a.begin());
} else {
for (std::size_t i = 0; i + 1 < count; ++i)
a[i] = (a[i + blk] >> rem) | (a[i + blk + 1] << (W - rem));
a[count - 1] = a.back() >> rem;
}
std::fill(a.begin() + count, a.end(), word_type{});
return *this;
}
friend dynamic_bitset operator&(dynamic_bitset arr, const dynamic_bitset &b) {
arr &= b;
return arr;
}
friend dynamic_bitset operator|(dynamic_bitset arr, const dynamic_bitset &b) {
arr |= b;
return arr;
}
friend dynamic_bitset operator^(dynamic_bitset arr, const dynamic_bitset &b) {
arr ^= b;
return arr;
}
friend dynamic_bitset operator<<(dynamic_bitset val, std::size_t k) {
val <<= k;
return val;
}
friend dynamic_bitset operator>>(dynamic_bitset val, std::size_t k) {
val >>= k;
return val;
}
friend dynamic_bitset operator~(dynamic_bitset val) {
val.flip();
return val;
}
friend bool operator==(const dynamic_bitset &,
const dynamic_bitset &) = default;
private:
static std::size_t word_count(std::size_t siz) { return (siz + W - 1) / W; }
void trim() {
if (!a.empty() && sz_ % W != 0) {
a.back() &= (word_type(1) << (sz_ % W)) - word_type(1);
}
}
void check_same_size(const dynamic_bitset &rhs) const {
assert(sz_ == rhs.sz_);
}
};
} // namespace noya
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <span>
#include <utility>
#include <vector>
namespace noya {
struct bipartite_decomposition_result {
std::vector<int> ml;
std::vector<int> mr;
std::vector<int> cvl;
std::vector<int> cvr;
std::vector<int> isl;
std::vector<int> isr;
bool ok = false;
std::vector<std::pair<int, int>> ec;
int matching_size() const {
return int(
std::count_if(ml.begin(), ml.end(), [](int val) { return val != -1; }));
}
};
namespace bipartite_decomposition_internal {
struct adjacency {
std::vector<std::size_t> off;
std::vector<int> to;
adjacency(int nl, int nr, const std::vector<std::pair<int, int>> &es)
: off(std::size_t(nl) + 1), to(es.size()) {
for (auto [l, r] : es) {
assert(0 <= l && l < nl && 0 <= r && r < nr);
++off[l + 1];
}
for (int l = 0; l < nl; ++l) off[l + 1] += off[l];
auto cur = off;
for (auto [l, r] : es) to[cur[l]++] = r;
}
explicit adjacency(const std::vector<dynamic_bitset> &rows)
: off(rows.size() + 1) {
for (std::size_t u = 0; u < rows.size(); ++u)
off[u + 1] = off[u] + rows[u].count();
to.reserve(off.back());
for (const auto &row : rows)
for (std::size_t v = row.find_first(); v < row.size(); v = row.find_next(v + 1))
to.push_back(int(v));
}
int size() const { return int(off.size() - 1); }
std::span<const int> operator[](int u) const {
return std::span(to).subspan(off[u], off[u + 1] - off[u]);
}
};
// Return final alternating-reachability distances from unmatched left vertices.
// Per-phase cursors scan each edge at most once during layered augmentation.
inline std::vector<int> match(const adjacency &g, int nr, std::vector<int> &ml,
std::vector<int> &mr) {
int nl = g.size();
ml.assign(nl, -1);
mr.assign(nr, -1);
std::vector<int> active;
for (int l = 0; l < nl; ++l) {
if (g[l].empty()) continue;
active.push_back(l);
for (int r : g[l]) {
if (mr[r] == -1) {
ml[l] = r;
mr[r] = l;
break;
}
}
}
std::vector<int> dis(nl), q;
std::vector<std::size_t> cur;
q.reserve(active.size());
while (true) {
// Isolates keep distance zero: they are unmatched and trivially reachable.
// Only O(E) active vertices need resetting/scanning in each matching phase.
for (int l : active) dis[l] = -1;
q.clear();
for (int l : active) {
if (ml[l] == -1) {
dis[l] = 0;
q.push_back(l);
}
}
int limit = -1;
for (std::size_t head = 0; head < q.size(); ++head) {
int l = q[head];
if (limit != -1 && dis[l] >= limit) continue;
for (int r : g[l]) {
int next = mr[r];
if (next == -1) {
limit = dis[l] + 1;
} else if (dis[next] == -1) {
dis[next] = dis[l] + 1;
q.push_back(next);
}
}
}
if (limit == -1) return dis;
cur.resize(nl);
for (int l : active) cur[l] = 0;
auto dfs = [&](auto &self, int l) -> bool {
auto adj = g[l];
while (cur[l] < adj.size()) {
int r = adj[cur[l]++], next = mr[r];
if (next == -1 ? dis[l] + 1 == limit
: dis[next] == dis[l] + 1 && dis[next] < limit && self(self, next)) {
ml[l] = r;
mr[r] = l;
return true;
}
}
dis[l] = -1;
return false;
};
for (int l : active) if (ml[l] == -1) dfs(dfs, l);
}
}
inline std::vector<int> topology(const adjacency &g) {
int n = g.size();
std::vector<std::size_t> deg(n);
std::vector<int> top;
top.reserve(n);
for (int u = 0; u < n; ++u)
for (int v : g[u]) ++deg[v];
for (int u = 0; u < n; ++u)
if (deg[u] == 0) top.push_back(u);
for (std::size_t at = 0; at < top.size(); ++at)
for (int v : g[top[at]])
if (--deg[v] == 0) top.push_back(v);
assert(int(top.size()) == n);
return top;
}
inline std::vector<std::vector<int>> chains(const std::vector<int> &ml,
const std::vector<int> &mr) {
std::vector<std::vector<int>> res;
for (int s = 0; s < int(mr.size()); ++s) {
if (mr[s] != -1) continue;
res.emplace_back();
for (int u = s; u != -1; u = ml[u]) res.back().push_back(u);
}
return res;
}
} // namespace bipartite_decomposition_internal
/// @brief Maximum matching, Konig minimum vertex cover and its complementary
/// maximum independent set. Also returns a minimum edge cover iff no vertex is
/// isolated (ok=true); otherwise ec is empty. Duplicate edges are allowed.
/// Vertex sets are sorted; matching/edge-cover witnesses need not be canonical.
inline bipartite_decomposition_result
bipartite_decomposition(int nl, int nr,
const std::vector<std::pair<int, int>> &es) {
assert(nl >= 0 && nr >= 0);
using namespace bipartite_decomposition_internal;
adjacency g(nl, nr, es);
bipartite_decomposition_result res;
auto dis = match(g, nr, res.ml, res.mr);
for (int l = 0; l < nl; ++l)
(dis[l] >= 0 ? res.isl : res.cvl).push_back(l);
// At maximum matching every reachable right vertex is matched, and its
// matched left partner is reachable too. No second graph traversal is needed.
for (int r = 0; r < nr; ++r)
(res.mr[r] >= 0 && dis[res.mr[r]] >= 0 ? res.cvr : res.isr).push_back(r);
std::vector<int> first(nr, -1);
for (auto [l, r] : es)
if (first[r] == -1) first[r] = l;
res.ok = true;
for (int l = 0; l < nl; ++l) res.ok &= !g[l].empty();
for (int r = 0; r < nr; ++r) res.ok &= first[r] != -1;
if (res.ok) {
res.ec.reserve(std::size_t(nl) + nr - res.cvl.size() - res.cvr.size());
for (int l = 0; l < nl; ++l)
if (res.ml[l] != -1) res.ec.emplace_back(l, res.ml[l]);
for (int l = 0; l < nl; ++l)
if (res.ml[l] == -1) res.ec.emplace_back(l, g[l][0]);
for (int r = 0; r < nr; ++r)
if (res.mr[r] == -1) res.ec.emplace_back(first[r], r);
}
return res;
}
/// @brief Minimum vertex-disjoint directed path cover of a DAG. Consecutive
/// vertices in a returned path are connected by an original edge. Duplicate
/// edges are allowed; the input must be acyclic.
inline std::vector<std::vector<int>>
minimum_dag_path_cover(int n, const std::vector<std::pair<int, int>> &es) {
assert(n >= 0);
using namespace bipartite_decomposition_internal;
adjacency g(n, n, es);
topology(g);
std::vector<int> ml, mr;
match(g, n, ml, mr);
return chains(ml, mr);
}
struct dag_order_decomposition_result {
std::vector<int> ac;
std::vector<std::vector<int>> cc;
};
/// @brief Maximum reachability antichain and minimum chain cover of a DAG.
/// Chain neighbors need only be reachable, unlike minimum_dag_path_cover.
/// Bitset closure takes O((n+E) ceil(n/64)); matching then takes O(n+R sqrt(n))
/// for R reachable ordered pairs. The input must be acyclic.
inline dag_order_decomposition_result
dag_order_decomposition(int n, const std::vector<std::pair<int, int>> &es) {
assert(n >= 0);
using namespace bipartite_decomposition_internal;
adjacency g(n, n, es);
auto top = topology(g);
std::vector<dynamic_bitset> reach(n, dynamic_bitset(n));
for (auto it = top.rbegin(); it != top.rend(); ++it) {
for (int v : g[*it]) {
if (reach[*it][v]) continue;
reach[*it].set(v);
reach[*it] |= reach[v];
}
}
adjacency comparable(reach);
std::vector<int> ml, mr;
auto dis = match(comparable, n, ml, mr);
dag_order_decomposition_result res;
for (int u = 0; u < n; ++u)
if (dis[u] >= 0 && (mr[u] == -1 || dis[mr[u]] < 0)) res.ac.push_back(u);
res.cc = chains(ml, mr);
return res;
}
} // namespace noya
// NOYA_BIPARTITE_DECOMPOSITION_HPP