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