hafnian.hpp

SECTIONMath INCLUDEnoya/hafnian.hpp

计算对称矩阵的 Hafnian；用于带权完全匹配等不带符号的两两配对计数。

\[ \displaystyle \operatorname{haf}(A)^2=\det(A) \]

Complexity: Time: O(n^2 2^(n/2)) field operations. Space: O(n^4) with the natural recursive polynomial states.

AC 记录：hafnian_of_matrix。

跳到代码 · GitHub ↗

Implementation

ReferenceSourceExpanded

当前头文件，省略 include guard；依赖见 #include。

/// @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 nd;

  explicit calculator(int nd_) : nd(nd_) {}

  void multiply_and_shift_add(polynomial &res, const polynomial &a,
                              const polynomial &b) const {
    for (int i = 0; i < nd; i++) {
      for (int j = 0; i + j + 1 < nd; j++) {
        res[i + j + 1] += a[i] * b[j];
      }
    }
  }

  polynomial eliminate(std::vector<std::vector<polynomial>> wgt) const {
    if (wgt.empty()) {
      polynomial res(nd);
      res[0] = T(1);
      return res;
    }

    std::vector<polynomial> las = std::move(wgt.back());
    wgt.pop_back();
    std::vector<polynomial> sl = std::move(wgt.back());
    wgt.pop_back();
    int rem = int(wgt.size());

    polynomial res = eliminate(wgt);
    for (T &cf : res) {
      cf = -cf;
    }

    for (int a = 0; a < rem; a++) {
      for (int b = 0; b < a; b++) {
        multiply_and_shift_add(wgt[a][b], las[a], sl[b]);
        multiply_and_shift_add(wgt[a][b], sl[a], las[b]);
      }
    }
    polynomial lin = eliminate(std::move(wgt));
    multiply_and_shift_add(res, las[rem], lin);
    for (int deg = 0; deg < nd; deg++) {
      res[deg] += lin[deg];
    }
    return res;
  }
};

} // 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>> &mat) {
  int sz = int(mat.size());
  assert(sz % 2 == 0);
  for (const auto &row : mat) {
    assert(int(row.size()) == sz);
  }
  int nd = sz / 2 + 1;
  std::vector<std::vector<std::vector<T>>> wgt(sz);
  for (int row = 0; row < sz; row++) {
    wgt[row].assign(row, std::vector<T>(nd));
    for (int col = 0; col < row; col++) {
      assert(mat[row][col] == mat[col][row]);
      wgt[row][col][0] = mat[row][col];
    }
  }
  return hafnian_internal::calculator<T>(nd).eliminate(std::move(wgt)).back();
}

} // namespace noya

#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 nd;

  explicit calculator(int nd_) : nd(nd_) {}

  void multiply_and_shift_add(polynomial &res, const polynomial &a,
                              const polynomial &b) const {
    for (int i = 0; i < nd; i++) {
      for (int j = 0; i + j + 1 < nd; j++) {
        res[i + j + 1] += a[i] * b[j];
      }
    }
  }

  polynomial eliminate(std::vector<std::vector<polynomial>> wgt) const {
    if (wgt.empty()) {
      polynomial res(nd);
      res[0] = T(1);
      return res;
    }

    std::vector<polynomial> las = std::move(wgt.back());
    wgt.pop_back();
    std::vector<polynomial> sl = std::move(wgt.back());
    wgt.pop_back();
    int rem = int(wgt.size());

    polynomial res = eliminate(wgt);
    for (T &cf : res) {
      cf = -cf;
    }

    for (int a = 0; a < rem; a++) {
      for (int b = 0; b < a; b++) {
        multiply_and_shift_add(wgt[a][b], las[a], sl[b]);
        multiply_and_shift_add(wgt[a][b], sl[a], las[b]);
      }
    }
    polynomial lin = eliminate(std::move(wgt));
    multiply_and_shift_add(res, las[rem], lin);
    for (int deg = 0; deg < nd; deg++) {
      res[deg] += lin[deg];
    }
    return res;
  }
};

} // 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>> &mat) {
  int sz = int(mat.size());
  assert(sz % 2 == 0);
  for (const auto &row : mat) {
    assert(int(row.size()) == sz);
  }
  int nd = sz / 2 + 1;
  std::vector<std::vector<std::vector<T>>> wgt(sz);
  for (int row = 0; row < sz; row++) {
    wgt[row].assign(row, std::vector<T>(nd));
    for (int col = 0; col < row; col++) {
      assert(mat[row][col] == mat[col][row]);
      wgt[row][col][0] = mat[row][col];
    }
  }
  return hafnian_internal::calculator<T>(nd).eliminate(std::move(wgt)).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 nd;

  explicit calculator(int nd_) : nd(nd_) {}

  void multiply_and_shift_add(polynomial &res, const polynomial &a,
                              const polynomial &b) const {
    for (int i = 0; i < nd; i++) {
      for (int j = 0; i + j + 1 < nd; j++) {
        res[i + j + 1] += a[i] * b[j];
      }
    }
  }

  polynomial eliminate(std::vector<std::vector<polynomial>> wgt) const {
    if (wgt.empty()) {
      polynomial res(nd);
      res[0] = T(1);
      return res;
    }

    std::vector<polynomial> las = std::move(wgt.back());
    wgt.pop_back();
    std::vector<polynomial> sl = std::move(wgt.back());
    wgt.pop_back();
    int rem = int(wgt.size());

    polynomial res = eliminate(wgt);
    for (T &cf : res) {
      cf = -cf;
    }

    for (int a = 0; a < rem; a++) {
      for (int b = 0; b < a; b++) {
        multiply_and_shift_add(wgt[a][b], las[a], sl[b]);
        multiply_and_shift_add(wgt[a][b], sl[a], las[b]);
      }
    }
    polynomial lin = eliminate(std::move(wgt));
    multiply_and_shift_add(res, las[rem], lin);
    for (int deg = 0; deg < nd; deg++) {
      res[deg] += lin[deg];
    }
    return res;
  }
};

} // 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>> &mat) {
  int sz = int(mat.size());
  assert(sz % 2 == 0);
  for (const auto &row : mat) {
    assert(int(row.size()) == sz);
  }
  int nd = sz / 2 + 1;
  std::vector<std::vector<std::vector<T>>> wgt(sz);
  for (int row = 0; row < sz; row++) {
    wgt[row].assign(row, std::vector<T>(nd));
    for (int col = 0; col < row; col++) {
      assert(mat[row][col] == mat[col][row]);
      wgt[row][col][0] = mat[row][col];
    }
  }
  return hafnian_internal::calculator<T>(nd).eliminate(std::move(wgt)).back();
}

} // namespace noya