浅谈生成树计数问题,以SPOJ HIGH, BZOJ 4894, BZOJ 1016为例

最新推荐文章于 2020-04-19 18:32:38 发布

meopass

最新推荐文章于 2020-04-19 18:32:38 发布

阅读量301

点赞数

分类专栏：矩阵树定理数学文章标签：线性代数矩阵矩阵树定理

本文链接：https://blog.csdn.net/meopass/article/details/78856162

版权

数学同时被 2 个专栏收录

17 篇文章 0 订阅

订阅专栏

矩阵树定理

1 篇文章 0 订阅

订阅专栏

关于矩阵树定理详见周冬：生成树的计数。

给定一个包含n个节点，m条边的无向图，问图中的生成树的种类数有多少。

点我，喵=w=？
这就是一个最基本的生成树问题，由此我们可以引出生成树计数的矩阵树定理。

矩阵树定理：
一个无向图G的生成树的个数为其基尔霍夫矩阵的任意n-1阶主子式的行列式的绝对值。
G的度数矩阵: $D[i][j] = (i==j)?deg[i]:0;$
G的邻接矩阵: $A[i][j]$ 代表i到j有 $A[i][j]$ 种方案。
G的基尔霍夫矩阵 $K = D - A$ 。

那么我们只需要计算出行列式的值即可。
先用搞死小圆 ~~小圆：我招谁惹谁了~~ ，将基尔矩阵化成上三角矩阵，那么其行列式的值为 $\prod_{i = 1}^{n} K[i][i]$ 。

那么我们这题就做完了！

#define others
#ifdef poj
#include <iostream>
#include <cstring>
#include <cmath>
#include <cstdio>
#include <algorithm>
#include <vector>
#include <string>
#endif // poj
#ifdef others
#include <bits/stdc++.h>
#endif // others
//#define file
#define all(x) x.begin(), x.end()
using namespace std;
const double eps = 1e-8;
const double pi = acos(-1.0);
int dcmp(double x) {
    if(fabs(x)<=eps) return 0;
    return (x>0)?1:-1;
};
typedef long long LL;

void file() {
    freopen("a.in", "r", stdin);
//    freopen("1.txt", "w", stdout);
}

namespace Solver {
    const int maxn = 15;
    const int mod = 31011;
    int n, m;
    double C[maxn][maxn];

    double gauss(int n, double a[maxn][maxn]) {
        n--;
        for(int i = 0 ; i < n ; i++){
            int r = i;
            for(int j = i + 1 ; j < n; j++)
                if(fabs(a[j][i]) > fabs(a[r][i]))
                    r = j;
            if(dcmp(a[r][i]) == 0) return 0;
            if(r != i){
                for(int j = 0 ; j < n ; j++)
                    swap(a[i][j] , a[r][j]);
            }
            for(int j = n; j >= i; j--)
                for(int k = i + 1; k < n; k++)
                    a[k][j] -= a[k][i]/a[i][i] * a[i][j];
        }
        double ans = 1;
        for(int i = 0; i < n; i++)
            ans *= a[i][i];
        return fabs(ans);
    }

    int solve() {
        int t;
        scanf("%d", &t);
        while(t--) {
            scanf("%d%d", &n, &m);
            memset(C, 0, sizeof C);
            for(int i = 1; i <= m; i++) {
                int u, v;
                scanf("%d%d", &u, &v);
                u--, v--;
                C[u][u] ++, C[v][v] ++;
                C[u][v] --, C[v][u] --;
            }
            printf("%.0f\n", gauss(n, C));
        }
    }
};

int main() {
//    file();
    Solver::solve();
    return 0;
}

上述代码使用了浮点数除法，事实上如果题目要求的方案数对素数取模，我们可以用整数逆元来代替。
需要特别注意的是，我们浮点搞死小圆中找到最大绝对值的系数来消圆本来是为了优化精度，如果在整数的模意义下，这种行为就没有意义了，此时我们只需要找到第i行及以下任意一个系数非零的元来消即可。

我们可以从这个例子入手：BZOJ4894

给出一个有向图的邻接矩阵，求图中的生成树方案数mod1e9+7。

