关于矩阵树定理，其实没有太多要讲的东西，不过你需要先了解行列式。

行列式

定义

首先要了解行列式是一个函数，并不是一个数。它的自变量是一个方阵 AA，它的因变量是一个标量(一个数)，记作 det⁡(A)\det (A) 或 ∣A∣|A|。例如下图：

A=[abcd]det⁡(A)=∣abcd∣=ad−bc\begin{array}{clr} A= \begin{bmatrix} a & b\\ c & d \end{bmatrix} & \det(A)=\begin{vmatrix} a & b\\ c & d \end{vmatrix}=ad-bc \end{array}

行列式最早被用来解线性方程组，一个线性方程组有唯一解当且仅当它对应的行列式不为零。

对于一个 n×nn\times n 矩阵 AA，它的行列式的计算方法如下：

det⁡(A)=∑P(−1)μ(P)∏i=1nAi,Pi\det(A)=\sum _{P}(-1) ^{\mu (P)}\prod _{i=1} ^{n}A _{i,P _{i}}

其中 PP 是一个 1∼n1\sim n 的排列，μ(P)\mu(P) 为排列 PP 的逆序对数。

即每行挑一个数乘起来，再用得到的排列的逆序对数有关的东西做系数相加，这样计算复杂度为 O(n!×n)O(n!\times n)。

性质

首先了解行列式的几条基本性质。

  • 单位矩阵 II 的行列式是 11；
  • ∣a+eb+fcd∣=∣abcd∣+∣efcd∣\begin{vmatrix}a+e & b+f\\c & d\end{vmatrix}=\begin{vmatrix}a & b\\c & d\end{vmatrix}+\begin{vmatrix}e & f\\c & d\end{vmatrix}；
  • 有某两行一样或成比例的矩阵，行列式为 00；
  • 有某一行或某一列是 00 的矩阵，行列式为 00；
  • 行列互换，行列式不变。

接着了解初等行变换与行列式之间的几条关系。

  • 交换矩阵两行，行列式变为其相反数。

    交换两行后，只会影响到逆序对系数。易证一个排列中交换两个数对逆序对的影响一定是奇数个，因此行列式要变号。

  • 某一行乘以 kk，行列式也乘 kk。

    根据公式就可以看出来，那一行的所有数都乘以 kk，所以连乘时也会变成原来的 kk 倍。每个不同的排列中，连乘都会变成原来的 kk 倍，故行列式值变为原来 kk 倍。

  • 矩阵的一行加上另一行的倍数，行列式不变。

    有基本性质第二三四条合起来即可证明。

计算

首先是著名的消元求行列式。

知道初等行变换对行列式的影响后，我们由于消元只用了初等行变换，那么我们在高斯消元的时候只用在每次交换两行时把答案变号，最后消成上三角矩阵后，行列式再乘以对角线乘积，如果消元不出来就返回 00。时间复杂度为 O(n3)O(n ^{3})。

然后是展开求行列式。依然需要了解一些基本定义。

主子式

在 nn 阶行列式中，选取行号，再选取与行号相同的列号，则行数和列数都为 ii 个的行列式即为 nn 阶行列式的 ii 阶主子式。

余子式

在 nn 阶行列式中，划去元 ai,ja _{i,j} 所在的第 ii 行和第 jj 列，剩下的元不改变原顺序所构成的 n−1n-1 阶行列式称为元 ai,ja _{i,j} 的余子式，把元 ai,ja _{i,j} 的余子式也计作 Mi,jM _{i,j}。

代数余子式

将某余子式 Mi,jM _{i,j} 的值再乘以 (−1)i+j(-1) ^{i+j} 得到一个数 Ai,jA _{i,j}，它即是元 ai,ja _{i,j} 的代数余子式。

一个行列式按行展开，那么行列式 DD 的计算方式为：

D=∑j=1nai,jAi,jD=\sum _{j=1} ^{n}a _{i,j}A _{i,j}

一个行列式按列展开，那么行列式 DD 的计算方式为：

D=∑i=1nai,jAi,jD=\sum _{i=1} ^{n}a _{i,j}A _{i,j}

主要用这两种方法来求行列式。对于一些特殊的行列式还有着特殊的求法(行列式的解法)。

接下来才是我们的重头戏——矩阵树定理。

我们先来看最最经典的情况，无向无权图求生成树个数。

