题解 ABC336G【16 Integers】

萌萌 BEST 定理练习题。赛时几乎做出来了，但写挂了，现在在火车上没事干就给补了。

考虑建图，图中共有 \(8\) 个节点，节点的编号是 \((\mathbb{Z}/2\mathbb{Z})^3\) 的每个元素。对于每个四元组 \((i,j,k,l)\in(\mathbb{Z}/2\mathbb{Z})^4\)，在图中连 \(X_{i,j,k,l}\) 条 \((i,j,k)\to(j,k,l)\) 的有向边。显然，每个合法的数列 \(\{A_i\}\) 恰好一一对应了图中的一条欧拉通路。

欧拉通路即去掉一条边的欧拉回路，对于每种可能的去掉边的方案，我们都把它连上并统计欧拉回路个数即可。边是无序的，所以需要除以一些阶乘。

// 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)

//By: OIer rui_er
#include <bits/stdc++.h>
#define rep(x, y, z) for(int x = (y); x <= (z); ++x)
#define per(x, y, z) for(int x = (y); x >= (z); --x)
#define debug(format...) fprintf(stderr, format)
#define fileIO(s) do {freopen(s".in", "r", stdin); freopen(s".out", "w", stdout);} while(false)
#define endl '\n'
using namespace std;
typedef long long ll;

mt19937 rnd(std::chrono::duration_cast<std::chrono::nanoseconds>(std::chrono::system_clock::now().time_since_epoch()).count());
int randint(int L, int R) {
    uniform_int_distribution<int> dist(L, R);
    return dist(rnd);
}

template<typename T> void chkmin(T& x, T y) {if(x > y) x = y;}
template<typename T> void chkmax(T& x, T y) {if(x < y) x = y;}

template<int mod>
inline unsigned int down(unsigned int x) {
	return x >= mod ? x - mod : x;
}

template<int mod>
struct Modint {
	unsigned int x;
	Modint() = default;
	Modint(unsigned int x) : x(x) {}
	friend istream& operator>>(istream& in, Modint& a) {return in >> a.x;}
	friend ostream& operator<<(ostream& out, Modint a) {return out << a.x;}
	friend Modint operator+(Modint a, Modint b) {return down<mod>(a.x + b.x);}
	friend Modint operator-(Modint a, Modint b) {return down<mod>(a.x - b.x + mod);}
	friend Modint operator*(Modint a, Modint b) {return 1ULL * a.x * b.x % mod;}
	friend Modint operator/(Modint a, Modint b) {return a * ~b;}
	friend Modint operator^(Modint a, int b) {Modint ans = 1; for(; b; b >>= 1, a *= a) if(b & 1) ans *= a; return ans;}
	friend Modint operator~(Modint a) {return a ^ (mod - 2);}
	friend Modint operator-(Modint a) {return down<mod>(mod - a.x);}
	friend Modint& operator+=(Modint& a, Modint b) {return a = a + b;}
	friend Modint& operator-=(Modint& a, Modint b) {return a = a - b;}
	friend Modint& operator*=(Modint& a, Modint b) {return a = a * b;}
	friend Modint& operator/=(Modint& a, Modint b) {return a = a / b;}
	friend Modint& operator^=(Modint& a, int b) {return a = a ^ b;}
	friend Modint& operator++(Modint& a) {return a += 1;}
	friend Modint operator++(Modint& a, int) {Modint x = a; a += 1; return x;}
	friend Modint& operator--(Modint& a) {return a -= 1;}
	friend Modint operator--(Modint& a, int) {Modint x = a; a -= 1; return x;}
	friend bool operator==(Modint a, Modint b) {return a.x == b.x;}
	friend bool operator!=(Modint a, Modint b) {return !(a == b);}
};

const ll N = 1e6 + 5, mod = 998244353;

typedef Modint<mod> mint;

ll n, fac[N], a[16], q[8][8], ideg[8], odeg[8], vis[8], A[8][8], id[8];
vector<ll> e[8];

void dfs(int u) {
    vis[u] = 1;
    for(int v : e[u]) if(!vis[v]) dfs(v);
}

ll det(int m) {
	ll now = 1;
	rep(i, 1, m) {
		rep(j, i + 1, m) {
			while(A[i][i]) {
				ll mt = (mod - A[j][i] / A[i][i]) % mod;
				rep(k, i, m) A[j][k] = (A[j][k] + mt * A[i][k] % mod) % mod;
				swap(A[i], A[j]);
                now = mod - now;
			}
			swap(A[i], A[j]);
            now = mod - now;
		}
	}
	rep(i, 1, m) now = now * A[i][i] % mod;
	return now;
}


int main() {
    ios::sync_with_stdio(false);
    cin.tie(0); cout.tie(0);
    rep(S, 0, 15) {
        cin >> a[S];
        n += a[S];
    }
    fac[0] = 1;
    rep(i, 1, n) fac[i] = fac[i - 1] * i % mod;
    rep(S, 0, 15) {
        ll S1 = S >> 1, S2 = S & 7;
        odeg[S1] += a[S];
        ideg[S2] += a[S];
        q[S1][S2] += a[S];
        rep(i, 1, a[S]) e[S1].push_back(S2);
    }
    mint ans = 0;
    rep(a, 0, 7) {
        rep(b, 0, 7) {
            memset(vis, 0, sizeof(vis));
            memset(A, 0, sizeof(A));
            memset(id, 0, sizeof(id));
            ++odeg[b]; ++ideg[a];
            e[b].push_back(a);
            dfs(b);
            bool ok = true;
            rep(u, 0, 7) if(ideg[u] != odeg[u] || !vis[u] && ideg[u]) ok = false;
            if(ok) {
                ll m = 0;
                rep(u, 0, 7) id[u] = -1;
                rep(u, 0, 7) if(vis[u]) id[u] = m++;
                rep(u, 0, 7) {
                    if(vis[u]) {
                        A[id[u]][id[u]] = odeg[u];
                        for(ll v : e[u]) A[id[u]][id[v]] = (A[id[u]][id[v]] - 1 + mod) % mod;
                    }
                }
                mint now = det(m - 1);
                rep(u, 0, 7) if(vis[u] && odeg[u]) now *= fac[odeg[u] - 1];
                rep(u, 0, 7) rep(v, 0, 7) now /= fac[q[u][v]];
                ans += now;
            }
            --odeg[b]; --ideg[a];
            e[b].pop_back();
        }
    }
    cout << ans << endl;
    return 0;
}