pfaffian.hpp¶
计算斜对称矩阵的 Pfaffian；用于平方等于行列式的配对计数公式。
\[
\displaystyle \operatorname{pf}(A)^2=\det(A)
\]
Complexity: Time: \(O(n^{3})\). Space: \(O(n^{2})\).
AC 记录：pfaffian_of_matrix。
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: O(n^3).
/// Space: O(n^2).
#include <algorithm>
#include <cassert>
#include <type_traits>
#include <vector>
namespace noya {
/// @brief Compute the Pfaffian of an even-order alternating matrix over a
/// field. Each step moves a nonzero entry into the next 2-by-2 pivot block;
/// the skew-symmetric Schur complement removes that pair, and the product of
/// pivot entries is the Pfaffian. Swapping a paired index changes its sign.
template <class T> T pfaffian(std::vector<std::vector<T>> mat) {
int sz = int(mat.size());
assert(sz % 2 == 0);
for (const auto &row : mat) {
assert(int(row.size()) == sz);
}
T res = T(1);
for (int a = 0; a < sz; a += 2) {
int piv = a + 1;
while (piv < sz && mat[a][piv] == T{}) {
piv++;
}
if (piv == sz) {
return T{};
}
if (piv != a + 1) {
std::swap(mat[piv], mat[a + 1]);
for (auto &row : mat) {
std::swap(row[piv], row[a + 1]);
}
res = -res;
}
T pv = mat[a][a + 1];
res *= pv;
if (a + 2 == sz) return res;
T inv = T(1) / pv;
for (int row = a + 2; row + 1 < sz; row++) {
T x = mat[a][row], y = mat[a + 1][row];
if constexpr (!std::is_floating_point_v<T>) {
x *= inv;
y *= inv;
}
for (int col = row + 1; col < sz; col++) {
if constexpr (std::is_floating_point_v<T>) {
// Scaling first can underflow before the other factor is applied.
mat[row][col] -= (x * mat[a + 1][col] - mat[a][col] * y) * inv;
} else {
mat[row][col] -= x * mat[a + 1][col] - mat[a][col] * y;
}
mat[col][row] = -mat[row][col];
}
}
}
return res;
}
} // namespace noya
// NOYA_PFAFFIAN_HPP
#ifndef NOYA_PFAFFIAN_HPP
#define NOYA_PFAFFIAN_HPP 1
/// @complexity Time: O(n^3).
/// Space: O(n^2).
#include <algorithm>
#include <cassert>
#include <type_traits>
#include <vector>
namespace noya {
/// @brief Compute the Pfaffian of an even-order alternating matrix over a
/// field. Each step moves a nonzero entry into the next 2-by-2 pivot block;
/// the skew-symmetric Schur complement removes that pair, and the product of
/// pivot entries is the Pfaffian. Swapping a paired index changes its sign.
template <class T> T pfaffian(std::vector<std::vector<T>> mat) {
int sz = int(mat.size());
assert(sz % 2 == 0);
for (const auto &row : mat) {
assert(int(row.size()) == sz);
}
T res = T(1);
for (int a = 0; a < sz; a += 2) {
int piv = a + 1;
while (piv < sz && mat[a][piv] == T{}) {
piv++;
}
if (piv == sz) {
return T{};
}
if (piv != a + 1) {
std::swap(mat[piv], mat[a + 1]);
for (auto &row : mat) {
std::swap(row[piv], row[a + 1]);
}
res = -res;
}
T pv = mat[a][a + 1];
res *= pv;
if (a + 2 == sz) return res;
T inv = T(1) / pv;
for (int row = a + 2; row + 1 < sz; row++) {
T x = mat[a][row], y = mat[a + 1][row];
if constexpr (!std::is_floating_point_v<T>) {
x *= inv;
y *= inv;
}
for (int col = row + 1; col < sz; col++) {
if constexpr (std::is_floating_point_v<T>) {
// Scaling first can underflow before the other factor is applied.
mat[row][col] -= (x * mat[a + 1][col] - mat[a][col] * y) * inv;
} else {
mat[row][col] -= x * mat[a + 1][col] - mat[a][col] * y;
}
mat[col][row] = -mat[row][col];
}
}
}
return res;
}
} // namespace noya
#endif // NOYA_PFAFFIAN_HPP
/// @complexity Time: O(n^3).
/// Space: O(n^2).
#include <algorithm>
#include <cassert>
#include <type_traits>
#include <vector>
namespace noya {
/// @brief Compute the Pfaffian of an even-order alternating matrix over a
/// field. Each step moves a nonzero entry into the next 2-by-2 pivot block;
/// the skew-symmetric Schur complement removes that pair, and the product of
/// pivot entries is the Pfaffian. Swapping a paired index changes its sign.
template <class T> T pfaffian(std::vector<std::vector<T>> mat) {
int sz = int(mat.size());
assert(sz % 2 == 0);
for (const auto &row : mat) {
assert(int(row.size()) == sz);
}
T res = T(1);
for (int a = 0; a < sz; a += 2) {
int piv = a + 1;
while (piv < sz && mat[a][piv] == T{}) {
piv++;
}
if (piv == sz) {
return T{};
}
if (piv != a + 1) {
std::swap(mat[piv], mat[a + 1]);
for (auto &row : mat) {
std::swap(row[piv], row[a + 1]);
}
res = -res;
}
T pv = mat[a][a + 1];
res *= pv;
if (a + 2 == sz) return res;
T inv = T(1) / pv;
for (int row = a + 2; row + 1 < sz; row++) {
T x = mat[a][row], y = mat[a + 1][row];
if constexpr (!std::is_floating_point_v<T>) {
x *= inv;
y *= inv;
}
for (int col = row + 1; col < sz; col++) {
if constexpr (std::is_floating_point_v<T>) {
// Scaling first can underflow before the other factor is applied.
mat[row][col] -= (x * mat[a + 1][col] - mat[a][col] * y) * inv;
} else {
mat[row][col] -= x * mat[a + 1][col] - mat[a][col] * y;
}
mat[col][row] = -mat[row][col];
}
}
}
return res;
}
} // namespace noya
// NOYA_PFAFFIAN_HPP