我们设 AA 为邻接矩阵(若 i,ji,j 之间存在边，那么 ai,j=aj,i=1a _{i,j}=a _{j,i}=1，其他为 00)，DD 为度数矩阵(Di,iD _{i,i} 为 ii 的度数，其他的为 00)，那么它的基尔霍夫(Kirchoff\text{Kirchoff}) 矩阵 K=D−AK=D-A。设 K′K' 为 KK 的任意一个 n−1n-1 阶主子式，那么 det⁡(K′)\det(K') 就是该图的生成树个数。

然后我们扩展一下，把它变成无向带权图。

求生成树的个数就把边权当成 11 来处理即可。但如果我们把 DD 的定义改一下，把 Di,iD _{i,i} 变成与 ii 点相连的边的边权之和，同理 Ai,jA _{i,j} 也变成 i,ji,j 之间的所有边的边权之和。这样得出来的 det⁡(K′)\det(K') 就是所有生成树边权乘积的总和。

继续扩展至有向图。

这个时候要分两种情况考虑，我们把生成树分成外向树和内向树，顾名思义，外向树就是树上的每条边都是从根向外的，内向树书上每条边都是由外向根的。因此外向树也被称作叶向树，内向树也被称作根向树。

这时候我们把 Ai,jA _{i,j} 只算 i→ji\to j 的有向边即可。对于度数矩阵则要分情况。如果要求外向树，那么 Di,iD _{i,i} 就是从外面到该点的边权总和，即入度边权总和；内向树则是从该点到外面的边权总和，即出度边权总和。同时既然是有向，那么根就是确定的。K′K' 必须要是 KK 去掉根那一行那一列的主子式。

例 1

模板题。

Luogu6178 Matrix-Tree 定理

/*
 * @Author: clorf 
 * @Date: 2020-09-14 08:50:18 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-14 09:16:22
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#define INF 1e9
using namespace std;
const int maxn=310;
const int mod=1e9+7;
const double Pi=acos(-1.0);
const double eps=1e-7;
//1.数组开小
//2.没用long long/爆了类型存储范围
//3.写之前式子没推清楚
//4.变量写错(想当然typo/没想清楚含义)
//5.写之前没想清楚复杂度
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(ch<'0'||ch>'9') {f|=(ch=='-');ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
int n,m,t;
long long a[maxn][maxn],ans=1,flag=1;
long long power(long long a,long long b,long long p)
{
    long long ans=1%p;
    for(;b;b>>=1)
    {
        if(b&1)
            ans=ans*a%p;
        a=a*a%p;
    }
    return ans%p;
}
long long inv(long long x)
{
    return power(x,mod-2,mod);
}
inline void gauss()
{
    int i=2,j=2;
    while(i<=n&&j<=n)
    {
        int maxx=i;
        for(int k=i+1;k<=n;k++)
            if(abs(a[k][j])>abs(a[maxx][j]))
                maxx=k;
        if(abs(a[maxx][j]))
        {
            if(maxx!=i)
                flag*=-1;
            swap(a[maxx],a[i]);
            for(int k=i+1;k<=n;k++)
            {
                long long num=(a[k][j]*inv(a[i][j]))%mod;
                for(int l=j;l<=n;l++)
                    a[k][l]=(a[k][l]+(mod-num)*a[i][l])%mod;
            }
            i++;
        }
        j++;
    }
}
int main()
{
    scanf("%d%d%d",&n,&m,&t);
    for(int i=1;i<=m;i++)
    {
        int x,y,z;
        scanf("%d%d%d",&x,&y,&z);
        if(!t)
        {   
            (a[x][x]+=z)%=mod;
            (a[y][y]+=z)%=mod;
            (a[x][y]-=z)%=mod;
            (a[y][x]-=z)%=mod;
        }
        else
        {
            (a[y][y]+=z)%=mod;
            (a[x][y]-=z)%=mod;
        }
    }
    gauss();
    for(int i=2;i<=n;i++)
        ans=(ans*a[i][i]%mod)%mod;
    ans=(ans*flag+mod)%mod;
    printf("%lld",ans);
    return 0;
}

例 2

另一道模板题。

Luogu4111 小 Z 的房间

解析：为了保证每条边只枚举到 11 次，加边时只枚举某个点上方和左方的点，然后就是跑生成树个数。注意模数是 10910 ^{9}，不是质数，不能直接求逆元。这里我们采用类似欧几里得算法，辗转相除的方式，用这行把另一行消后(消元就相当于欧几里得中的取模)就交换两行，直到某一行为 00。

/*
 * @Author: clorf 
 * @Date: 2020-09-14 21:58:02 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-14 22:13:31
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#define INF 1e9
using namespace std;
const int maxn=15;
const int mod=1e9;
const double Pi=acos(-1.0);
const double eps=1e-7;
//1.数组开小
//2.没用long long/爆了类型存储范围
//3.写之前式子没推清楚
//4.变量写错(想当然typo/没想清楚含义)
//5.写之前没想清楚复杂度
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(ch<'0'||ch>'9') {f|=(ch=='-');ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
int n,m,cnt,num[maxn][maxn];
long long ans=1,a[maxn*maxn][maxn*maxn];
char s[maxn][maxn];
inline void add(int u,int v)
{
    a[u][u]++;
    a[v][v]++;
    a[u][v]--;
    a[v][u]--;
}
inline void gauss(int m,int n)
{
    int i=2,j=2;
    while(i<=m&&j<=n)
    {
        int maxx=i;
        for(int k=i+1;k<=m;k++)
            if(abs(a[k][j])>abs(a[maxx][j]))
                maxx=k;
        if(abs(a[maxx][j]))
        {
            if(maxx!=i)
                ans=-ans;
            swap(a[maxx],a[i]);
            for(int k=i+1;k<=m;k++)
            {
                while(a[k][j])
                {
                    int num=a[i][j]/a[k][j];
                    for(int l=j;l<=n;l++)
                        a[i][l]=(a[i][l]+(mod-num)*a[k][l])%mod;
                    swap(a[k],a[i]);
                    ans=-ans;
                }//辗转相除消元
            }
            i++;
        }
        j++;
    }
}
int main()
{
    scanf("%d%d",&n,&m);
    for(int i=1;i<=n;i++)
        scanf("%s",s[i]+1);
    for(int i=1;i<=n;i++)
        for(int j=1;j<=m;j++)
            if(s[i][j]=='.')
                num[i][j]=++cnt;
    for(int i=1;i<=n;i++)  
        for(int j=1;j<=m;j++)
            if(num[i][j])
            {
                if(num[i-1][j])
                    add(num[i][j],num[i-1][j]);
                if(num[i][j-1])
                    add(num[i][j],num[i][j-1]);
            }
    gauss(cnt,cnt);
    for(int i=2;i<=cnt;i++)
        ans=(ans*a[i][i])%mod;
    ans=(ans+mod)%mod;
    printf("%lld",ans);
    return 0;
}

例 3

Luogu2144 轮状病毒

解析：我们看到 n≤100n\le 100，100100 阶的行列式肯定要用高精度，加上高斯消元的复杂度后约为 O(n4)∼O(n5)O(n ^{4})\sim O (n ^{5})，需要卡一些常数才能卡过。有没有其他优秀的方法呢(比如能跑过 n≤100000n\le 100000 的那种)？我们先把 KK 矩阵画出来：

∣n−1−1−1⋯−1−1−13−10⋯0−1−1−13−1⋯00−10−13⋯00⋮⋮⋮⋮⋱⋮⋮−1000⋯3−1−1−100⋯−13∣n\begin{vmatrix} n & -1 & -1 & -1 & \cdots & -1 & -1\\ -1 & 3 & -1 & 0 & \cdots & 0 & -1\\ -1 & -1 & 3 & -1 & \cdots & 0 & 0\\ -1 & 0 & -1 & 3 & \cdots & 0 & 0\\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\ -1 & 0 & 0 & 0 & \cdots & 3 & -1\\ -1 & -1 & 0 & 0 & \cdots & -1 & 3 \end{vmatrix} _{n}

删掉第 11 行第 11 列得到 K′K'：

∣3−10⋯0−1−13−1⋯000−13⋯00⋮⋮⋮⋱⋮⋮000⋯3−1−100⋯−13∣n−1\begin{vmatrix} 3 & -1 & 0 & \cdots & 0 & -1\\ -1 & 3 & -1 & \cdots & 0 & 0\\ 0 & -1 & 3 & \cdots & 0 & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & 0 & \cdots & 3 & -1\\ -1 & 0 & 0 & \cdots & -1 & 3 \end{vmatrix} _{n-1}

然后我们把 K′K' 按第一列展开：

det⁡(K′)=3det⁡∣3−1⋯00−13⋯00⋮⋮⋱⋮⋮00⋯3−100⋯−13∣n−2+det⁡∣−10⋯0−1−13⋯00⋮⋮⋱⋮⋮00⋯3−100⋯−13∣n−2+(−1)n−1det⁡∣−10⋯0−13−1⋯00⋮⋮⋱⋮⋮00⋯−1000⋯3−1∣n−2\begin{array}{clr} \det(K') & =3\det \begin{vmatrix} 3 & -1 & \cdots & 0 & 0\\ -1 & 3 & \cdots & 0 & 0\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & \cdots & 3 & -1\\ 0 & 0 & \cdots & -1 & 3 \end{vmatrix} _{n-2}+\det \begin{vmatrix} -1 & 0 & \cdots & 0 & -1\\ -1 & 3 & \cdots & 0 & 0\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & \cdots & 3 & -1\\ 0 & 0 & \cdots & -1 & 3 \end{vmatrix} _{n-2}\\ & + (-1) ^{n-1}\det \begin{vmatrix} -1 & 0 & \cdots & 0 & -1\\ 3 & -1 & \cdots & 0 & 0\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & \cdots & -1 & 0\\ 0 & 0 & \cdots & 3 & -1 \end{vmatrix} _{n-2} \end{array}

我们把这 33 个行列式记作 an−2,bn−2,cn−2a _{n-2},b _{n-2},c _{n-2}。接下来我们发现 aa 是个具有特殊性质的行列式——三对角行列式，即只有主对角线和主对角线左右的两条对角线，每条对角线上的值都相同，其他地方为 00。这种行列式可以通过递推来求解，我们现在开始推 aa 的递推式。

把 ana _{n} 按第 11 行展开：

det⁡an=3det⁡∣3−1⋯00−13⋯00⋮⋮⋱⋮⋮00⋯3−100⋯−13∣n−1+det⁡∣−1−1⋯0003⋯00⋮⋮⋱⋮⋮00⋯3−100⋯−13∣n−1=3det⁡an−1+det⁡∣−1−1⋯0003⋯00⋮⋮⋱⋮⋮00⋯3−100⋯−13∣n−1\begin{aligned} \det a _{n} & =3\det\begin{vmatrix} 3 & -1 & \cdots & 0 & 0\\ -1 & 3 & \cdots & 0 & 0\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & \cdots & 3 & -1\\ 0 & 0 & \cdots & -1 & 3 \end{vmatrix} _{n-1}+\det\begin{vmatrix} -1 & -1 & \cdots & 0 & 0\\ 0 & 3 & \cdots & 0 & 0\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & \cdots & 3 & -1\\ 0 & 0 & \cdots & -1 & 3 \end{vmatrix} _{n-1}\\ & =3\det a _{n-1}+\det\begin{vmatrix} -1 & -1 & \cdots & 0 & 0\\ 0 & 3 & \cdots & 0 & 0\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & \cdots & 3 & -1\\ 0 & 0 & \cdots & -1 & 3 \end{vmatrix} _{n-1}\\ \end{aligned}

然后又得到了一个 n−1n-1 的主子式，继续按第一列展开它：

=3det⁡an−1−det⁡∣3⋯00⋮⋱⋮⋮0⋯3−10⋯−13∣n−2=3det⁡an−1−det⁡an−2\begin{aligned} & = 3\det a _{n-1}-\det\begin{vmatrix} 3 & \cdots & 0 & 0\\ \vdots & \ddots & \vdots & \vdots\\ 0 & \cdots & 3 & -1\\ 0 & \cdots & -1 & 3\\ \end{vmatrix} _{n-2}\\ & =3\det a _{n-1}-\det a _{n-2} \end{aligned}

如此，我们就得到了 ana _{n} 的递推式：

det⁡an=3det⁡an−1−det⁡an−2\det a _{n}=3\det a_{n-1}-\det a _{n-2}

用同样的方法，我们可以得到 bn,cnb _{n},c _{n} 与 ana _{n} 的关系，不过要注意系数 (−1)k(-1) ^{k} 不要漏掉。

bn−2=−an−3−1(−1)n−1cn−2=−an−3−1\begin{aligned} b _{n-2} & = -a _{n-3}-1\\ (-1) ^{n-1}c _{n-2} & = -a _{n-3}-1 \end{aligned}

带回去得到：

det⁡(K′)=3an−2−2an−3−2\det(K')=3a _{n-2}-2a _{n-3}-2

结合 ana _{n} 的递推式，自然而然猜测 det⁡(K′)\det (K') 也是二阶线性递推，待定系数可得：

det⁡(Kn′)=3det⁡(Kn−1′)−det⁡(Kn−2′)+2\det (K' _{n})=3\det(K' _{n-1})-\det(K' _{n-2})+2

于是直接递推即可，注意要写高精度。

/*
 * @Author: clorf 
 * @Date: 2020-09-14 09:17:52 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-14 21:54:01
 */
#include<cstdio>
#include<cstring>
#include<iostream>
#include<stdexcept>
#include<cmath>
#define INF 1e9
using namespace std;
const int maxn=100010;
const int M=10000,P=4;
const int mod=1e9+7;
const double Pi=acos(-1.0);
const double eps=1e-7; 
//1.数组开小
//2.没用long long/爆了类型存储范围
//3.写之前式子没推清楚
//4.变量写错(想当然typo/没想清楚含义)
//5.写之前没想清楚复杂度
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(ch<'0'||ch>'9') {f|=(ch=='-');ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
int n;
struct bign
{
    int x[20001],len;
    bign()
    {
        len=1;
        memset(x,0,sizeof(x));
    }
    bign(int x)
    {
        *this=x;
    }
    bign operator =(int n)
    {
        char s[20001];
        sprintf(s,"%d",n);
        *this=s;
    }
    bign operator =(char *s)
    {
        bign a;
        a.len=strlen(s);
        for(int i=0;i<a.len;i++)
            a.x[a.len-i]=s[i]-'0';
        *this=a;
    }
    friend bool operator <(bign a,bign b)
    {
        if(a.len<b.len) 
            return true;
        if(a.len==b.len)
        {
            int i=a.len;
            while(a.x[i]==b.x[i]&&i>=1)
                i--;
            if(a.x[i]<b.x[i])
                return true;
        }
        return false;
    }
    friend bool operator >(bign a,bign b)
    {
        return b<a;
    }
    friend bool operator <=(bign a,bign b)
    {
        return !(a>b);
    }
    friend bool operator >=(bign a,bign b)
    {
        return !(a<b);
    }
    friend bool operator ==(bign a,bign b)
    {
        return (a>=b)&&(b>=a);
    }
    friend bool operator !=(bign a,bign b)
    {
        return (a<b)||(b<a);
    }
    friend bign operator +(bign a,bign b)
    {
        int i=1,jw=0;
        bign c,aa=a,bb=b;
        while((i<=aa.len+1)||(i<=bb.len+1))
        {
            c.x[i]=aa.x[i]+bb.x[i]+jw;
            jw=c.x[i]/10;
            c.x[i]%=10;
            i++;
        }
        c.len=i;
        while(c.x[c.len]==0)
            c.len--;
        return c;
    }
    friend bign operator -(bign a,bign b)
    {
        if(a<b)
            swap(a,b);
        int i=1;
        bign c,aa=a,bb=b;
        while(i<=aa.len+1||i<=bb.len+1)
        {
            if(aa.x[i]<bb.x[i])
            {
                aa.x[i+1]--;
                aa.x[i]+=10;
            }
            c.x[i]=aa.x[i]-bb.x[i];
            i++;
        }
        c.len=i;
        while(c.x[c.len]==0)
            c.len--;
        return c;
    }
    friend bign operator *(bign a,int b)
    {
        bign c;
        int g=0,l=a.len;
        for(int i=1;i<=l;i++)
        {
            c.x[i]=a.x[i]*b+g;
            g=c.x[i]/10;
            c.x[i]%=10;
        }
        while(g>0)
        {
            l++;
            c.x[l]=g%10;
            g/=10;
        }
        while(l>1&&c.x[l]==0)
            l--;
        c.len=l;
    }
}; 
void cin_bign(bign &a)
{
    char s[20001];
    scanf("%s",s);
    a=s;
}
void cout_bign(bign a)
{
    bool ok=false;
    for(int i=a.len;i>0;i--)
    {
        printf("%d",a.x[i]);
        ok=1;
    }
    if(!ok) printf("0");
    printf("\n");
}
int main()
{
    scanf("%d",&n);
    bign x,y,z;
    if(n==1)
    {
        printf("1");
        return 0;
    }
    if(n==2)
    {
        printf("5");
        return 0;
    }
    x=5;
    y=1;
    for(int i=3;i<=n;i++)
    {
        z=x+x+x-y+2;
        y=x;
        x=z;
    }
    cout_bign(z);
    return 0;
}