基尔霍夫矩阵树定理

例题

很显然就是让你求一个有\(n+1\)个点的图以某种形式连接的生成树计数,所以可以直接跑生成树计数的模板也就是Kirchhoff矩阵树(基尔霍夫矩阵)
之前一直知道这东西可以跑最小生成树计数,但是没有系统学过.这次来系统学一下.

基尔霍夫矩阵树定理

定义介绍

无向图定义

设\(G\)是一个有\(n\)个顶点的无向图,定义度数矩阵\(D(G)\)为:

\[D[i][i](G)=deg(i),D[i][j](G)=0,i\neq j \]

设\(\# e(i,j)\)为点\(i\)与点\(j\)相连的边数,并定义邻接矩阵\(A\)为:

\[A[i][j](G)=A[j][i](G)=\# e(i,j),i \neq j \]

定义\(Laplace\)矩阵(亦称为\(Kirchhoff\)矩阵)\(L\)为:

\[L(G)=D(G)-A(G) \]

记图\(G\)的生成树个数为\(t(G)\)

有向图定义

设\(G\)是一个有\(n\)个顶点的有向图,定义出度矩阵\(D^{out}(G)\)为:

\[D^{out}[i][i]=deg^{out}(i),D^{out}[i][j]=0,i \neq j \]

类似地定义入度矩阵\(D^{in}(G)\)
设\(\# e(i,j)\)为点\(i\)指向点\(j\)的边数,并定义邻接矩阵\(A\)为:

\[A[i][j]= \# e(i,j),i \neq j \]

定义出度\(Laplace\)矩阵\(L^{out}\)为:

\[L^{out}(G)=D^{out}(G)-A(G) \]

定义入度\(Laplace\)矩阵\(L^{in}\)为:

\[L^{in}(G)=D^{in}(G)-A(G) \]

记图\(G\)的以\(r\)为根的所有根向树形图个数为\(t^{root}(G,r)\),所谓根向树形图,是说这张图的基图是一棵树,所有的边全部指向父亲.

记图\(G\)的以\(r\)为根的所有叶向树形图个数为\(t^{leaf}(G,r)\),所谓叶向树形图,是说这张图的基图是一棵树,所有的边全部指向儿子.

定理叙述

矩阵树定理具有多种形式,其中用的较多的是定理1,3,4

定理1(矩阵树定理,无向图行列式形式)

对于任意的\(i\)都有

\[t(G)=det \ L(G) \begin{pmatrix}1,2,...,i-1,i+1,...,n\\1,2,...,i-1,i+1,...,n\end{pmatrix} \]

其中记号\(L(G) \begin{pmatrix}1,2,...,i-1,i+1,...,n\\1,2,...,i-1,i+1,...,n\end{pmatrix}\)代表矩阵\(L(G)\)的第\(1,2,...,i-1,i+1,...,n\)行与第\(1,2,...,i-1,i+1,...,n\)列组成的子矩阵.也就是说,无向图的\(Laplce\)矩阵具有这样的性质,他所有的\(n-1\)阶子矩阵都相等.

定理2(矩阵树定理,无向图特征值形式)

设\(\lambda[1], \lambda[2],...,\lambda[n-1]\)为\(L(G)\)的\(n-1\)个非零特征值,那么有

\[t(G)=\frac{1}{n}\lambda[1] \lambda[2]...\lambda[n-1] \]

定理3(矩阵树定理,有向图根向形式)

对于任意的\(k\),都有

\[t^{root}(G,k)=det \ L^{out}(G)\begin{pmatrix}1,2,...,k-1,k+1,...,n\\1,2,...,k-1,k+1,...,n\end{pmatrix} \]

如果要统计一张图所有的根向树形图,只要枚举所有的根\(k\)并对\(t^{root}(G,k)\)求和即可.

定理4(矩阵树定理,有向图叶形式)

对于任意的\(k\),都有

\[t^{leaf}(G,k)=det \ L^{in}(G)\begin{pmatrix}1,2,...,k-1,k+1,...,n\\1,2,...,k-1,k+1,...,n\end{pmatrix} \]

如果要统计一张图所有的叶向树形图,只需要枚举所有的根\(k\)并对\(t^{root}(G,k)\)求和即可

---

所以知道了这东西怎么求之后,我们只需要求出对应的\(n\)个\(n-1\)阶的行列式即可,所以我们直接用高斯消元求解多项式即可.而高斯消元我也没有系统的写过,所以顺便补一发高斯消元即可.

高斯消元

例题

思路

初中老师说过：我们可以加减消元或者代入消元,但是我们需要在程序里实现的时候，需要一种有规律可循的算法。所以我们选择加减消元，但用代入消元回带。
整体思路就是我们可以先在某一个式子里,用这个式子的\(x\)消去其他式子里的\(x\),然后在剩下的两个式子里在选择一个式子里的\(y\),用这个\(y\)去消掉最后剩下式子里的\(y\),那么最后一个方程组里就只有一个未知数\(z\)了.倘若\(z\)的系数是\(1\),那么我们把\(z\)回带到第二个式子中求出\(y\),再把\(y,z\)带入第一个式子就可以最终得出答案来了.
好像这个方法再数学逻辑上讲是特别麻烦的，但是却是一个通用性强的方法
但是我们最终还是要写代码来实现这个算法的,我们可以先用一个\(n \times (n+1)\)的矩阵来记录每一个式子的系数以及结果.
那么首先我们需要用第一列中（所有的\(x\)中）系数最大的来消其他两个式子。而这里为了方便起见，我们将这个选中的系数置为\(1\)，方便上例中地不断带 回原式的操作（这样在回带的时候就可以不考虑原本的系数了）。
由于最多也只能用\(double\)型存储，所以必然会有精度误差。但如果我们每次都选用最大系数的来消掉其他系数，就可以最大程度地来减小误差。
在置为\(1\)之后，我们需要来用这个式子去消其他的式子（别忘了每个式子的结果也要消）。那么在最后，我们只需要将这个矩阵的最右下角（也就是最后一个元的实际值）不断回带即可。

代码

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;
double map[111][111];
double ans[111];
double eps=1e-7;
int main(){
    int n;
    cin>>n;
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n+1;j++)
            scanf("%lf",&map[i][j]);
    for(int i=1;i<=n;i++){
        int r=i;
        for(int j=i+1;j<=n;j++)
            if(fabs(map[r][i])<fabs(map[j][i]))
                r=j;//find_the_biggest_number_of_the_first_column（at present) 
        if(fabs(map[r][i])<eps){
            printf("No Solution");
            return 0;
        }
        if(i!=r)swap(map[i],map[r]);//对换一行或一列,属于找最大当前系数的其中一步。（这样就可以只处理当前行的系数啦！） 
        double div=map[i][i];
        for(int j=i;j<=n+1;j++)
            map[i][j]/=div;
        for(int j=i+1;j<=n;j++){
            div=map[j][i];
            for(int k=i;k<=n+1;k++)
                map[j][k]-=map[i][k]*div;
        }
    }
    ans[n]=map[n][n+1];
    for(int i=n-1;i>=1;i--){
        ans[i]=map[i][n+1];
        for(int j=i+1;j<=n;j++)
            ans[i]-=(map[i][j]*ans[j]);
    }//回带操作
    for(int i=1;i<=n;i++)
        printf("%.2lf\n",ans[i]);
}

行列式求解思路

高斯消元没什么好讲的，就说一下要注意的细节吧。

1、高斯消元交换两行时答案要除以\((−1)\).(性质二)
2、假设处理到第\(i\)个方程，一般的高斯消元是用第\(i\)个方程减去第\(j(j>i\))个方程,所得的答案作为新的第\(j\)个方程，但求行列式的时候要用第\(j\)行减去第\(i\)行，所得答案作为新的第\(j\)行。因为方式一相当于某一行乘\((−1)\),另一行加到这一行上，这并不符合性质六。

然后答案就是主对角线的乘积。因为虽然高斯消元后的行列式为一个倒三角形，但可以按列来消，最后就只剩下主对角线了。或者这样说，行列式为一个倒三角，用最原始的算值方法，最后一行一定要选第n个数，倒数第二行要选第n-1个数，以此类推，就把主对角线都乘起来了，如果不这样选，答案就是0，对最终答案没用贡献。

#include <cstdio>
#include <cmath>
#include <algorithm>
#include <cstdlib>
#include <cstring>
#include <ctime>
#include <deque>
#include <queue>
#include <vector>
#include <map>
using namespace std;

const int mod=65521;
const int maxn=300;
typedef long long LL;

int n, m;
LL cnt[maxn][maxn];
LL ans, dv;

LL gcd(LL b, LL c)
{
	if (c==0) return b;
	else return gcd(c, b%c);
}
void init()
{
	for (int i=1; i<n; ++i)
		for (int j=1; j<n; ++j)
		{
			if (i==j) cnt[i][j]=(min(m, n-i)+min(m, i-1))%mod;
			else
				if (abs(i-j)<=m) cnt[i][j]=-1;
				else cnt[i][j]=0;
		}
}
void solve()
{
	for (int i=1, j=1; i<n && j<n; ++i, ++j)
	{
		int id=i;
		for (int k=i; k<n; ++k)
			if (cnt[k][j]!=0) { id=k; break; }
		if (cnt[id][j]==0) { --i; continue; }
		if (id!=i)
		{
			ans*=-1;
			for (int k=j; j<n; ++k) swap(cnt[id][k], cnt[i][k]);
		}
		for (int k=i+1; k<n; ++k)
		if (cnt[k][j]!=0)
		{
			int tmp1=gcd(abs(cnt[i][j]), abs(cnt[k][j]));
			int tmp2=cnt[i][j]/tmp1;
			tmp1=cnt[k][j]/tmp1;
			dv=dv*tmp2%mod;
			for (int p=j; p<n; ++p)
				cnt[k][p]=(cnt[k][p]*tmp2-cnt[i][p]*tmp1)%mod;
		}
	}
}
LL POW(LL b, LL c)
{
	LL s=1, cur=b;
	while (c)
	{
		if (c & 1) s=s*cur%mod;
		cur=cur*cur%mod;
		c>>=1;
	}
	return s;
}
int main()
{
	freopen("count.in", "r", stdin);
	freopen("count.out", "w", stdout);
	scanf("%d%d", &m, &n);
	init();
	ans=dv=1;
	solve();
	for (int i=1; i<n; ++i) ans=ans*cnt[i][i]%mod;
	printf("%I64d\n", (ans*POW(dv, mod-2)%mod+mod)%mod);
	return 0;
}

因为知道了自己是多么的菜,所以才要更加努力去追求那个永远也不可能实现的梦想