library

This documentation is automatically generated by online-judge-tools/verification-helper

Project maintained by tko919 Hosted on GitHub Pages — Theme by mattgraham

:warning: Pfaffian
(Math/pfaffian.hpp)

Depends on

Code

#pragma once
#include "Math/matrix.hpp"

template <typename T> T Pfaffian(Matrix<T> &a) {
    int n = a.h;
    assert(n == a.w and n % 2 == 0);
    T ret = 1;
    for (int i = 1; i < n; i += 2) {
        rep(j, i, n) if (a[j][i - 1] != 0) {
            swap(a[j], a[i]);
            for (auto &v : a.val)
                swap(v[j], v[i]);
            if (i != j)
                ret = -ret;
            break;
        }
        if (a[i][i - 1] == 0)
            return 0;
        ret *= -a[i][i - 1];
        T inv = T(1) / a[i][i - 1];
        rep(j, i + 1, n) {
            T z = a[j][i - 1] * inv;
            rep(k, 0, n) a[j][k] -= z * a[i][k];
        }
        inv = T(1) / a[i - 1][i];
        rep(j, i + 1, n) {
            T z = a[i - 1][j] * inv;
            rep(k, 0, n) a[k][j] -= z * a[k][i];
        }
    }
    return ret;
}

/**
 * @brief Pfaffian
 */
#line 2 "Math/matrix.hpp"

template <class T> struct Matrix {
    int h, w;
    vector<vector<T>> val;
    T det;
    Matrix() {}
    Matrix(int n) : h(n), w(n), val(vector<vector<T>>(n, vector<T>(n))) {}
    Matrix(int n, int m)
        : h(n), w(m), val(vector<vector<T>>(n, vector<T>(m))) {}
    vector<T> &operator[](const int i) {
        return val[i];
    }
    Matrix &operator+=(const Matrix &m) {
        assert(h == m.h and w == m.w);
        rep(i, 0, h) rep(j, 0, w) val[i][j] += m.val[i][j];
        return *this;
    }
    Matrix &operator-=(const Matrix &m) {
        assert(h == m.h and w == m.w);
        rep(i, 0, h) rep(j, 0, w) val[i][j] -= m.val[i][j];
        return *this;
    }
    Matrix &operator*=(const Matrix &m) {
        assert(w == m.h);
        Matrix<T> res(h, m.w);
        rep(i, 0, h) rep(j, 0, m.w) rep(k, 0, w) res.val[i][j] +=
            val[i][k] * m.val[k][j];
        *this = res;
        return *this;
    }
    Matrix operator+(const Matrix &m) const {
        return Matrix(*this) += m;
    }
    Matrix operator-(const Matrix &m) const {
        return Matrix(*this) -= m;
    }
    Matrix operator*(const Matrix &m) const {
        return Matrix(*this) *= m;
    }
    Matrix pow(ll k) {
        Matrix<T> res(h, h), c = *this;
        rep(i, 0, h) res.val[i][i] = 1;
        while (k) {
            if (k & 1)
                res *= c;
            c *= c;
            k >>= 1;
        }
        return res;
    }
    vector<int> gauss(int c = -1) {
        det = 1;
        if (val.empty())
            return {};
        if (c == -1)
            c = w;
        int cur = 0;
        vector<int> res;
        rep(i, 0, c) {
            if (cur == h)
                break;
            rep(j, cur, h) if (val[j][i] != 0) {
                swap(val[cur], val[j]);
                if (cur != j)
                    det *= -1;
                break;
            }
            det *= val[cur][i];
            if (val[cur][i] == 0)
                continue;
            rep(j, 0, h) if (j != cur) {
                T z = val[j][i] / val[cur][i];
                rep(k, i, w) val[j][k] -= val[cur][k] * z;
            }
            res.push_back(i);
            cur++;
        }
        return res;
    }
    Matrix inv() {
        assert(h == w);
        Matrix base(h, h * 2), res(h, h);
        rep(i, 0, h) rep(j, 0, h) base[i][j] = val[i][j];
        rep(i, 0, h) base[i][h + i] = 1;
        base.gauss(h);
        det = base.det;
        rep(i, 0, h) rep(j, 0, h) res[i][j] = base[i][h + j] / base[i][i];
        return res;
    }
    bool operator==(const Matrix &m) {
        assert(h == m.h and w == m.w);
        rep(i, 0, h) rep(j, 0, w) if (val[i][j] != m.val[i][j]) return false;
        return true;
    }
    bool operator!=(const Matrix &m) {
        assert(h == m.h and w == m.w);
        rep(i, 0, h) rep(j, 0, w) if (val[i][j] == m.val[i][j]) return false;
        return true;
    }
    friend istream &operator>>(istream &is, Matrix &m) {
        rep(i, 0, m.h) rep(j, 0, m.w) is >> m[i][j];
        return is;
    }
    friend ostream &operator<<(ostream &os, Matrix &m) {
        rep(i, 0, m.h) {
            rep(j, 0, m.w) os << m[i][j]
                              << (j == m.w - 1 and i != m.h - 1 ? '\n' : ' ');
        }
        return os;
    }
};

/**
 * @brief Matrix
 */
#line 3 "Math/pfaffian.hpp"

template <typename T> T Pfaffian(Matrix<T> &a) {
    int n = a.h;
    assert(n == a.w and n % 2 == 0);
    T ret = 1;
    for (int i = 1; i < n; i += 2) {
        rep(j, i, n) if (a[j][i - 1] != 0) {
            swap(a[j], a[i]);
            for (auto &v : a.val)
                swap(v[j], v[i]);
            if (i != j)
                ret = -ret;
            break;
        }
        if (a[i][i - 1] == 0)
            return 0;
        ret *= -a[i][i - 1];
        T inv = T(1) / a[i][i - 1];
        rep(j, i + 1, n) {
            T z = a[j][i - 1] * inv;
            rep(k, 0, n) a[j][k] -= z * a[i][k];
        }
        inv = T(1) / a[i - 1][i];
        rep(j, i + 1, n) {
            T z = a[i - 1][j] * inv;
            rep(k, 0, n) a[k][j] -= z * a[k][i];
        }
    }
    return ret;
}

/**
 * @brief Pfaffian
 */
Back to top page