best_theorem.hpp¶
用 BEST 定理计算欧拉有向图中欧拉回路数量，结合有向矩阵树定理处理入树计数。
Complexity: Time: O(V^3 + E) field operations. Space: O(V^2 + E).
AC 记录：counting_eulerian_circuits。
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: O(V^3 + E) field operations.
/// Space: O(V^2 + E).
#include "noya/matrix_tree.hpp"
#include <cassert>
#include <tuple>
#include <utility>
#include <vector>
namespace noya {
/// @brief Count Euler circuits of an Eulerian directed multigraph by the BEST
/// theorem. Delete the root row and column of the outgoing Laplacian to count
/// in-arborescences rooted there, then multiply by `(dou(v) - 1)!` for
/// every active vertex. This orders the remaining labeled outgoing edges after
/// the arborescence chooses the last exit from each vertex. With
/// `fix=true`, the first outgoing labeled edge at `rt` is
/// prescribed; otherwise all choices for that first edge are counted.
template <class T>
T best_euler_circuit_count(int n, int rt,
const std::vector<std::pair<int, int>> &es,
bool fix = true) {
assert(n > 0);
assert(0 <= rt && rt < n);
if (es.empty()) {
return T(1);
}
std::vector<int> din(n);
std::vector<int> dou(n);
for (auto [u1, to] : es) {
assert(0 <= u1 && u1 < n);
assert(0 <= to && to < n);
dou[u1]++;
din[to]++;
}
for (int u = 0; u < n; u++) {
if (din[u] != dou[u]) {
return T(0);
}
}
if (dou[rt] == 0) {
return T(0);
}
std::vector<int> id(n, -1);
int cnt = 0;
for (int u = 0; u < n; u++) {
if (dou[u] > 0) {
id[u] = cnt++;
}
}
std::vector<std::tuple<int, int, T>> es1;
es1.reserve(es.size());
for (auto [u1, to] : es) {
es1.emplace_back(id[u1], id[to], T(1));
}
T res = in_arborescence_count(cnt, id[rt], es1);
for (int u = 0; u < n; u++) {
for (int fac = 2; fac < dou[u]; fac++) {
res *= T(fac);
}
}
if (!fix) {
res *= T(dou[rt]);
}
return res;
}
} // namespace noya
#ifndef NOYA_BEST_THEOREM_HPP
#define NOYA_BEST_THEOREM_HPP 1
/// @complexity Time: O(V^3 + E) field operations.
/// Space: O(V^2 + E).
#include "noya/matrix_tree.hpp"
#include <cassert>
#include <tuple>
#include <utility>
#include <vector>
namespace noya {
/// @brief Count Euler circuits of an Eulerian directed multigraph by the BEST
/// theorem. Delete the root row and column of the outgoing Laplacian to count
/// in-arborescences rooted there, then multiply by `(dou(v) - 1)!` for
/// every active vertex. This orders the remaining labeled outgoing edges after
/// the arborescence chooses the last exit from each vertex. With
/// `fix=true`, the first outgoing labeled edge at `rt` is
/// prescribed; otherwise all choices for that first edge are counted.
template <class T>
T best_euler_circuit_count(int n, int rt,
const std::vector<std::pair<int, int>> &es,
bool fix = true) {
assert(n > 0);
assert(0 <= rt && rt < n);
if (es.empty()) {
return T(1);
}
std::vector<int> din(n);
std::vector<int> dou(n);
for (auto [u1, to] : es) {
assert(0 <= u1 && u1 < n);
assert(0 <= to && to < n);
dou[u1]++;
din[to]++;
}
for (int u = 0; u < n; u++) {
if (din[u] != dou[u]) {
return T(0);
}
}
if (dou[rt] == 0) {
return T(0);
}
std::vector<int> id(n, -1);
int cnt = 0;
for (int u = 0; u < n; u++) {
if (dou[u] > 0) {
id[u] = cnt++;
}
}
std::vector<std::tuple<int, int, T>> es1;
es1.reserve(es.size());
for (auto [u1, to] : es) {
es1.emplace_back(id[u1], id[to], T(1));
}
T res = in_arborescence_count(cnt, id[rt], es1);
for (int u = 0; u < n; u++) {
for (int fac = 2; fac < dou[u]; fac++) {
res *= T(fac);
}
}
if (!fix) {
res *= T(dou[rt]);
}
return res;
}
} // namespace noya
#endif // NOYA_BEST_THEOREM_HPP
#include <algorithm>
#include <cassert>
#include <tuple>
#include <utility>
#include <vector>
/// @complexity Time: O(V^3 + E) field operations.
/// Space: O(V^2 + E).
/// @complexity Time: O(V^3 + E) field operations.
/// Space: O(V^2).
/// @complexity Time: O(rows * columns * min(rows,columns)) elimination; O(n^3) square determinant/inverse.
/// Space: O(rows * columns).
namespace noya {
/// @brief Exact zero predicate used by elimination routines by default.
template <class T> struct exact_zero {
bool operator()(const T &val) const { return val == T{}; }
};
/// @brief Consistency flag, one solution, nullspace basis, and pivot columns.
template <class T> struct linear_system_solution {
bool ok = false;
std::vector<T> sol;
std::vector<std::vector<T>> ker;
std::vector<int> pc;
};
namespace linear_algebra_internal {
template <class T> int column_count(const std::vector<std::vector<T>> &mat) {
if (mat.empty()) {
return 0;
}
int cs = int(mat[0].size());
for (const auto &row : mat) {
assert(int(row.size()) == cs);
}
return cs;
}
} // namespace linear_algebra_internal
/// @brief Compute matrix rank over a field.
template <class T, class IsZero = exact_zero<T>>
int matrix_rank(std::vector<std::vector<T>> mat, IsZero iz = {}) {
int rs = int(mat.size());
int cs = linear_algebra_internal::column_count(mat);
int ran = 0;
for (int col = 0; col < cs && ran < rs; col++) {
int piv = ran;
while (piv < rs && iz(mat[piv][col])) {
piv++;
}
if (piv == rs) {
continue;
}
std::swap(mat[piv], mat[ran]);
for (int row = ran + 1; row < rs; row++) {
if (iz(mat[row][col])) {
continue;
}
T r = mat[row][col] / mat[ran][col];
for (int j = col; j < cs; j++) {
mat[row][j] -= r * mat[ran][j];
}
}
ran++;
}
return ran;
}
/// @brief Compute the determinant of a square matrix over a field.
template <class T, class IsZero = exact_zero<T>>
T determinant(std::vector<std::vector<T>> mat, IsZero iz = {}) {
int n = int(mat.size());
assert(linear_algebra_internal::column_count(mat) == n);
T res = T(1);
for (int col = 0; col < n; col++) {
int piv = col;
while (piv < n && iz(mat[piv][col])) {
piv++;
}
if (piv == n) {
return T{};
}
if (piv != col) {
std::swap(mat[piv], mat[col]);
res = -res;
}
T pv = mat[col][col];
res *= pv;
for (int row = col + 1; row < n; row++) {
if (iz(mat[row][col])) {
continue;
}
T r = mat[row][col] / pv;
for (int j = col; j < n; j++) {
mat[row][j] -= r * mat[col][j];
}
}
}
return res;
}
/// @brief Solve A*x=b and return one solution plus a basis of the nullspace.
template <class T, class IsZero = exact_zero<T>>
linear_system_solution<T> solve_linear(std::vector<std::vector<T>> mat,
std::vector<T> rhs, IsZero iz = {}) {
int rs = int(mat.size());
assert(int(rhs.size()) == rs);
int cs = linear_algebra_internal::column_count(mat);
std::vector<int> pc;
int ran = 0;
for (int col = 0; col < cs && ran < rs; col++) {
int piv = ran;
while (piv < rs && iz(mat[piv][col])) {
piv++;
}
if (piv == rs) {
continue;
}
std::swap(mat[piv], mat[ran]);
std::swap(rhs[piv], rhs[ran]);
T inv = T(1) / mat[ran][col];
for (int j = col; j < cs; j++) {
mat[ran][j] *= inv;
}
rhs[ran] *= inv;
for (int row = 0; row < rs; row++) {
if (row == ran || iz(mat[row][col])) {
continue;
}
T r = mat[row][col];
for (int j = col; j < cs; j++) {
mat[row][j] -= r * mat[ran][j];
}
rhs[row] -= r * rhs[ran];
}
pc.push_back(col);
ran++;
}
for (int row = ran; row < rs; row++) {
bool az = true;
for (int col = 0; col < cs; col++) {
az &= iz(mat[row][col]);
}
if (az && !iz(rhs[row])) {
return {};
}
}
linear_system_solution<T> res;
res.ok = true;
res.sol.assign(cs, T{});
res.pc = pc;
std::vector<bool> ip(cs);
for (int row = 0; row < ran; row++) {
int col = pc[row];
ip[col] = true;
res.sol[col] = rhs[row];
}
for (int fc = 0; fc < cs; fc++) {
if (ip[fc]) {
continue;
}
std::vector<T> bv(cs, T{});
bv[fc] = T(1);
for (int row = 0; row < ran; row++) {
bv[pc[row]] = -mat[row][fc];
}
res.ker.push_back(std::move(bv));
}
return res;
}
} // namespace noya
namespace noya {
/// @brief Count weighted spanning trees of an undirected graph over a field
/// using the Matrix-Tree theorem in O(n^3) time.
template <class T>
T spanning_tree_count(int n, const std::vector<std::tuple<int, int, T>> &es) {
assert(n >= 0);
if (n <= 1) {
return T(1);
}
std::vector<std::vector<T>> mat(n - 1, std::vector<T>(n - 1));
for (const auto &[a, b, w] : es) {
assert(0 <= a && a < n);
assert(0 <= b && b < n);
if (a == b) {
continue;
}
if (a < n - 1) {
mat[a][a] += w;
}
if (b < n - 1) {
mat[b][b] += w;
}
if (a < n - 1 && b < n - 1) {
mat[a][b] -= w;
mat[b][a] -= w;
}
}
return determinant(std::move(mat));
}
/// @brief Count weighted directed spanning trees whose edges point toward root
/// over a field, using the directed Matrix-Tree theorem in O(n^3) time.
template <class T>
T in_arborescence_count(int n, int rt,
const std::vector<std::tuple<int, int, T>> &es) {
assert(n > 0);
assert(0 <= rt && rt < n);
if (n == 1) {
return T(1);
}
std::vector<std::vector<T>> lap(n, std::vector<T>(n));
for (const auto &[u, to, w] : es) {
assert(0 <= u && u < n);
assert(0 <= to && to < n);
if (u == to) {
continue;
}
lap[u][u] += w;
lap[u][to] -= w;
}
std::vector<std::vector<T>> mat;
mat.reserve(n - 1);
for (int row = 0; row < n; row++) {
if (row == rt) {
continue;
}
mat.emplace_back();
mat.back().reserve(n - 1);
for (int col = 0; col < n; col++) {
if (col != rt) {
mat.back().push_back(lap[row][col]);
}
}
}
return determinant(std::move(mat));
}
} // namespace noya
namespace noya {
/// @brief Count Euler circuits of an Eulerian directed multigraph by the BEST
/// theorem. Delete the root row and column of the outgoing Laplacian to count
/// in-arborescences rooted there, then multiply by `(dou(v) - 1)!` for
/// every active vertex. This orders the remaining labeled outgoing edges after
/// the arborescence chooses the last exit from each vertex. With
/// `fix=true`, the first outgoing labeled edge at `rt` is
/// prescribed; otherwise all choices for that first edge are counted.
template <class T>
T best_euler_circuit_count(int n, int rt,
const std::vector<std::pair<int, int>> &es,
bool fix = true) {
assert(n > 0);
assert(0 <= rt && rt < n);
if (es.empty()) {
return T(1);
}
std::vector<int> din(n);
std::vector<int> dou(n);
for (auto [u1, to] : es) {
assert(0 <= u1 && u1 < n);
assert(0 <= to && to < n);
dou[u1]++;
din[to]++;
}
for (int u = 0; u < n; u++) {
if (din[u] != dou[u]) {
return T(0);
}
}
if (dou[rt] == 0) {
return T(0);
}
std::vector<int> id(n, -1);
int cnt = 0;
for (int u = 0; u < n; u++) {
if (dou[u] > 0) {
id[u] = cnt++;
}
}
std::vector<std::tuple<int, int, T>> es1;
es1.reserve(es.size());
for (auto [u1, to] : es) {
es1.emplace_back(id[u1], id[to], T(1));
}
T res = in_arborescence_count(cnt, id[rt], es1);
for (int u = 0; u < n; u++) {
for (int fac = 2; fac < dou[u]; fac++) {
res *= T(fac);
}
}
if (!fix) {
res *= T(dou[rt]);
}
return res;
}
} // namespace noya