HDU 6064 RXD and numbers(BEST定理)

/**
HDU 6064
题意：有一个序列A，每个值都在[1,m]之间，且A1=An=1，对每个1 <= x <= m，A中一定有一个位置出现过x， 现有一个矩阵D，
其中D中(x,y)的位置的值表示在序列中Ai=x, Ai+1=y中i的个数，现在只知道矩阵D，问A的有效序列个数，很显然n的值是D矩阵
中所有的值加起来在加1

思路：Ai=x,Ai+1=y，表示y是紧接着在x之后出现，联想到图上的一条有向路径x->y，走过一个节点之后紧接着走下一个节点，又有A1=An=1，也就是
说，从标号为1的点出发又回到了1这个点， 这里联想到有向图的欧拉回路，正好满足这个模型，换言之就是求这个模型建立起来之后以1为起点的
欧拉回路的数量， 可用best定理求解， 最后注意一些欧拉回路不存在的条件，一种是作为无向图时图不连通，一种是有向图中某个（某些）点的
入度不等于出度,还有就是计数的时候重边会计数多次，用乘法原理，最后的结果除去所有重边数量的阶乘的乘积
**/

#include<cstdio>
#include<cstring>
#include<cmath>
#include<vector>
#include<iostream>
#include<queue>
#include<map>
#include<stack>
#include<set>
#include<algorithm>
typedef long long ll;
const int maxn = 510;
const ll mod = 998244353;
using namespace std;

ll K[maxn][maxn], D[maxn][maxn];
ll C[maxn][maxn], F[maxn * maxn];
ll O[maxn];
int m, kase = 1;
int pa[maxn];

int findset(int x) { return x == pa[x] ? x : pa[x] = findset(pa[x]); }
void unit(int x, int y) { pa[findset(x)] = findset(y); }

ll qmod(ll x, ll n) {
    ll ans = 1;
    while(n) {
        if(n & 1) ans = ans * x % mod;
        x = x * x % mod;
        n >>= 1;
    }
    return ans;
}

ll Gauss(int n) {
    for(int i = 1; i < n; i++) {
        bool found = false;
        for(int j = i; j < n; j++) {
            if(!K[j][i]) continue;
            for(int k = i; k < n; k++) swap(K[i][k], K[j][k]);
            found = true; break;
        }
        if(!found) return 0;

        for(int j = i + 1; j < n; j++) {
            ll inv = qmod(K[i][i], mod - 2), x = K[j][i] * inv % mod;
            if(!K[j][i]) continue; K[j][i] = 0;
            for(int k = i + 1; k < n; k++) {
                K[j][k] = (K[j][k] - (K[i][k] * x % mod) + mod * 2) % mod;
            }
        }
    }
    ll ans = 1;
    for(int i = 1; i < n; i++) ans = (ans * K[i][i]) % mod;
    return ans;
}

int main() {
    F[0] = 1;
    for(ll i = 1; i < maxn * maxn; i++) F[i] = F[i - 1] * i % mod;
    while(scanf("%d", &m) != EOF) {
        memset(D, 0, sizeof D);
        memset(O, 0, sizeof O);
        ll ks = 1;
        for(int i = 0; i < m; i++) pa[i] = i;
        for(int i = 0; i < m; i++) {
            for(int j = 0; j < m; j++) {
                scanf("%lld", &C[i][j]);
                D[j][j] += C[i][j];
                O[i] += C[i][j];
                if(C[i][j]) unit(i, j);
                ks = ks * F[C[i][j]] % mod;
            }
        }
        for(int i = 0; i < m; i++) {
            for(int j = 0; j < m; j++)
                K[i][j] = D[i][j] - C[i][j];
        }
        ll ans = 1, dd = 1;
        for(int i = 0; i < m; i++) {
            if(D[i][i] != O[i]) dd = 0;
            if(i && findset(i - 1) != findset(i)) dd = 0;
        }
        if(!dd) { printf("Case #%d: %lld\n", kase++, 0); continue; }
        for(int i = 0; i < m; i++) ans = ans * F[D[i][i] - 1] % mod;
        ans = ans * D[0][0] % mod;
        ans = ans * Gauss(m) % mod;
        ll inv = qmod(ks, mod - 2);
        ans = ans * inv % mod;
        printf("Case #%d: %lld\n", kase++, ans);

    }
    return 0;
}