给出两幅 \(n(\leq 400)\) 个点的无向图 \(G_1 ,G_2\)，对于 \(G_1\) 的每一颗生成树，它的权值定义为有多少条边在 \(G_2\) 中出现。求 \(G_1\) 所有生成树的权值和。

Solution

很容易想到，设 \(G_1\) 中每条边的权值，这条边在 \(G_2\) 中出现则权值为 \(1\)，否则权值为 \(0\)。

现在就真的是求所有生成树的边权和的权值和了。

然而标准的 Matrix-Tree Theorem 求的是生成树的边权积的和。

现在我们定义每条只出现在 \(G_1\) 中的边边权为 \(1\)，同时出现在 \(G_1,G_2\) 中的边权为 \(x\)，则基尔霍夫矩阵的每个元素嗾使一个多项式，记为 \(B(x)\)。

\(\det B(x)\) 是一个 \(n-1\) 次多项式 \(f(x) = \sum a_i x^i\)，那么其中 \(a_i\) 就是使用了 \(i\) 条公共变的生成树个数。

于是答案就是 \(f'(1)=\sum ia_i=\det B(1) \cdot \sum_i \sum_j (B^{-1}(1))_{i,j}\cdot B'(1)_{i,j}\)

于是用 Gauss 消元法求行列式和逆矩阵即可

#include <bits/stdc++.h>

using namespace std;

#define int long long

const int N = 405;

const int mod = 998244353;

namespace mat {

int f[N][N<<1],a[N][N],n;

inline void exgcd(int a,int b,int &x,int &y) {

    if(!b) {

        x=1,y=0;

        return;

    }

    exgcd(b,a%b,x,y);

    int t=x;

    x=y,y=t-(a/b)*y;

}

inline int inv(int a,int b) {

    int x,y;

    return exgcd(a,b,x,y),(x%b+b)%b;

}

int getdet() {

    int det=1;

    int flag=0;

    for(int i=1; i<=n; i++) {

        for(int j=i+1; j<=n; j++) {

            int x=i,y=j;

            while(a[y][i]!=0) {

                int t=a[x][i]*inv(a[y][i],mod)%mod;

                for(int k=i; k<=n; k++) (a[x][k]-=t*a[y][k]%mod)%=mod;

                swap(x,y);

            }

            if(x!=i) {

                for(int k=1; k<=n; k++) {

                    swap(a[x][k],a[i][k]);

                }

                flag^=1;

            }

        }

        if(a[i][i]==0)  return 0;

        det=det*a[i][i]%mod;

    }

    if(flag) det=-det;

    det%=mod; det+=mod; det%=mod;

    return det;

}

int solve() {

    int m=n*2;

    for(int i=1;i<=n;i++) {

        for(int j=1;j<=n;j++) a[i][j]=f[i][j], f[i][j+n]=0;

    }

    int ret= getdet();

    for(int i=1; i<=n; ++i) {

        f[i][n+i]=1;

    }

    for(int i=1; i<=n; ++i) {

        for(int j=i; j<=n; j++)

            if(f[j][i]) {

                for(int k=1; k<=m; k++)

                    swap(f[i][k],f[j][k]);

                break;

            }

        if(!f[i][i]) {

            return 0;

        }

        int r=inv(f[i][i],mod);

        for(int j=i; j<=m; ++j)

            f[i][j]=f[i][j]*r%mod;

        for(int j=1; j<=n; ++j)

            if(j!=i) {

                r=f[j][i];

                for(int k=i; k<=m; ++k)

                    f[j][k]=(f[j][k]-r*f[i][k]%mod+mod)%mod;

            }

    }

    return ret;

}

}

int n,b[N][N],bd[N][N];

char g1[N][N],g2[N][N];

signed main() {

    cin>>n;

    for(int i=1;i<=n;i++) {

        cin>>g1[i]+1;

    }

    for(int i=1;i<=n;i++) {

        cin>>g2[i]+1;

    }

    for(int i=1;i<=n;i++) {

        for(int j=1;j<i;j++) {

            if(g1[i][j]=='1') b[i][j]=-1, b[j][i]=-1, b[i][i]++, b[j][j]++;

            if(g1[i][j]=='1' && g2[i][j]=='1') bd[i][j]=-1, bd[j][i]=-1, bd[i][i]++, bd[j][j]++;

        }

    }

    for(int i=1;i<n;i++) {

        for(int j=1;j<n;j++) {

            mat::f[i][j]=b[i][j];

        }

    }

    --n;

    mat::n=n;

    int det=mat::solve();

    int ans=0;

    /*for(int i=1;i<=n;i++) {

        for(int j=1;j<=n;j++) cout<<b[i][j]<<" ";

        cout<<endl;

    }

    for(int i=1;i<=n;i++) {

        for(int j=1;j<=n;j++) cout<<bd[i][j]<<" ";

        cout<<endl;

    }

    for(int i=1;i<=n;i++) {

        for(int j=1;j<=n;j++) cout<<mat::f[i][j+n]<<" ";

        cout<<endl;

    }*/

    for(int i=1;i<=n;i++) {

        for(int j=1;j<=n;j++) {

            ans+=mat::f[i][j+n]*bd[i][j];

            ans%=mod;

            ans+=mod;

            ans%=mod;

        }

    }

    //cout<<ans<<" "<<det<<endl;

    cout<<((ans*det)%mod+mod)%mod;

}