【算法】Matrix - Tree 矩阵树定理 & 题目总结

　　最近集中学习了一下矩阵树定理，自己其实还是没有太明白原理（证明）类的东西，但想在这里总结一下应用中的一些细节，矩阵树定理的一些引申等等。

　　首先，矩阵树定理用于求解一个图上的生成树个数。实现方式是：\(A\)为邻接矩阵，\(D\)为度数矩阵，则基尔霍夫（Kirchhoff）矩阵即为：\(K = D - A\)。具体实现中，记 \(a\) 为Kirchhoff矩阵，则若存在 \(E(u, v)\) ，则\(a[u][u] ++, a[v][v] ++, a[u][v] --, a[v][u] --\) 。即\(a[i][i]\) 为 \(i\) 点的度数，\(a[i][j]\) 为 \(i, j\)之间边的条数的相反数。

　　这样构成的矩阵的行列式的值，就为生成树的个数。而求解行列式的快速方法为使用高斯消元进行消元消处上三角矩阵，则有对角线上的值的乘积 = 行列式的值。一般而言求解生成树个数的题目数量会非常庞大，需要取模处理。取模处理中，不能出现小数，于是使用辗转相除法：（其中因为消的是行列式，所以与消方程有所不同。交换两行行列式的值变号，且消元只能将一行的数 * k 之后加到别的行上。）

int Gauss()

{

    int ans = ;

    for(int i = ; i < tot; i ++)

    {

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

            while(f[j][i])

            {

                int t = f[i][i] / f[j][i];

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

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

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

                ans = - ans;

            }

        ans = (ans * f[i][i]) % mod;

    }

    return (ans + mod) % mod;

}

　　变元矩阵树定理：求所有生成树的总边积的和。和矩阵树的求法相同，不过行列式中 \(a[i][i]\) 记录的是总的边权和，\(a[i][j]\) 记录 \(i, j\) 之间边权的相反数。

　　以下为几道题目：

　　　　1.HEOI2015 小Z的房间　　　　2.SHOI2016 黑暗前的幻想乡

　　　　3.SDOI2014 重建　　　　　　 4.JSOI2008 最小生成树计数

　　1.HEOI2015 小Z的房间（妥妥的模板题一个）

#include <bits/stdc++.h>

using namespace std;

#define maxn 90

#define int long long

#define mod 1000000000

int n, m, f[maxn][maxn];

int tot, Map[maxn][maxn];

int read()

{

    int x = , k = ;

    char c;

    c = getchar();

    while(c < '' || c > '') { if(c == '-') k = -; c = getchar(); }

    while(c >= '' && c <= '') x = x *  + c - '', c = getchar();

    return x * k;

}

void add(int x, int y)

{

    if(x > y) return;

    f[x][x] ++, f[y][y] ++;

    f[x][y] --, f[y][x] --;

}

int Gauss()

{

    int ans = ;

    for(int i = ; i < tot; i ++)

    {

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

            while(f[j][i])

            {

                int t = f[i][i] / f[j][i];

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

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

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

                ans = - ans;

            }

        ans = (ans * f[i][i]) % mod;

    }

    return (ans + mod) % mod;

}

signed main()

{

    n = read(), m = read();

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

    {

        char c;

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

        {

            cin >> c;

            if(c == '.') Map[i][j] = ++ tot;

        }

    }

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

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

        {

            int tem, u;

            if(!(u = Map[i][j])) continue;

            if(tem = Map[i - ][j]) add(u, tem);

            if(tem = Map[i + ][j]) add(u, tem);

            if(tem = Map[i][j - ]) add(u, tem);

            if(tem = Map[i][j + ]) add(u, tem);

        }

    printf("%lld\n", Gauss());

    return ;

}

　　2.SHOI2016黑暗前的幻想乡

　　容斥+矩阵树定理。与模板的不同之处在于每一家公司都要参与修建，则合法方案数 = 总的方案数 - 有一个公司未修建的方案数 + 有两个公司未修建的方案数……暴力重构矩阵求解即可。

#include <bits/stdc++.h>

using namespace std;

#define ll long long

const int mod = ;

int n;

ll g[][];

vector<pair<int , int > > q[];

int read()

{

    int x = , k = ;

    char c;

    c = getchar();

    while(c < '' || c > '') { if(c == '-') k = -; c = getchar(); }

    while(c >= '' && c <= '') x = x *  + c - '', c = getchar();

    return x * k;

}

int Gauss()

{

    ll ans = ;

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

    {

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

            while(g[j][i])

            {

                ll t = g[i][i] / g[j][i];

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

                    g[i][k] = (g[i][k] - g[j][k] * t) % mod;

                swap(g[i], g[j]);

                ans = -ans;

            }

        ans = (ans * g[i][i]) % mod;

        if(!ans) return ;

    }

    return (ans + mod) % mod;

}

int main()

{

    n = read();

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

    {

        int m = read();

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

        {

            int x = read(), y = read();

            q[i].push_back(make_pair(x, y));

        }

    }

    int ans = , CNST =  << (n - );

    for(int i = ; i < CNST; i ++)

    {

        int cnt = ; memset(g, , sizeof(g));

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

            if(i & ( << (j - )))

            {

                for(int k = ; k < q[j].size(); k ++)

                {

                    int x = q[j][k].first, y = q[j][k].second;

                    g[x][x] ++, g[y][y] ++;

                    g[x][y] --, g[y][x] --;

                }

                cnt ++;

            }

        if((n - cnt) & ) ans = (ans + Gauss()) % mod;

        else ans = (ans - Gauss() + mod) % mod;

    }

    printf("%d\n", ans);

    return ;

}