除去刚才我们说的需要注意的地方，这题还有一个要点是有向图。
对于有向图的基尔霍夫矩阵，其度数矩阵D只记录出度，并且我们计算的时候，只能计算去除根所在的行和列得到的 $n-1$ 阶主子式的行列式。

#define others
#ifdef poj
#include <iostream>
#include <cstring>
#include <cmath>
#include <cstdio>
#include <algorithm>
#include <vector>
#include <string>
#endif // poj
#ifdef others
#include <bits/stdc++.h>
#endif // others
//#define file
#define all(x) x.begin(), x.end()
using namespace std;
const double eps = 1e-8;
const double pi = acos(-1.0);
int dcmp(double x) {
    if(fabs(x)<=eps) return 0;
    return (x>0)?1:-1;
};
typedef long long LL;

void file() {
    freopen("a.in", "r", stdin);
//    freopen("1.txt", "w", stdout);
}

namespace Solver {
    const int maxn = 330;
    const int mod = 1e9+7;
    int n, m;
    LL Pow(LL a, LL b) {
        LL res = 1;
        while(b) {
            if(b & 1) res *= a, res %= mod;
            b >>= 1;
            a *= a; a %= mod;
        }
        return res;
    }
    LL C[maxn][maxn], D[maxn][maxn];

    LL gauss(int n, LL a[maxn][maxn]) {
        for(int i = 1 ; i < n ; i++){
            int r = i;
            for(int j = i ; j < n; j++)
                if(a[j][i]!=0) {
                    r = j;
                    break;
                }
            if(r != i){
                for(int j = 1 ; j < n ; j++)
                    swap(a[i][j] , a[r][j]);
            }
            LL inv = Pow(a[i][i], mod-2);
            for(int j = n - 1; j >= i; j--)
                for(int k = i + 1; k < n; k++) {
                    a[k][j] -= ((a[k][i]*inv)%mod * a[i][j])%mod;
                    if(a[k][j] <0) a[k][j] += mod;
                }
        }
        LL ans = 1;
        for(int i = 1; i < n; i++)
            ans *= a[i][i], ans %= mod;;
        return ans;
    }

    int solve() {
        scanf("%d", &n);
        for(int i = 0; i < n; i++)
            for(int j = 0; j < n; j++) {
                char c; scanf(" %c", &c);
                if(c == '1')
                    C[i][j] = mod-1, D[j][j]++;
            }
        for(int i = 0; i < n; i++)
            for(int j = 0; j < n; j++) {
                C[i][j] += D[i][j];
                C[i][j] %= mod;
            }
        printf("%lld\n", gauss(n, C));
    }
};

int main() {
//    file();
    Solver::solve();
    return 0;
}

BZOJ1016

给出n个点，m条带权边，问原图中最小生成树的方案数是多少。

我们需要知道一个性质，最小生成树每个权值能起的作用是一定的。
换句话说，即使不同的最小生成树方案，对于某个权值即使所用的边的数目是相同的，且这些边连接了相同的连通块。
做法1：
先跑一次最小生成树，获得每个边权在最小生成树中用了几次。
由于每个权值的作用是独立的，且相同权值的边数不超过10，我们先枚举每个权值对答案的贡献，在枚举当前权值的时候，把当前权值从最小生成树中抠掉，然后二进制枚举出当前权值的每条边用或不用，看能否组成最小生成树。乘法原理统计答案即可。
需要注意边权很大，需要离散化。

