题目
给16个非负整数，x[i∈(0,1)][j∈(0,1)][k∈(0,1)][l∈(0,1)]
求长为n+3的01串的方案数，满足长度为4的ijkl（2*2*2*2，16种情况）串恰为x[i][j][k][l]个
答案对998244353取模
思路来源
https://www.cnblogs.com/tzcwk/p/matrix-tree-best-theroem.html
矩阵树定理 - OI Wiki
知识点总结
矩阵树定理
对于一张无向图G，
记D为其度数矩阵，满足：
1. D[i][i]=i的度数
2. D[i][j]=0(i≠j）
记A为其邻接矩阵，满足：
1. A[i][j]=i与j之间边的条数，如果有重边则算作多条边。
记基尔霍夫矩阵（拉普拉斯矩阵）K=D−A，
则去掉第k行第k列得到的矩阵行列式即为G生成树的个数。
BEST 定理
对于一个欧拉图（有向图）而言，
从x出发回到x的欧拉回路的个数为：
其中，为任一个点的根向树形图的数量，用矩阵树定理求得，
后半部分，即对于每个点，再乘以每个点的度数减1的阶乘
根向树形图
根向树形图是一棵树，所有边都往根的方向指
一些技巧
求欧拉路径的数量
二层枚举欧拉路径的起点、终点，钦定加一条终点到起点的边，转化为求欧拉回路
仅要求所有边经过一次（部分点可以孤立）
忽略掉孤立点后，对剩下的点离散化后，重新建矩阵，求矩阵的秩
需要视题目决定是否需要乘deg[v]（从v出发回到v）
少乘度数的是带循环同构的边序列
思路来源
官方题解
题解
1. 把形如abcd的出现次数，转化为abc->bcd有向边的边的条数，
转化为成边数后，即求欧拉路径条数，
枚举欧拉路径起点终点后，强制加一条边，转成欧拉回路后，套用best定理
2. 对于枚举起点为i，终点为j的情况，先强制加一条j->i的边，
统计每个点x的入度in[x]、出度out[x]
① 忽略孤立点，即in[x]=out[x]=0
② 若in[x]和out[x]不相等，则无解
③ 若x非孤立点，且与i不联通，则无解
否则，
④忽略孤立点，将非孤立点重新编号建边
⑤按定义构建基尔霍夫矩阵K=D-A
即a[k][k]+=in[k]（D矩阵），a[k][l]-=b[k][l]（A矩阵）
rep(k,0,7){ // 基尔霍夫矩阵
a[to[k]][to[k]]=in[k];
rep(l,0,7){
if(!to[l])continue;
a[to[k]][to[l]]-=b[k][l];
//printf("%d ",a[k][l]);
}
}⑥求矩阵的秩，得到生成树数量，依次乘以(deg[i]-1)!得到欧拉回路数量
⑦欧拉回路是一个环，每个环被统计一次，
固定新增的那条边在串最后指向串最前，即可唯一对应一个串
但是注意到，当abc->bcd有两条相同的边x1、x2时，二者的顺序会被欧拉回路视为不同的方案
而在序列中，abcd是一个唯一的序列，被重复计算
所以，需要除掉完全相同的边的顺序，类似可重集全排列的方案数
代码
// Problem: G - 16 Integers
// Contest: AtCoder - AtCoder Beginner Contest 336
// URL: https://atcoder.jp/contests/abc336/tasks/abc336_g
// Memory Limit: 1024 MB
// Time Limit: 2000 ms
//
// Powered by CP Editor (https://cpeditor.org)
#include<bits/stdc++.h>
using namespace std;
#define rep(i,a,b) for(int i=(a);i<=(b);++i)
#define per(i,a,b) for(int i=(a);i>=(b);--i)
typedef long long ll;
typedef double db;
typedef pair<ll,int> P;
#define fi first
#define se second
#define pb push_back
#define dbg(x) cerr<<(#x)<<":"<<x<<" ";
#define dbg2(x) cerr<<(#x)<<":"<<x<<endl;
#define SZ(a) (int)(a.size())
#define sci(a) scanf("%d",&(a))
#define scll(a) scanf("%lld",&(a))
#define pt(a) printf("%d",a);
#define pte(a) printf("%d\n",a)
#define ptlle(a) printf("%lld\n",a)
#define debug(...) fprintf(stderr, __VA_ARGS__)
const int N=1e6+10,M=16,K=8,mod=998244353;
int t,fac[N],x[M],b[K][K],in[K],out[K],par[K];
int find(int x){
return par[x]==x?x:par[x]=find(par[x]);
}
void merge(int x,int y){
x=find(x),y=find(y);
if(x==y)return;
par[y]=x;
}
int modpow(int x,int n,int mod){
int res=1;
for(;n;n>>=1,x=1ll*x*x%mod){
if(n&1)res=1ll*res*x%mod;
}
return res;
}
// 求解行列式时模数不是质数，没法求逆元，这时只能利用辗转消除法进行高斯消元
// 矩阵的秩
int Gauss(vector<vector<int> >&a,int n){
//printf("x:%d\n",x);
int ans=1;
//swap(a[0],a[x]);
for(int i=2;i<=n;i++){
for(int j=i;j<=n;j++){
if(!a[i][i] && a[j][i]){
swap(a[i],a[j]);
ans=mod-ans;
break;
}
}
ans=1ll*ans*a[i][i]%mod;
for(int j=1;j<=n;j++){
if(i==j || !a[j][i])continue;
int t=1ll*a[j][i]*modpow(a[i][i],mod-2,mod)%mod;
for(int k=1;k<=n;k++){
a[j][k]=(a[j][k]-1ll*t*a[i][k]%mod+mod)%mod;
}
}
}
return ans;
}
void sol(){
rep(i,0,15){
b[i>>1][i&7]=x[i];
}
int res=0;
rep(i,0,7){
rep(j,0,7){
b[j][i]++;//欧拉路径->欧拉回路
bool ok=1;
memset(in,0,sizeof in);
memset(out,0,sizeof out);
rep(k,0,7)par[k]=k;
rep(k,0,7){
rep(l,0,7){
in[l]+=b[k][l];
out[k]+=b[k][l];
if(b[k][l])merge(k,l);
//if(b[k][l])printf("i:%d j:%d k:%d l:%d b:%d\n",i,j,k,l,b[k][l]);
}
}
rep(k,0,7){
if(!in[k] && !out[k])continue;//孤立点，只是要求遍历所有边时，可忽略
if(in[k]!=out[k])ok=0;//判出入度
if(find(k)!=find(i))ok=0;//判连通
//printf("i:%d j:%d k:%d in:%d out:%d\n",i,j,k,in[k],out[k]);
}
if(!ok){
b[j][i]--;
continue;
}
int bs=1;//deg[i]
rep(k,0,7){
if(!in[k])continue;
bs=1ll*bs*fac[in[k]-1]%mod;//\prod (deg[i]-1)
}
vector<vector<int>>a(K+1,vector<int>(K+1,0));
vector<int>to(K+1,0);
int id=0;
rep(k,0,7){
if(in[k])to[k]=++id;
}
rep(k,0,7){ // 基尔霍夫矩阵
a[to[k]][to[k]]=in[k];
rep(l,0,7){
if(!to[l])continue;
a[to[k]][to[l]]-=b[k][l];
//printf("%d ",a[k][l]);
}
}
int rank=Gauss(a,id);//求基尔霍夫矩阵的秩
//printf("i:%d j:%d bs:%d rk:%d\n",i,j,bs,rank);
bs=1ll*bs*rank%mod;
rep(k,0,15){
bs=1ll*bs*modpow(fac[x[k]],mod-2,mod)%mod;//去重
}
res=(res+bs)%mod;
b[j][i]--;
}
}
printf("%d\n",res);
}
void init(){
fac[0]=1;
rep(i,1,N-1)fac[i]=1ll*fac[i-1]*i%mod;
}
int main(){
init();
rep(i,0,15)sci(x[i]);
sol();
return 0;
}
