hafnian.hpp¶
Compute the hafnian of an even-order symmetric matrix over a commutative ring. A paired-vertex elimination keeps every remaining edge weight as a polynomial in a marker z. Eliminating vertices u,v adds the two possible paths i-u-v-j to each surviving edge; an include/exclude recurrence cancels structures that reuse a vertex. After n/2 eliminations, the coefficient of z^(n/2) contains exactly the weighted perfect matchings.
Verified by hafnian_of_matrix.
\[
\displaystyle \operatorname{haf}(A)^2=\det(A)
\]
Implementation¶
#ifndef NOYA_HAFNIAN_HPP
#define NOYA_HAFNIAN_HPP 1
/// @complexity Time: O(n^2 2^(n/2)) field operations.
/// Space: O(n^4) with the natural recursive polynomial states.
#include <cassert>
#include <utility>
#include <vector>
namespace noya {
namespace hafnian_internal {
template <class T> struct calculator {
using polynomial = std::vector<T>;
int degree_count;
explicit calculator(int degree_count_) : degree_count(degree_count_) {}
void multiply_and_shift_add(polynomial &result, const polynomial &first,
const polynomial &second) const {
for (int i = 0; i < degree_count; i++) {
for (int j = 0; i + j + 1 < degree_count; j++) {
result[i + j + 1] += first[i] * second[j];
}
}
}
polynomial eliminate(
std::vector<std::vector<polynomial>> weights) const {
if (weights.empty()) {
polynomial result(degree_count);
result[0] = T(1);
return result;
}
std::vector<polynomial> last = std::move(weights.back());
weights.pop_back();
std::vector<polynomial> second_last = std::move(weights.back());
weights.pop_back();
int remaining = int(weights.size());
polynomial result = eliminate(weights);
for (T &coefficient : result) {
coefficient = -coefficient;
}
for (int first = 0; first < remaining; first++) {
for (int second = 0; second < first; second++) {
multiply_and_shift_add(weights[first][second], last[first],
second_last[second]);
multiply_and_shift_add(weights[first][second], second_last[first],
last[second]);
}
}
polynomial linked = eliminate(std::move(weights));
multiply_and_shift_add(result, last[remaining], linked);
for (int degree = 0; degree < degree_count; degree++) {
result[degree] += linked[degree];
}
return result;
}
};
} // namespace hafnian_internal
/// @brief Compute the hafnian of an even-order symmetric matrix over a
/// commutative ring. A paired-vertex elimination keeps every remaining edge
/// weight as a polynomial in a marker z. Eliminating vertices u,v adds the two
/// possible paths i-u-v-j to each surviving edge; an include/exclude recurrence
/// cancels structures that reuse a vertex. After n/2 eliminations, the
/// coefficient of z^(n/2) contains exactly the weighted perfect matchings.
template <class T> T hafnian(const 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);
}
int degree_count = size / 2 + 1;
std::vector<std::vector<std::vector<T>>> weights(size);
for (int row = 0; row < size; row++) {
weights[row].assign(row, std::vector<T>(degree_count));
for (int column = 0; column < row; column++) {
assert(matrix[row][column] == matrix[column][row]);
weights[row][column][0] = matrix[row][column];
}
}
return hafnian_internal::calculator<T>(degree_count)
.eliminate(std::move(weights))
.back();
}
} // namespace noya
#endif // NOYA_HAFNIAN_HPP
#include <cassert>
#include <utility>
#include <vector>
/// @complexity Time: O(n^2 2^(n/2)) field operations.
/// Space: O(n^4) with the natural recursive polynomial states.
namespace noya {
namespace hafnian_internal {
template <class T> struct calculator {
using polynomial = std::vector<T>;
int degree_count;
explicit calculator(int degree_count_) : degree_count(degree_count_) {}
void multiply_and_shift_add(polynomial &result, const polynomial &first,
const polynomial &second) const {
for (int i = 0; i < degree_count; i++) {
for (int j = 0; i + j + 1 < degree_count; j++) {
result[i + j + 1] += first[i] * second[j];
}
}
}
polynomial eliminate(
std::vector<std::vector<polynomial>> weights) const {
if (weights.empty()) {
polynomial result(degree_count);
result[0] = T(1);
return result;
}
std::vector<polynomial> last = std::move(weights.back());
weights.pop_back();
std::vector<polynomial> second_last = std::move(weights.back());
weights.pop_back();
int remaining = int(weights.size());
polynomial result = eliminate(weights);
for (T &coefficient : result) {
coefficient = -coefficient;
}
for (int first = 0; first < remaining; first++) {
for (int second = 0; second < first; second++) {
multiply_and_shift_add(weights[first][second], last[first],
second_last[second]);
multiply_and_shift_add(weights[first][second], second_last[first],
last[second]);
}
}
polynomial linked = eliminate(std::move(weights));
multiply_and_shift_add(result, last[remaining], linked);
for (int degree = 0; degree < degree_count; degree++) {
result[degree] += linked[degree];
}
return result;
}
};
} // namespace hafnian_internal
/// @brief Compute the hafnian of an even-order symmetric matrix over a
/// commutative ring. A paired-vertex elimination keeps every remaining edge
/// weight as a polynomial in a marker z. Eliminating vertices u,v adds the two
/// possible paths i-u-v-j to each surviving edge; an include/exclude recurrence
/// cancels structures that reuse a vertex. After n/2 eliminations, the
/// coefficient of z^(n/2) contains exactly the weighted perfect matchings.
template <class T> T hafnian(const 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);
}
int degree_count = size / 2 + 1;
std::vector<std::vector<std::vector<T>>> weights(size);
for (int row = 0; row < size; row++) {
weights[row].assign(row, std::vector<T>(degree_count));
for (int column = 0; column < row; column++) {
assert(matrix[row][column] == matrix[column][row]);
weights[row][column][0] = matrix[row][column];
}
}
return hafnian_internal::calculator<T>(degree_count)
.eliminate(std::move(weights))
.back();
}
} // namespace noya