library
This documentation is automatically generated by online-judge-tools/verification-helper
Hafnian of matrix
(Math/hafnian.hpp)
- View this file on GitHub
- Last update: 2024-04-26 03:18:17+09:00
- Include:
#include "Math/hafnian.hpp"
Depends on
Verified with
Code
#pragma once
#include "Math/matrix.hpp"
#include "Convolution/subset.hpp"
template<typename T>T Hafnian(Matrix<T>& a){
int n=a.h;
assert(n%2==0);
n>>=1;
vector<T> cycle(1<<n);
vector dp(1<<n,vector<T>(n*2));
rep(base,0,n)dp[1<<base][base*2]=1;
rep(mask,0,1<<n){
int top=-1;
rep(i,0,n)if(mask>>i&1)top=i;
rep(base,0,n*2)if(mask>>(base>>1)&1){
T add=dp[mask][base];
rep(nxt,0,top)if(!(mask>>nxt&1)){
dp[mask|(1<<nxt)][nxt*2+1]+=add*a[base][nxt*2];
dp[mask|(1<<nxt)][nxt*2]+=add*a[base][nxt*2+1];
}
}
}
rep(mask,1,1<<n){
int top=-1;
rep(i,0,n)if(mask>>i&1)top=i;
rep(base,0,n*2)cycle[mask]+=dp[mask][base]*a[base][top*2+1];
}
SubsetConvolution<T> buf;
auto ret=buf.exp(cycle);
return ret[(1<<n)-1];
}
/**
* @brief Hafnian of matrix
*/#line 2 "Math/hafnian.hpp"
#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){
if(val.empty())return {};
if(c==-1)c=w;
int cur=0; vector<int> res; det=1;
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 2 "Convolution/subset.hpp"
template <typename T, int LG = 20> struct SubsetConvolution {
using POL = array<T, LG + 1>;
vector<int> bpc;
SubsetConvolution() : bpc(1 << LG) {
rep(i, 1, 1 << LG) bpc[i] = bpc[i - (i & -i)] + 1;
}
void zeta(vector<POL> &a) {
int n = topbit(SZ(a));
rep(d, 0, n) {
rep(i, 0, 1 << n) if (i >> d & 1) {
const int pc = bpc[i];
rep(j, 0, pc) a[i][j] += a[i ^ (1 << d)][j];
}
}
}
void mobius(vector<POL> &a) {
int n = topbit(SZ(a));
rep(d, 0, n) {
rep(i, 0, 1 << n) if (i >> d & 1) {
const int pc = bpc[i];
rep(j, pc, n + 1) a[i][j] -= a[i ^ (1 << d)][j];
}
}
}
vector<T> mult(vector<T> &a, vector<T> &b) {
assert(a.size() == b.size());
int n = SZ(a), m = topbit(n);
vector<POL> A(n), B(n);
rep(i, 0, n) {
A[i][bpc[i]] = a[i];
B[i][bpc[i]] = b[i];
}
zeta(A);
zeta(B);
rep(i, 0, n) {
POL c = {};
rep(j, 0, m + 1) rep(k, 0, m + 1 - j) c[j + k] += A[i][j] * B[i][k];
swap(A[i], c);
}
mobius(A);
vector<T> ret(n);
rep(i, 0, n) ret[i] = A[i][bpc[i]];
return ret;
}
vector<T> TransposeMult(vector<T> &a, vector<T> &b) {
auto ret = a;
reverse(ALL(ret));
ret = mult(ret, b);
reverse(ALL(ret));
return ret;
}
vector<T> exp(vector<T> &f) {
int n = topbit(SZ(f));
vector<T> ret(1 << n);
ret[0] = 1;
rep(i, 0, n) {
vector<T> a = {ret.begin(), ret.begin() + (1 << i)};
vector<T> b = {f.begin() + (1 << i), f.begin() + (2 << i)};
a = mult(a, b);
copy(ALL(a), ret.begin() + (1 << i));
}
return ret;
}
vector<T> CompositionEGF(vector<T> &s, vector<T> &f) {
int n = topbit(SZ(s));
f.resize(n + 1);
vector<T> dp(1);
dp[0] = f.back();
rrep(d, 0, n) {
vector<T> ndp(1 << (n - d));
ndp[0] = f[d];
rep(i, 0, n - d) {
vector<T> a = {dp.begin(), dp.begin() + (1 << i)};
vector<T> b = {s.begin() + (1 << i), s.begin() + (2 << i)};
a = mult(a, b);
copy(ALL(a), ndp.begin() + (1 << i));
}
swap(dp, ndp);
}
return dp;
}
vector<T> Composition(vector<T> &s, vector<T> &f) {
int n = topbit(SZ(s));
T c = s[0];
s[0] = 0;
// taylor shift
vector<T> pw(n + 1), g(n + 1);
pw[0] = 1;
rep(i, 0, SZ(f)) {
rep(j, 0, n + 1) g[j] += f[i] * pw[j];
rrep(j, 0, n + 1) {
if (j != n)
pw[j + 1] += pw[j];
pw[j] *= c;
}
}
// to EGF
T fact = 1;
rep(i, 1, n + 1) {
fact *= i;
g[i] *= fact;
}
return CompositionEGF(s, g);
}
vector<T> TransposeCompositionEGF(vector<T> &s, vector<T> &t) {
int n = topbit(SZ(s));
vector<T> dp = t, ret(n + 1);
ret[0] = dp[0];
rep(d, 0, n) {
vector<T> ndp(1 << (n - d - 1), 0);
rrep(i, 0, n - d) {
vector<T> a = {dp.begin() + (1 << i), dp.begin() + (2 << i)};
vector<T> b = {s.begin() + (1 << i), s.begin() + (2 << i)};
a = TransposeMult(a, b);
rep(k, 0, SZ(a)) ndp[k] += a[k];
}
swap(dp, ndp);
ret[d + 1] = dp[0];
}
return ret;
}
vector<T> TransposeComposition(vector<T> &s, vector<T> &t, int m) {
int n = topbit(SZ(s));
T c = s[0];
s[0] = 0;
auto g = TransposeCompositionEGF(s, t);
vector<T> coe(m);
T pw = 1;
rep(i, 0, m) {
coe[i] = pw * Fact<T>(i, 1);
pw *= c;
}
vector<T> f(m);
rep(i, 0, m) rep(j, 0, n + 1) if (i + j < m) {
f[i + j] += coe[i] * g[j];
}
rep(i, 0, m) f[i] *= Fact<T>(i);
return f;
}
};
/**
* @brief Subset Convolution
*/
#line 5 "Math/hafnian.hpp"
template<typename T>T Hafnian(Matrix<T>& a){
int n=a.h;
assert(n%2==0);
n>>=1;
vector<T> cycle(1<<n);
vector dp(1<<n,vector<T>(n*2));
rep(base,0,n)dp[1<<base][base*2]=1;
rep(mask,0,1<<n){
int top=-1;
rep(i,0,n)if(mask>>i&1)top=i;
rep(base,0,n*2)if(mask>>(base>>1)&1){
T add=dp[mask][base];
rep(nxt,0,top)if(!(mask>>nxt&1)){
dp[mask|(1<<nxt)][nxt*2+1]+=add*a[base][nxt*2];
dp[mask|(1<<nxt)][nxt*2]+=add*a[base][nxt*2+1];
}
}
}
rep(mask,1,1<<n){
int top=-1;
rep(i,0,n)if(mask>>i&1)top=i;
rep(base,0,n*2)cycle[mask]+=dp[mask][base]*a[base][top*2+1];
}
SubsetConvolution<T> buf;
auto ret=buf.exp(cycle);
return ret[(1<<n)-1];
}
/**
* @brief Hafnian of matrix
*/