AtCoder Beginner Contest 336 G. 16 Integers（图计数 欧拉路径转欧拉回路 矩阵树定理 best定理）

题目

给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;
}