Tifa's CP Library

:heavy_check_mark: src/lalg/mat/pfaffian/lib.hpp

Depends on

Verified with

Code

#pragma once

#include "../../../util/alias/num/lib.hpp"

namespace tifa_libs {

template <class Mat, class Is0>
requires requires(Is0 is0, TPN Mat::val_t t) {
  { is0(t) } -> std::same_as<bool>;
}
CEXP auto pfaffian(Mat mat, Is0 is0) NE {
  using T = TPN Mat::val_t;
  cu32 n = mat.row();
  assert(n == mat.col() && !(n & 1));
  bool neg = false;
  for (u32 i = 0; i < n; i += 2) {
    u32 p = i + 1;
    flt_ (u32, j, i + 2, n)
      if (!is0(mat(i, j))) {
        p = j;
        break;
      }
    if (p != i + 1) {
      neg ^= 1;
      flt_ (u32, j, 0, i + 1) swap(mat(j, i + 1), mat(j, p));
      flt_ (u32, j, i + 2, p) {
        swap(mat(i + 1, j), mat(j, p));
        mat(i + 1, j) = -mat(i + 1, j), mat(j, p) = -mat(j, p);
      }
      mat(i + 1, p) = -mat(i + 1, p);
      flt_ (u32, j, p + 1, n) swap(mat(i + 1, j), mat(p, j));
    }
    if (is0(mat(i, i + 1))) return T{};
    T const _ = mat(i, i + 1).inv();
    flt_ (u32, j, i + 2, n)
      if (mat(i, j)) {
        T const _2 = -mat(i, j) * _;
        flt_ (u32, k, 0, i + 1) mat(k, j) += _2 * mat(k, i + 1);
        flt_ (u32, k, i + 2, j) mat(k, j) -= _2 * mat(i + 1, k);
        flt_ (u32, k, j + 1, n) mat(j, k) += _2 * mat(i + 1, k);
      }
  }
  T ans = 1;
  for (u32 i = 0; i < n; i += 2) ans *= mat(i, i + 1);
  retif_((neg), -ans, ans);
}

}  // namespace tifa_libs
#line 2 "src/lalg/mat/pfaffian/lib.hpp"

#line 2 "src/util/alias/num/lib.hpp"

#line 2 "src/util/util/lib.hpp"
// https://github.com/Tiphereth-A/CP-lib
#include <bits/extc++.h>
// clang-format off
namespace tifa_libs {

#define CEXP constexpr
#define CEXPE constexpr explicit
#define CR const&
#define CP const*
#define PC *const
#define CPC const*const
#define TPN typename
#define NE noexcept
#define CNE const noexcept
#define ND [[nodiscard]]
#define cT_(...) std::conditional_t<sizeof(__VA_ARGS__) <= sizeof(size_t) * 2, __VA_ARGS__, __VA_ARGS__ CR>
// NOLINTNEXTLINE(misc-const-correctness)
#define flt_(T, i, l, r, ...) for (T i = (l), i##e = (r)__VA_OPT__(, ) __VA_ARGS__; i < i##e; ++i)
#define retif_(cond, if_true, ...) if cond return if_true __VA_OPT__(; else return __VA_ARGS__)
#ifdef ONLINE_JUDGE
#undef assert
#define assert(x) 42
#endif

using namespace std::ranges;
using namespace std::literals;

template <class T>
CEXP T abs(T x) NE { retif_((x < 0), -x, x); }

}  // namespace tifa_libs
// clang-format on
#line 4 "src/util/alias/num/lib.hpp"
// clang-format off
namespace tifa_libs {

#define mk0_(w, t) using w = t; using c##w = const t
#define mk_(w, t) mk0_(w, t); CEXP w operator""_##w(unsigned long long x) NE { return (w)x; }
mk_(i8, int8_t) mk_(u8, uint8_t) mk_(i16, int16_t) mk_(u16, uint16_t) mk_(i32, int32_t) mk_(u32, uint32_t) mk_(i64, int64_t) mk_(u64, uint64_t) mk_(isz, ssize_t) mk_(usz, size_t) mk_(chr, char) mk_(schr, signed char) mk_(uchr, unsigned char) mk_(sint, signed) mk_(uint, unsigned);
mk0_(i128, __int128_t); mk0_(u128, __uint128_t); mk0_(f32, float); mk0_(f64, double); mk0_(f128, long double);
#undef mk0_
#undef mk_

}  // namespace tifa_libs
// clang-format on
#line 4 "src/lalg/mat/pfaffian/lib.hpp"

namespace tifa_libs {

template <class Mat, class Is0>
requires requires(Is0 is0, TPN Mat::val_t t) {
  { is0(t) } -> std::same_as<bool>;
}
CEXP auto pfaffian(Mat mat, Is0 is0) NE {
  using T = TPN Mat::val_t;
  cu32 n = mat.row();
  assert(n == mat.col() && !(n & 1));
  bool neg = false;
  for (u32 i = 0; i < n; i += 2) {
    u32 p = i + 1;
    flt_ (u32, j, i + 2, n)
      if (!is0(mat(i, j))) {
        p = j;
        break;
      }
    if (p != i + 1) {
      neg ^= 1;
      flt_ (u32, j, 0, i + 1) swap(mat(j, i + 1), mat(j, p));
      flt_ (u32, j, i + 2, p) {
        swap(mat(i + 1, j), mat(j, p));
        mat(i + 1, j) = -mat(i + 1, j), mat(j, p) = -mat(j, p);
      }
      mat(i + 1, p) = -mat(i + 1, p);
      flt_ (u32, j, p + 1, n) swap(mat(i + 1, j), mat(p, j));
    }
    if (is0(mat(i, i + 1))) return T{};
    T const _ = mat(i, i + 1).inv();
    flt_ (u32, j, i + 2, n)
      if (mat(i, j)) {
        T const _2 = -mat(i, j) * _;
        flt_ (u32, k, 0, i + 1) mat(k, j) += _2 * mat(k, i + 1);
        flt_ (u32, k, i + 2, j) mat(k, j) -= _2 * mat(i + 1, k);
        flt_ (u32, k, j + 1, n) mat(j, k) += _2 * mat(i + 1, k);
      }
  }
  T ans = 1;
  for (u32 i = 0; i < n; i += 2) ans *= mat(i, i + 1);
  retif_((neg), -ans, ans);
}

}  // namespace tifa_libs
Back to top page