CSU 1805: Three Capitals(BEST定理)

$给你A,B,G三个点,AB之间有a条边，AG之间有b条边，BG之间有c条边$
$问你从A出发然后回到A，走过所有的边，欧拉回路的个数$
$求有向图的欧拉回路个数，是BEST定理$
$ec(G)=t_s(G)\cdot deg(s)! \cdot \prod_{v\in V,\ v\ne s} (deg(v)-1)!,\ t_s(G):=以s为根的外向树的个数$
$如何计算有向图的外向树个数，这要用到Matrix Tree定理$
$基尔霍夫矩阵K=度数矩阵D-邻接矩阵A$
$无向图的度数矩阵就是每个点自己的度数$
$有向图的度数矩阵就是每个点自己的入度$
$邻接矩阵是表示u->v边的个数的矩阵$
$重边：按照边数计算，自环：不计入度数$
$无向图生成树计数：c=|K的任意1个n-1阶主子式|$
$有向图外向树计数：c=|去掉根所在的那阶得到的主子式|$

---

$以上都是copyTaosama大佬的博客的$
$这题只要枚举A点，在a条边里的入边个数，因为欧拉回路，对于一个点：入度=出度$
$所以可以计算出其他点的入度和出度，然后构造基尔霍夫矩阵，然后套一下公式$

---

代码：

#include <map>
#include <set>
#include <stack>
#include <queue>
#include <cmath>
#include <string>
#include <vector>
#include <cstdio>
#include <cctype>
#include <cstring>
#include <sstream>
#include <cstdlib>
#include <iostream>
#include <algorithm>
#pragma comment(linker,"/STACK:102400000,102400000")

using namespace std;
#define   MAX           100010
#define   MAXN          1000005
#define   maxnode       5
#define   sigma_size    30
#define   lson          l,m,rt<<1
#define   rson          m+1,r,rt<<1|1
#define   lrt           rt<<1
#define   rrt           rt<<1|1
#define   middle        int m=(r+l)>>1
#define   LL            long long
#define   ull           unsigned long long
#define   mem(x,v)      memset(x,v,sizeof(x))
#define   lowbit(x)     (x&-x)
#define   pii           pair<int,int>
#define   bits(a)       __builtin_popcount(a)
#define   mk            make_pair
#define   limit         10000

//const int    prime = 999983;
const int    INF   = 0x3f3f3f3f;
const LL     INFF  = 0x3f3f;
const double pi    = acos(-1.0);
const double inf   = 1e18;
const double eps   = 1e-6;
const LL     mod   = 1e9+7;
const ull    mx    = 133333331;

/*****************************************************/
inline void RI(int &x) {
      char c;
      while((c=getchar())<'0' || c>'9');
      x=c-'0';
      while((c=getchar())>='0' && c<='9') x=(x<<3)+(x<<1)+c-'0';
 }
/*****************************************************/

bool isFreeX[5];
LL A[5][5];
LL gauss(int n, int m) {//求行列式
    for(int i = 0; i < m; ++i) isFreeX[i] = false;
    LL ret = 1, neg = 0;
    int r = 1, c = 1;//求n-1阶的行列式，去掉第一阶，所以从1开始
    for(; r < n && c < m; ++r, ++c) {
        int p = r;
        for(; p < n; ++p) if(A[p][c]) break;
        if(p == n) {--r; isFreeX[c] = true; continue;}
        if(p != r) {
            neg ^= 1;
            for(int i = c; i <= m; ++i) swap(A[p][i], A[r][i]);
        }

        //eliminate coefficient
        for(int i = r + 1; i < n; ++i) {
            while(A[i][c]) {
                LL delta = A[i][c] / A[r][c];
                for(int j = c; j <= m; ++j) {
                    A[i][j] += mod - delta * A[r][j] % mod;
                    A[i][j] %= mod;
                }
                if(!A[i][c]) break;
                neg ^= 1;
                for(int j = c; j <= m; ++j) swap(A[r][j], A[i][j]);
            }
        }
    }
    if(r != n) return 0;
    //0-r-1求n阶行列式，1-r-1求n-1阶行列式
    for(int i = 1; i < r; ++i) ret = ret * A[i][i] % mod;
    if(neg) ret = (-ret + mod) % mod;
    return ret;
}

bool check(int x,int y){
    if(x%2) return false;
    x/=2;
    if(x>=0&&x<=y) return true;
    return false;
}

int deg[5];
LL fact[MAX];

void init(){
    fact[0]=1;
    for(int i=1;i<=100000;i++) fact[i]=fact[i-1]*i%mod;
}

LL qpow(LL a,LL n){
    LL ans=1;
    while(n){
        if(n&1) ans=ans*a%mod;
        a=a*a%mod;
        n>>=1;
    }
    return ans;
}
LL C(int n,int m){
    return (fact[n]*qpow(fact[m],mod-2)%mod)*qpow(fact[n-m],mod-2)%mod;
}
int main(){
    //freopen("in.txt","r",stdin);
    int a,b,c;
    init();
    while(cin>>a>>b>>c){
        LL ans=0;
        for(int x=0;x<=a;x++){
            int y=2*x+c-a;
            if(!check(y,c)) continue;
            y/=2;
            int z=2*y+b-c;
            if(!check(z,b)) continue;
            z/=2;
            if(x+b-z!=a-x+z) continue;
            deg[0]=x+b-z;
            deg[1]=y+a-x;
            deg[2]=z+c-y;
            for(int i=0;i<3;i++) A[i][i]=deg[i];
            A[0][1]=-(a-x);A[0][2]=-z;
            A[1][0]=-x;A[1][2]=-(c-y);
            A[2][0]=-(b-z);A[2][1]=-y;
            LL ret=(C(a,x)*C(c,y)%mod)*C(b,z)%mod;
            ret=ret*gauss(3,3)%mod;
            for(int i=0;i<3;i++) ret=ret*fact[deg[i]-1]%mod;
            ret=ret*deg[0]%mod;
        //cout<<ret<<endl;
            ans+=ret;
            if(ans>=mod) ans-=mod;
        }   
        cout<<ans<<endl;
    }
    return 0;
}