namespace solver {
    const int maxn = 1111;
    const int mod = 31011;
    namespace DSU {
        int fa[maxn], sz[maxn];
        stack<pair<int*, int> > stk;
        void init() {
            for(int i = 0; i < maxn; i++) fa[i] = i, sz[i] = 1;
        }
        int find(int x) {
            return x == fa[x]? fa[x] : find(fa[x]);
        };
        void merge(int x, int y, int on) {
            int u = find(x), v = find(y);
            if(sz[u] > sz[v]) swap(u, v);
            if(on) stk.push({&fa[u], fa[u]}), stk.push({&sz[v], sz[u]});
            fa[u] = v;
            sz[v] += sz[u];
        }
        void goback() {
            while(!stk.empty()) (*stk.top().first) = stk.top().second, stk.pop();
        }
    }
    int n, m;
    struct A {
        int u, v, w;
        bool operator < (const A & b) const {
            return w < b.w;
        }
    };
    int bit_count(int x) {
        return x == 0?0:x%2+bit_count(x/2);
    }
    vector<A> G[maxn];
    int cnt[maxn];
    int solve() {
        scanf("%d%d", &n, &m);
        vector<A> V, used;
        map<int, int> mid;
        int id = 0;
        int res = 1;
        for(int i = 1; i <= m; i++) {
            int u, v, w;
            scanf("%d%d%d", &u, &v, &w);
            if(mid.count(w) == 0)
                mid[w] = ++id;
            G[mid[w]].push_back({u, v, w});
            V.push_back({u, v, w});
        }
        sort(all(V));
        DSU::init();
        int tt = 0;
        for(int i = 0; i < V.size(); i++) {
            int f1 = DSU::find(V[i].u), f2 = DSU::find(V[i].v);
            if(f1 != f2) {
                used.push_back(V[i]);
                DSU::merge(f1, f2, 0);
                cnt[mid[V[i].w]]++;
                tt++;
            }
        }
        if(tt != n - 1) return 0 * puts("0");

        for(int i = 0; i < maxn; i++) {
            if(cnt[i] == 0) continue;
            DSU::init();
            int tmp = 0;
            for(int j = 0; j < used.size(); j++) {
                if(mid[used[j].w] == i) continue;
                A &buf = used[j];
                if(DSU::find(buf.u) != DSU::find(buf.v))
                    DSU::merge(buf.u, buf.v, 0), tmp++;
            }
            int len = G[i].size();
            int ptmp = tmp, ans = 0;
            for(int j = 0; j < (1 << len); j++) {
                if(bit_count(j) == cnt[i]) {
                    for(int k = 0; k < len; k++)
                        if(j & (1 << k)) {
                            A &buf = G[i][k];
                            if(DSU::find(buf.u) != DSU::find(buf.v))
                                DSU::merge(buf.u, buf.v, 1), tmp++;
                        }
                    if(tmp == n - 1)
                        ans++;
                    tmp = ptmp;
                    DSU::goback();
                }
            }
            res *= ans;
            res %= mod;
        }
        printf("%d\n", res);
    }
}

做法2：
依旧是枚举边权，统计每个边权的贡献，不过我们不再是二进制枚举，而是把当前权值的边单独拿出来，剩下的各个连通块用并查集缩点之后建新图，这样我们就可以用矩阵树定理了！
算完之后依旧是乘法原理统计答案。
需要注意：31011不是质数！并且这题方案数很少，浮点数的精度是够用的。

#define others
#ifdef poj
#include <iostream>
#include <cstring>
#include <cmath>
#include <cstdio>
#include <algorithm>
#include <vector>
#include <string>
#endif // poj
#ifdef others
#include <bits/stdc++.h>
#endif // others
//#define file
#define all(x) x.begin(), x.end()
using namespace std;
const double eps = 1e-8;
const double pi = acos(-1.0);
int dcmp(double x) {
    if(fabs(x)<=eps) return 0;
    return (x>0)?1:-1;
};
typedef long long LL;

void file() {
    freopen("a.in", "r", stdin);
//    freopen("1.txt", "w", stdout);
}

namespace Solver {
    const int maxn = 1111;
    const int mod = 31011;
    int n;
    double C[111][111];

    double gauss(int n, double a[111][111]) {
        n--;
        for(int i = 0 ; i < n ; i++){
            int r = i;
            for(int j = i + 1 ; j < n; j++)
                if(fabs(a[j][i]) > fabs(a[r][i]))
                    r = j;
            if(r != i){
                for(int j = 0 ; j <= n ; j++)
                    swap(a[i][j] , a[r][j]);
            }
            for(int j = n; j >= i; j--)
                for(int k = i + 1; k < n; k++)
                    a[k][j] -= a[k][i]/a[i][i] * a[i][j];
        }
        double ans = 1;
        for(int i = 0; i < n; i++)
            ans *= a[i][i];
        return fabs(ans);
    }

