pfaffian.hpp¶
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.
Verified by pfaffian_of_matrix.
\[
\displaystyle \operatorname{pf}(A)^2=\det(A)
\]
Implementation¶
#ifndef NOYA_PFAFFIAN_HPP
#define NOYA_PFAFFIAN_HPP 1
/// @complexity Time: O(n^3).
/// Space: O(n^2).
#include <algorithm>
#include <cassert>
#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>> matrix) {
int size = int(matrix.size());
assert(size % 2 == 0);
for (const auto &row : matrix) {
assert(int(row.size()) == size);
}
T result = T(1);
for (int first = 0; first < size; first += 2) {
int pivot = first + 1;
while (pivot < size && matrix[first][pivot] == T{}) {
pivot++;
}
if (pivot == size) {
return T{};
}
if (pivot != first + 1) {
std::swap(matrix[pivot], matrix[first + 1]);
for (auto &row : matrix) {
std::swap(row[pivot], row[first + 1]);
}
result = -result;
}
T pivot_value = matrix[first][first + 1];
result *= pivot_value;
T inverse = T(1) / pivot_value;
for (int row = first + 2; row < size; row++) {
for (int column = row + 1; column < size; column++) {
matrix[row][column] -=
(matrix[first][row] * matrix[first + 1][column] -
matrix[first][column] * matrix[first + 1][row]) *
inverse;
matrix[column][row] = -matrix[row][column];
}
}
}
return result;
}
} // namespace noya
#endif // NOYA_PFAFFIAN_HPP
#include <algorithm>
#include <cassert>
#include <vector>
/// @complexity Time: O(n^3).
/// Space: O(n^2).
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>> matrix) {
int size = int(matrix.size());
assert(size % 2 == 0);
for (const auto &row : matrix) {
assert(int(row.size()) == size);
}
T result = T(1);
for (int first = 0; first < size; first += 2) {
int pivot = first + 1;
while (pivot < size && matrix[first][pivot] == T{}) {
pivot++;
}
if (pivot == size) {
return T{};
}
if (pivot != first + 1) {
std::swap(matrix[pivot], matrix[first + 1]);
for (auto &row : matrix) {
std::swap(row[pivot], row[first + 1]);
}
result = -result;
}
T pivot_value = matrix[first][first + 1];
result *= pivot_value;
T inverse = T(1) / pivot_value;
for (int row = first + 2; row < size; row++) {
for (int column = row + 1; column < size; column++) {
matrix[row][column] -=
(matrix[first][row] * matrix[first + 1][column] -
matrix[first][column] * matrix[first + 1][row]) *
inverse;
matrix[column][row] = -matrix[row][column];
}
}
}
return result;
}
} // namespace noya