best_theorem.hpp¶
用 BEST 定理统计有标号边的有向欧拉回路；直接构造矩阵树余子式，选择非自环高度数根并共享阶乘计算，支持固定首边或指定起点的全部首边。
Complexity: Time: \(O(V + E + A^{3})\), where \(A\) is the number of active vertices. Space: \(O(V + A^{2})\). Field operations are assumed constant time.
AC 记录：counting_eulerian_circuits。
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: O(V + E + A^3), where A is the number of active vertices.
/// Space: O(V + A^2). Field operations are assumed constant time.
#include "noya/math/linear_algebra.hpp"
#include <algorithm>
#include <cassert>
#include <utility>
#include <vector>
namespace noya {
/// @brief Count Euler circuits of an Eulerian directed multigraph by the BEST
/// theorem over a field. Parallel edges are labeled; loops affect factorials
/// but not the Laplacian, and isolated vertices are omitted. An empty edge set
/// returns 1; otherwise an inactive root or a non-Eulerian graph returns 0.
/// With `fix=true`, one outgoing labeled edge at `rt` is prescribed first;
/// otherwise all choices of the first edge at the original `rt` 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);
std::vector<int> deg(n);
for (auto [u1, to] : es) {
assert(0 <= u1 && u1 < n);
assert(0 <= to && to < n);
dou[u1]++;
din[to]++;
if (u1 != to) {
deg[u1]++;
}
}
if (din != dou || dou[rt] == 0) {
return T(0);
}
int root = rt;
for (int u = 0; u < n; u++) {
if (deg[u] > deg[root]) {
root = u;
}
}
// Eulerian in-tree counts are root-independent. Deleting a high-degree
// vertex can reduce elimination fill; loops must not influence this choice.
T res = fix ? T(1) : T(dou[rt]);
std::vector<int> id = std::move(din);
int cnt = 0, used = 0;
for (int u = 0; u < n; u++) {
int d = dou[u];
id[u] = -1;
if (d > 0 && u != root) {
id[u] = cnt++;
}
if (d > 1) {
dou[used++] = d;
}
}
// Reuse the degree buffer and share factorial work without an O(max degree)
// table. Degrees 0 and 1 contribute no factor.
dou.resize(used);
std::sort(dou.begin(), dou.end());
T fac(1);
int next = 2;
for (int d : dou) {
while (next < d) {
fac *= T(next++);
}
res *= fac;
}
if (res == T(0)) {
return res;
}
std::vector<std::vector<T>> mat(cnt, std::vector<T>(cnt));
for (int u = 0; u < n; u++) {
if (id[u] >= 0) {
mat[id[u]][id[u]] = T(deg[u]);
}
}
for (auto [u, to] : es) {
if (u != to && id[u] >= 0 && id[to] >= 0) {
mat[id[u]][id[to]] -= T(1);
}
}
return res * determinant(std::move(mat));
}
} // namespace noya
// NOYA_BEST_THEOREM_HPP
#ifndef NOYA_BEST_THEOREM_HPP
#define NOYA_BEST_THEOREM_HPP 1
/// @complexity Time: O(V + E + A^3), where A is the number of active vertices.
/// Space: O(V + A^2). Field operations are assumed constant time.
#include "noya/math/linear_algebra.hpp"
#include <algorithm>
#include <cassert>
#include <utility>
#include <vector>
namespace noya {
/// @brief Count Euler circuits of an Eulerian directed multigraph by the BEST
/// theorem over a field. Parallel edges are labeled; loops affect factorials
/// but not the Laplacian, and isolated vertices are omitted. An empty edge set
/// returns 1; otherwise an inactive root or a non-Eulerian graph returns 0.
/// With `fix=true`, one outgoing labeled edge at `rt` is prescribed first;
/// otherwise all choices of the first edge at the original `rt` 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);
std::vector<int> deg(n);
for (auto [u1, to] : es) {
assert(0 <= u1 && u1 < n);
assert(0 <= to && to < n);
dou[u1]++;
din[to]++;
if (u1 != to) {
deg[u1]++;
}
}
if (din != dou || dou[rt] == 0) {
return T(0);
}
int root = rt;
for (int u = 0; u < n; u++) {
if (deg[u] > deg[root]) {
root = u;
}
}
// Eulerian in-tree counts are root-independent. Deleting a high-degree
// vertex can reduce elimination fill; loops must not influence this choice.
T res = fix ? T(1) : T(dou[rt]);
std::vector<int> id = std::move(din);
int cnt = 0, used = 0;
for (int u = 0; u < n; u++) {
int d = dou[u];
id[u] = -1;
if (d > 0 && u != root) {
id[u] = cnt++;
}
if (d > 1) {
dou[used++] = d;
}
}
// Reuse the degree buffer and share factorial work without an O(max degree)
// table. Degrees 0 and 1 contribute no factor.
dou.resize(used);
std::sort(dou.begin(), dou.end());
T fac(1);
int next = 2;
for (int d : dou) {
while (next < d) {
fac *= T(next++);
}
res *= fac;
}
if (res == T(0)) {
return res;
}
std::vector<std::vector<T>> mat(cnt, std::vector<T>(cnt));
for (int u = 0; u < n; u++) {
if (id[u] >= 0) {
mat[id[u]][id[u]] = T(deg[u]);
}
}
for (auto [u, to] : es) {
if (u != to && id[u] >= 0 && id[to] >= 0) {
mat[id[u]][id[to]] -= T(1);
}
}
return res * determinant(std::move(mat));
}
} // namespace noya
#endif // NOYA_BEST_THEOREM_HPP
/// @complexity Time: O(V + E + A^3), where A is the number of active vertices.
/// Space: O(V + A^2). Field operations are assumed constant time.
/// @complexity Time: O(nm min(n,m)) elimination and O(n^3) square determinant.
/// Solving also takes O(m(m-r)) to construct the nullspace basis.
/// Space: O(nm) for rank/determinant; O(nm + m(m-r)) for solving.
/// n and m are the row and column counts, and r is the rank.
/// Field operations are counted as O(1).
#include <algorithm>
#include <cassert>
#include <type_traits>
#include <utility>
#include <vector>
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 {
// Direct division avoids overflowing the reciprocal of a tiny floating pivot.
template <class T>
T pivot_ratio(const T &val, const T &piv, T &inv, bool &ready) {
if constexpr (std::is_floating_point_v<T>) return val / piv;
if (!ready) {
inv = T(1) / piv;
ready = true;
}
return val * inv;
}
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]);
T inv{};
bool ready = false;
for (int row = ran + 1; row < rs; row++) {
if (iz(mat[row][col])) {
continue;
}
T r = linear_algebra_internal::pivot_ratio(
mat[row][col], mat[ran][col], inv, ready);
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;
T inv{};
bool ready = false;
for (int row = col + 1; row < n; row++) {
if (iz(mat[row][col])) {
continue;
}
T r = linear_algebra_internal::pivot_ratio(mat[row][col], pv, inv, ready);
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
#include <algorithm>
#include <cassert>
#include <utility>
#include <vector>
namespace noya {
/// @brief Count Euler circuits of an Eulerian directed multigraph by the BEST
/// theorem over a field. Parallel edges are labeled; loops affect factorials
/// but not the Laplacian, and isolated vertices are omitted. An empty edge set
/// returns 1; otherwise an inactive root or a non-Eulerian graph returns 0.
/// With `fix=true`, one outgoing labeled edge at `rt` is prescribed first;
/// otherwise all choices of the first edge at the original `rt` 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);
std::vector<int> deg(n);
for (auto [u1, to] : es) {
assert(0 <= u1 && u1 < n);
assert(0 <= to && to < n);
dou[u1]++;
din[to]++;
if (u1 != to) {
deg[u1]++;
}
}
if (din != dou || dou[rt] == 0) {
return T(0);
}
int root = rt;
for (int u = 0; u < n; u++) {
if (deg[u] > deg[root]) {
root = u;
}
}
// Eulerian in-tree counts are root-independent. Deleting a high-degree
// vertex can reduce elimination fill; loops must not influence this choice.
T res = fix ? T(1) : T(dou[rt]);
std::vector<int> id = std::move(din);
int cnt = 0, used = 0;
for (int u = 0; u < n; u++) {
int d = dou[u];
id[u] = -1;
if (d > 0 && u != root) {
id[u] = cnt++;
}
if (d > 1) {
dou[used++] = d;
}
}
// Reuse the degree buffer and share factorial work without an O(max degree)
// table. Degrees 0 and 1 contribute no factor.
dou.resize(used);
std::sort(dou.begin(), dou.end());
T fac(1);
int next = 2;
for (int d : dou) {
while (next < d) {
fac *= T(next++);
}
res *= fac;
}
if (res == T(0)) {
return res;
}
std::vector<std::vector<T>> mat(cnt, std::vector<T>(cnt));
for (int u = 0; u < n; u++) {
if (id[u] >= 0) {
mat[id[u]][id[u]] = T(deg[u]);
}
}
for (auto [u, to] : es) {
if (u != to && id[u] >= 0 && id[to] >= 0) {
mat[id[u]][id[to]] -= T(1);
}
}
return res * determinant(std::move(mat));
}
} // namespace noya
// NOYA_BEST_THEOREM_HPP