    namespace DSU {
        int fa[maxn], sz[maxn];
        stack<pair<int*, int> > stk;
        void init() {
            for(int i = 0; i < maxn; i++) fa[i] = i, sz[i] = 1;
        }
        int find(int x) {
            return x == fa[x]? fa[x] : find(fa[x]);
        };
        void merge(int x, int y, int on) {
            int u = find(x), v = find(y);
            if(sz[u] > sz[v]) swap(u, v);
            if(on) stk.push({&fa[u], fa[u]}), stk.push({&sz[v], sz[u]});
            fa[u] = v;
            sz[v] += sz[u];
        }
        void goback() {
            while(!stk.empty()) (*stk.top().first) = stk.top().second, stk.pop();
        }
    }
    int m;
    struct A {
        int u, v, w;
        bool operator < (const A & b) const {
            return w < b.w;
        }
    };
    int bit_count(int x) {
        return x == 0?0:x%2+bit_count(x/2);
    }
    vector<A> G[maxn];
    int cnt[maxn];
    vector<int> G_id;
    int get(int x) {
        return lower_bound(all(G_id), x) - G_id.begin();
    }
    int solve() {
        scanf("%d%d", &n, &m);
        vector<A> V, used;
        map<int, int> mid;
        int id = 0;
        int res = 1;
        for(int i = 1; i <= m; i++) {
            int u, v, w;
            scanf("%d%d%d", &u, &v, &w);
            if(mid.count(w) == 0)
                mid[w] = ++id;
            G[mid[w]].push_back({u, v, w});
            V.push_back({u, v, w});
        }
        sort(all(V));
        DSU::init();
        int tt = 0;
        for(int i = 0; i < V.size(); i++) {
            int f1 = DSU::find(V[i].u), f2 = DSU::find(V[i].v);
            if(f1 != f2) {
                used.push_back(V[i]);
                DSU::merge(f1, f2, 0);
                cnt[mid[V[i].w]]++;
                tt++;
            }
        }
        if(tt != n - 1) return 0 * puts("0");

        for(int i = 0; i < maxn; i++) {
            if(cnt[i] == 0) continue;
            DSU::init();
            int tmp = 0;
            for(int j = 0; j < used.size(); j++) {
                if(mid[used[j].w] == i) continue;
                A &buf = used[j];
                if(DSU::find(buf.u) != DSU::find(buf.v))
                    DSU::merge(buf.u, buf.v, 0), tmp++;
            }
            int len = G[i].size();
            G_id.clear();
            for(int p = 1; p <= n; p++) G_id.push_back(DSU::find(p));
            sort(all(G_id)); G_id.erase(unique(all(G_id)), G_id.end());
            memset(C, 0, sizeof C);
            for(int j = 0; j < len; j++) {
                A &buf = G[i][j];
                int u = get(DSU::find(buf.u)), v = get(DSU::find(buf.v));
                C[u][u] ++, C[v][v]++;
                C[u][v] --, C[v][u] --;
            }
            LL ans = 0;
            ans = gauss(G_id.size(), C)+0.0000001;
            res *= ans;
            res %= mod;
        }
        printf("%d\n", res);
    }
};
//g++ "C:\Users\meopass\Desktop\code\solution.cpp" -o solution && solution
int main() {
//    file();
    Solver::solve();
    return 0;
}

呼…暂时就到这里了。

meopass

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
浅谈生成树计数问题,以SPOJ HIGH, BZOJ 4894, BZOJ 1016为例

关于矩阵树定理详见周冬：生成树的计数。给定一个包含n个节点，m条边的无向图，问图中的生成树的种类数有多少。点我，喵=w=？这就是一个最基本的生成树问题，由此我们可以引出生成树计数的矩阵树定理。矩阵树定理：一个无向图G的生成树的个数为其基尔霍夫矩阵的任意n-1阶主子式的行列式的绝对值。 G的度数矩阵: D[i][j]=(i==j)?deg[i]:0;D[i][j] = (i==j)?deg[
复制链接

扫一扫