　　3.SDOI2014重建

　　化式子 + 变元矩阵树定理。将概率的式子写出来变形即可得到矩阵树定理求 \(\prod \frac{p(u, v)}{1 - p(u, v)}\)

#include <bits/stdc++.h>

using namespace std;

#define maxn 100

#define db double

#define eps 0.000001

int n;

db ans = 1.0, a[maxn][maxn];

db Gauss(int n)

{

    db ans = 1.0;

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

    {

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

        {

            int t = i;

            if(fabs(a[j][i]) > fabs(a[t][i])) t = j;

            if(t != i) swap(a[t], a[i]), ans = -ans;

        }

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

        {

            db t = a[j][i] / a[i][i];

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

                a[j][k] -= t * a[i][k];

        }

        ans *= a[i][i];

    }

    return fabs(ans);

}

int main()

{

    scanf("%d", &n);

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

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

            scanf("%lf", &a[i][j]);

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

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

        {

            db t = fabs(1.0 - a[i][j]) < eps ? eps : (1.0 - a[i][j]);

            if(i < j) ans *= t;

            a[i][j] = a[i][j] / t;

        }

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

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

            if(i != j) { a[i][i] += a[i][j], a[i][j] = -a[i][j]; }

    printf("%.10lf\n", Gauss(n - ) * ans);

    return ;

}

　　4.JSOI2008最小生成树计数

　　这题虽然最早年，然而也最强啊……个人认为这位博主解释得很好了 Z-Y-Y-S的博客

　　两个性质 mark 一下：

#include<iostream>

#include<cstdio>

#include<cstring>

#include<algorithm>

#include<cmath>

using namespace std;

#define maxn 200

#define mod 31011

int n, m, ans = , tmp[maxn];

int sum, fa[maxn], set[maxn];

int a[maxn][maxn];

struct edge

{

    int u, v, w;

}E[maxn * ], e[maxn * ];

int read()

{

    int x = , k = ;

    char c;

    c = getchar();

    while(c < '' || c > '') { if(c == '-') k = -; c = getchar(); }

    while(c >= '' && c <= '') x = x *  + c - '', c = getchar();

    return x * k;

}

bool cmp(edge a, edge b) { return a.w < b.w; }

int find(int x) { return set[x] == x ? x : set[x] = find(set[x]); }

int find2(int x) { return fa[x] == x ? x : fa[x] = find2(fa[x]); }

int Gauss(int n)

{

    int ans = ;

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

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

            a[i][j] = (a[i][j] + mod) % mod;

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

    {

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

            while(a[j][i])

            {

                int t = a[i][i] / a[j][i];

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

                    a[i][k] = (a[i][k] - 1ll * t * a[j][k] % mod + mod) % mod;

                swap(a[i], a[j]); ans = - ans;

            }

        ans = 1ll * ans * a[i][i] % mod;

    }

    return (ans + mod) % mod;

}

void Cal(int S, int T)

{

    int cnt = ;

    for(int i = S; i <= T; i ++)

    {

        e[i] = E[i];

        int p = find(e[i].u), q = find(e[i].v);

        e[i].u = p, e[i].v = q;

        if(p == q) continue;

        tmp[++ cnt] = p, tmp[++ cnt] = q;

    }

    sort(tmp + , tmp +  + cnt);

    cnt = unique(tmp + , tmp + cnt + ) - tmp - ;

    memset(a, , sizeof(a));

    for(int i = ; i <= cnt; i ++) fa[i] = i;

    for(int i = S; i <= T; i ++)

    {

        if(e[i].u == e[i].v) continue;

        int p = find(e[i].u), q = find(e[i].v);

        if(p != q) -- sum, set[p] = q;

        int u = lower_bound(tmp + , tmp + cnt + , e[i].u) - tmp;

        int v = lower_bound(tmp + , tmp + cnt + , e[i].v) - tmp;

        a[u][u] ++, a[v][v] ++;

        a[u][v] --, a[v][u] --;

        p = find2(u), q = find2(v);

        if(p != q) fa[p] = q;

    }

    for(int i = ; i <= cnt; i ++)

        if(find2(i) != find2(i - ))

        {

            int p = find2(i), q = find2(i - );

            a[p][p] ++, a[q][q] ++;

            a[p][q] --, a[q][p] --;

            fa[p] = q;

        }

    ans = 1ll * ans * Gauss(cnt - ) % mod;

}

int main()

{

    n = read(), m = read();

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

        E[i].u = read(), E[i].v = read(), E[i].w = read();

    sort(E + , E +  + m, cmp);

    for(int i = ; i <= n; i ++) set[i] = i;

    sum = n;

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

    {

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

            if(E[j].w != E[i].w) break;

        if(j - i > ) Cal(i, j - );

        else

        {

            int p = find(E[i].u), q = find(E[i].v);

            if(p != q) set[p] = q;

            sum --;

        }

    }

    if(sum > ) printf("");

    else printf("%d\n", ans);

    return ;

}

巴特西

【算法】Matrix - Tree 矩阵树定理 & 题目总结

最新文章

热门文章