先从上往下消成一个上三角矩阵
再从下往上代入求值
1.出现类似0=1形式:无解
2.出现0=0形式:多组解
3.其他情况:一组解
typedef double db[N][N];
int n;
bool Gauss(db A){
for(int i=1;i<n;i++){
int r=i;
for(int j=i+1;j<=n;j++)
if(fabs(A[j][i]-A[r][i])>eps) r=j; //为保证精度,选取最大的A[r][i]作除数
for(int j=i;j<=n+1;j++) swap(A[i][j],A[r][j]);
if(fabs(A[i][i])<=eps) continue; //可能多解也可能无解
for(int j=i+1;j<=n;j++){
double x=A[j][i]/A[i][i];
for(int k=i;k<=n+1;k++)
A[j][k]-=A[i][k]*x;
}
}
for(int i=n;i>0;i--){
for(int j=n;j>i;j--) A[i][n+1]-=A[i][j]*A[j][n+1];
if(fabs(A[i][i])<=eps) return false; //可能多解也可能无解
A[i][n+1]/=A[i][i];
}
return true;
}
对于 \(n \times n\) 的矩阵 \(A\) , 若存在矩阵 \(B\) 满足 \(A \times B = E\) (\(E\) 为单位矩阵),则称 \(B\) 为 \(A\) 的逆矩阵
令矩阵 \(P=[A|E]\) ,对 \(P\) 进行一些线性变换使之变成 \([E|B]\) 的形式,则 \(A \times B = E\)
首先有一个可证的性质,对矩阵 \(M\) 进行线性变换 等价于 \(M\) 左乘一个矩阵 \(N\) (证明略)
然后就比较显然了,\(A \times B = E,E \times B = B\)
先从上往下将 \(P\) 的 \(A\) 部分消成上三角矩阵
再从下往上将 \(P\) 的 \(A\) 部分消成单位矩阵
(对大质数取模版)
#define xzy 1000000007
typedef long long ll;
typedef int Mat[N][N];
int Pow_mod(int x,int y){
int ret=1;
while(y){
if(y&1) ret=((ll)ret*x)%xzy;
x=((ll)x*x)%xzy;
y>>=1;
}
return ret;
}
bool Gauss_inv(Mat A,int n,int m){ // m=2n
for(int i=1;i<=n;i++){
int r=i;
for(int j=i+1;j<=n;j++)
if(A[j][i]>A[r][i]) r=j;
for(int j=i;j<=m;j++) swap(A[i][j],A[r][j]);
if(!A[i][i]) return false;
int y=Pow_mod(A[i][i],xzy-2);
for(int j=i+1;j<=n;j++){
int x=(ll)A[j][i]*y%xzy;
for(int k=i;k<=m;k++)
A[j][k]=(A[j][k]-(ll)A[i][k]*x%xzy+xzy)%xzy;
}
}
for(int i=n;i>0;i--){
for(int j=i+1;j<=n;j++){
for(int k=n+1;k<=m;k++)
A[i][k]=(A[i][k]-(ll)A[i][j]*A[j][k]%xzy+xzy)%xzy;
A[i][j]=0;
}
int x=Pow_mod(A[i][i],xzy-2);
for(int j=n+1;j<=m;j++) A[i][j]=((ll)A[i][j]*x)%xzy;
A[i][i]=1;
}
return true;
}
略(会用就行了)
交换行列式的两行,行列式取相反数
将行列式某一行元素乘以同一个数后加到另一行对应的元素上去,行列式不变
上三角矩阵或下三角矩阵行列式为正对角线上元素乘积
(对合数取模——利用扩展欧几里得)
#define xzy 1000000000
typedef int Mat[N][N];
typedef long long ll;
int Gauss(Mat A,int n){
int ret=1;
for(int j=1;j<=n;j++){
for(int i=j+1;i<=n;i++)
while(A[i][j]){ //扩欧 gcd(a,b)= (b==0? a : gcd(b,a%b))
int t=A[j][j]/A[i][j]; //A[i][j]相当于b,A[j][j]相当于a
for(int k=j;k<=n;k++){
A[j][k]=(A[j][k]-(ll)A[i][k]*t%xzy+xzy)%xzy;
swap(A[i][k],A[j][k]);
}
ret*=-1;
}
ret=((ll)ret*A[j][j]%xzy+xzy)%xzy;
}
return ret;
}
度数矩阵 \(D\):
\[
D[i][j]= \begin{cases}
度数,&i==j \0,&i!=j
\end{cases}
\]
邻接矩阵 \(A\):
\[
A[i][j]= \begin{cases}
1,&edge(i,j)=true \0,&else
\end{cases}
\]
基尔霍夫(\(kirchhoff\))矩阵 \(K\): \(K=D-A\)
余子式:一个矩阵 \(C\) 的余子式 \(M[i,j]\) 表示 \(C\) 去掉第 \(i\) 行与第 \(j\) 列后得到的矩阵的行列式
无向图的生成树数量,等于基尔霍夫矩阵任意余子式 \(M[i,i]\) (!!!注意是 \("i,i"\))的行列式
如果图不连通,则其任意余子式 \(M[i,i]\) 的行列式均为0
(模板题:\(bzoj4031\) \([HEOI2015]\) 小 \(Z\) 的房间)
#include<cstdio>
#include<iostream>
#include<algorithm>
#define xzy 1000000000
using namespace std;
const int N = 85;
typedef int Mat[N][N];
typedef long long ll;
int Gauss(Mat A,int n){
int ret=1;
for(int j=1;j<=n;j++){
for(int i=j+1;i<=n;i++)
while(A[i][j]){
int t=A[j][j]/A[i][j];
for(int k=j;k<=n;k++){
A[j][k]=(A[j][k]-(ll)A[i][k]*t%xzy+xzy)%xzy;
swap(A[i][k],A[j][k]);
}
ret*=-1;
}
ret=((ll)ret*A[j][j]%xzy+xzy)%xzy;
}
return ret;
}
Mat a;
char ch[10][10];
int dir[4][2]={-1,0,1,0,0,-1,0,1};
int tot,id[10][10];
int main()
{
int n,m;
scanf("%d%d",&n,&m);
for(int i=1;i<=n;i++) scanf("%s",ch[i]+1);
for(int i=1;i<=n;i++)
for(int j=1;j<=m;j++)
if(ch[i][j]=='.') id[i][j]=++tot;
for(int i=1;i<=n;i++)
for(int j=1;j<=m;j++){
if(ch[i][j]=='*') continue;
for(int k=0;k<4;k++)
if(id[i+dir[k][0]][j+dir[k][1]]){
a[id[i][j]][id[i][j]]++;
a[id[i][j]][id[i+dir[k][0]][j+dir[k][1]]]=-1;
}
}
printf("%d\n",Gauss(a,tot-1));
return 0;
}
入度矩阵 \(D\):
\[
D[i][j]= \begin{cases}
入度,&i==j \0,&i!=j
\end{cases}
\]
出度矩阵 \(D\):
\[
D[i][j]= \begin{cases}
出度,&i==j \0,&i!=j
\end{cases}
\]
邻接矩阵 \(A\):同上
内向生成树:生成树上的边为儿子指向父亲
外向生成树:生成树上的边为父亲指向儿子
扩展基尔霍夫矩阵\(K\): \(K=D-A\)
以 \(root\) 为根的生成树个数为 \(K\) 的余子式 \(M[root,root]\) 的行列式
其中:
若为内向生成树,则 \(D\) 为出度矩阵
若为外向生成树,则 \(D\) 为入度矩阵
(模板题:\(bzoj5297\) \([Cqoi2018]\) 社交网络)
#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cmath>
#define xzy 10007
using namespace std;
const int N = 255;
typedef long long ll;
typedef int Mat[N][N];
int Gauss(Mat A,int n){
int ret=1;
for(int j=1;j<=n;j++){
for(int i=j+1;i<=n;i++)
while(A[i][j]){
int t=A[j][j]/A[i][j];
for(int k=j;k<=n;k++){
A[j][k]=(A[j][k]-(ll)A[i][k]*t%xzy+xzy)%xzy;
swap(A[i][k],A[j][k]);
}
ret*=-1;
}
ret=((ll)ret*A[j][j]%xzy+xzy)%xzy;
}
return ret;
}
Mat a;
int main()
{
int n,m,x,y;
scanf("%d%d",&n,&m);
for(int i=0;i<m;i++){
scanf("%d%d",&x,&y);//y->x
x--; y--;
a[x][x]++; a[y][x]--;
}
printf("%d\n",Gauss(a,n-1));
return 0;
}
原文:https://www.cnblogs.com/lindalee/p/10702908.html