推荐阅读： 矩阵树定理(+行列式) - command_block 的博客。
行列式
定义
这个东西一般用于求解图的生成树个数（矩阵树定理）。
称一个大小为 \(n\times n\) 的矩阵 \(A\) 的行列式为 \(\det(A)\)（或 \(|A|\)）。
\[\det(A)=\sum_{p\texttt{是一个大小为n排列}}(-1)^{F(p)}\prod_{i=1}^{n}A[i][p_i] \]其中 \(F(p)\) 表示排列 \(p\) 的逆序对数。\(\prod_{i=1}^{n}A[i][p_i]\) 可以理解为每一行选一个，每次选的列不与其他行重复，将选中的数乘起来。
性质
- 单位矩阵 \(I\) 的行列式为 \(1\)，上三角、下三角矩阵的行列式为对角线的值的乘积。
- 交换矩阵的任意两行，行列式变号。（交换一个排列的两个值，逆序对数的奇偶性一定发生变化）。
- 如果矩阵的某一行均为 \(0\) 则行列式为值为 \(0\)。（每行选且仅选一个）
- 若矩阵 \(A\) 的第 \(i\) 行和第 \(j\) 行相同，行列式的值为零。（交换第 \(i\)、\(j\) 行后行列式变号，但实际上矩阵 \(A\) 没有发生变化即行列式不变号，所以行列式值为零）。
- 若矩阵 \(A\) 的第 \(i\) 行乘以 \(k\)，行列式的值也乘以 \(k\)。（每一行都会选一个）
- \( \det(\begin{matrix}a+e &b+f\\c &d\\\end{matrix})= \det(\begin{matrix}a &b\\c &d\\\end{matrix})+ \det(\begin{matrix}e &f\\c &d\\\end{matrix}) \)（结合行列式定义自行理解）
- \( \det(\begin{matrix}a+kc &b+kd\\c &d\\\end{matrix})=\det(\begin{matrix}a &b\\c &d\\\end{matrix}) \)（\(\det(\begin{matrix}kc &kd\\c &d\\\end{matrix})=\det(\begin{matrix}c &d\\c &d\\\end{matrix})=0\)）
高斯消元求解行列式
根据性质 \(6\)，我们可以用高斯消元求解行列式，将原矩阵转化为上三角矩阵、顺便记录交换行的次数的奇偶性即可。
至于如何用第 \(i\) 行将 \(a_{j,i}\) 变为零（\(n\ge j>i\)）。可以用辗转相除法。
这样既可以求出任意一个矩阵在模 \(p\) 下的行列式的值了（\(p\) 可以不为质数）。
Code:
typedef long long LL;
#define arr2D vector<vector<int>>
#define arr vector<int>
int solve(arr2D &a,int n,int mo){
bool Swap=false; LL ans=1;
for(int i=1;i<=n;i++){
int tmp=i;
for(int j=i;j<=n;j++)if(a[j][i]){tmp=j;break;}
if(!a[tmp][i]){return 0;}
for(int j=tmp+1;j<=n;j++)if(a[j][i]&&a[j][i]<tmp)tmp=j;
if(i!=tmp){a[i].swap(a[tmp]);Swap^=1;}
for(int j=i+1;j<=n;j++){
if(a[j][i]>a[i][i]){a[i].swap(a[j]),Swap^=1;}
while(a[j][i]){
int d=a[i][i]/a[j][i];
for(int k=i;k<=n;k++)a[i][k]=(a[i][k]+(LL)(mo-d)*a[j][k]%mo)%mo;
a[i].swap(a[j]),Swap^=1;
}
}
ans=(LL)1LL*ans*a[i][i]%mo;
}
if(Swap)ans=(mo-ans)%mo;
return ans;
}
矩阵树定理
定义
一般用于求解图的生成树问题。（若无特殊说明，默认图为无向图，图中一定没有自环）
度数矩阵 \(D\) 满足 \(d_{i,i}=du_i\) 其他的 \(d_{i,j}=0,i\not=j\)（在无向图中 \(du_i\) 表示与 \(i\) 相连的边数），邻接矩阵 \(A\) 表示连边情况。其基尔霍夫矩阵为 \(K=D-A\)。则该图的生成树个数为矩阵 \(K\) 删去第 \(r\) 行、\(r\) 列后的矩阵的行列式的值。
如果矩阵树定理只能求生成树个数就太废了，实际上矩阵树定理求的所有生成树的边权乘积之和（求生成树个数时所有边权均为 \(1\)）。
因此将度数矩阵 \(D\) 中的 \(d_{i,i}\) 改为与 \(i\) 相连的边权之和，邻接矩阵 \(A\) 中的 \(a_{i,j}\) 改为 \(i\) 和 \(j\) 之间的边权之和。其基尔霍夫矩阵为 \(K=D-A\)。则该图的所有生成树的边权乘积之和为矩阵 \(K\) 删去第 \(r\) 行、\(r\) 列后的矩阵的行列式的值。
因此在用矩阵树定理求解生成树个数时可以支持重边（边 \((i,j)\) 的权值为 \((i,j)\) 在原图中出现的次数）。
有向图的矩阵树定理
邻接矩阵 \(A\) 为有向图的连边情况。
度数矩阵要分类讨论：
- 度数矩阵 \(D\) 的 \(d_{i,i}=\sum{j=1}^{n}A_{i,j}\) 时求的是外向树生成树个数。
- 度数矩阵 \(D\) 的 \(d_{i,i}=\sum{j=1}^{n}A_{j,i}\) 时求的是内向树生成树个数。
若以 \(k\) 为根则其生成树个数为基尔霍夫矩阵删去第 \(k\) 行、\(k\) 列的矩阵的行列式的值。
例题
直接容斥计算即可，\(\sum_{i}^{}(-1)^{n-P(i)}F(i)\)，\(P(i)\) 表示集合 \(i\) 的元素个数，\(F(i)\) 表示包含集合 \(i\) 内的所有建筑商可以修建的边的图的生成树个数。
#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef vector<int> arr;
typedef vector<arr> arr2D;
constexpr int mo=1e9+7;
int sol(arr2D a,int n){
LL ans=1;bool Swap=false;
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
a[i][j]%=mo;
for(int i=1;i<=n;i++){
int tmp=i;
for(int j=i;j<=n;j++)if(a[j][i]){tmp=j;break;}
if(!a[tmp][i])return 0;
for(int j=tmp+1;j<=n;j++)if(a[j][i]&&a[j][i]<a[tmp][i])tmp=j;
if(i!=tmp)a[i].swap(a[tmp]),Swap^=1;
for(int j=i+1;j<=n;j++){
if(a[j][i]<a[i][i])a[i].swap(a[j]),Swap^=1;
while(a[j][i]){
int d=a[i][i]/a[j][i];
for(int k=i;k<=n;k++)
a[i][k]=(a[i][k]+(LL)1LL*(mo-d)*a[j][k])%mo;
a[i].swap(a[j]),Swap^=1;
}
}
ans=(LL)ans*a[i][i]%mo;
}
if(Swap)ans=(mo-ans)%mo;
return ans;
}
arr2D Val;
vector<pair<int,int>>e[55];
int n;
int Sol(int S){
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
Val[i][j]=0;
int t=0;
for(int i=0;i<n-1;i++)if((S>>i)&1){
t++;
for(auto j:e[i])
Val[j.first][j.first]++,
Val[j.second][j.second]++,
Val[j.first][j.second]--,
Val[j.second][j.first]--;
}
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
Val[i][j]=(mo+Val[i][j]%mo)%mo;
int ans=sol(Val,n-1);
if((t&1)!=((n-1)&1))ans=(mo-ans)%mo;
return ans;
}
int main(){
cin>>n;
Val.resize(n+1,arr(n+1));
for(int i=0;i<n-1;i++){
int m;cin>>m;e[i].resize(m);
for(auto& j:e[i])
cin>>j.first>>j.second;
}
int ans=0;
for(int i=1;i<(1<<n-1);i++)
(ans+=Sol(i))%=mo;
cout<<ans<<'\n';
return 0;
}
边权改为 \(\frac{p_{i,j}}{1-p_{i,j}}\)，注意 \(p_{i,j}=1\) 的情况，手动加个误差即可。
#include <bits/stdc++.h>
using namespace std;
typedef vector<double> arr;
typedef vector<arr> arr2D;
constexpr double eps=1e-10;
double sol(arr2D &a,int n){
double ans=1.;
for(int i=1;i<=n;i++){
for(int j=i+1;j<=n;j++)if(fabs(a[j][i])>eps){
double d=a[j][i]/a[i][i];
for(int k=i;k<=n;k++)a[j][k]-=d*a[i][k];
}
ans*=a[i][i];
}
return ans;
}
arr2D a;
int n;
int main(){
cin>>n;
a.resize(n+1,arr(n+1,0.));
double pi=1.;
for(int i=1;i<=n;i++){
double sum=0;
for(int j=1;j<=n;j++){
cin>>a[i][j];a[i][j]+=eps;
if(i<j)pi*=(1.-a[i][j]);
if(i^j)a[i][j]=a[i][j]/(1.-a[i][j]);
sum+=a[i][j],a[i][j]=-a[i][j];
}
a[i][i]=sum;
}
printf("%.6lf\n",sol(a,n-1)*pi);
return 0;
}