Luogu P4336 [SHOI2016]黑暗前的幻想乡（容斥，矩阵树定理）

Posted 2021-07-29 繁凡さん

整理的算法模板合集： ACM模板

点我看算法全家桶系列！！！

实际上是一个全新的精炼模板整合计划

---

Luogu P4336 [SHOI2016]黑暗前的幻想乡（容斥，矩阵树定理）

Problem

$n\le 17,P=10^9+7$

Solution

题目就是要求由 $n-1$ 个公司每个公司一条边建成的生成树的方案数。

显然求生成树的方案数，我们直接用矩阵树定理计算即可。

现在考虑满足 每个公司都只负责一条边的方案数 如何计算。

显然有：$n\le 17+计数问题=\text{容斥}$

我们可以先用矩阵树定理计算 $n-1$ 个公司包含的所有边集的生成树的个数，显然这里算出来的生成树，不一定 $n-1$ 个公司都参与了建设，我们用容斥原理减去不合法的即可。

显然就是枚举有多少个公司没有参与建设，然后利用容斥原理奇加偶减，答案减去容斥原理计算出的结果，变成奇减偶加，即：减去 $1$ 个公司没有参与，由剩下的 $n-2$ 个公司建成的生成树的个数，加上 $2$ 个公司没有参与，由剩下的 $n-3$ 个公司建成的生成树的个数，减去 $3$ 个公司没有参与，由剩下的 $n-4$ 个公司建成的生成树的个数 $\cdots$

我们可以通过二进制枚举来枚举具体选择了那几个公司的边，处理出此时的基尔霍夫矩阵，然后直接利用矩阵树定理高斯消元 $O(n^3)$ 计算代数余子式的行列式即可。

注意我们矩阵树定理计算的时候是可以处理重边的，所以如果一条边被多个公司覆盖，就把它当成重边即可，这样都当成重边算，最后减下来是没有重边的。

Time

$O(2^n\times n^3\log P)$

Code

// Problem: P4336 [SHOI2016]黑暗前的幻想乡
// Contest: Luogu
// URL: https://www.luogu.com.cn/problem/P4336
// Memory Limit: 250 MB
// Time Limit: 2000 ms
// 
// Powered by CP Editor (https://cpeditor.org)

#include <bits/stdc++.h>
#define int long long
using namespace std;

const int N = 100 + 7, M = 1e5 + 7, mod = 1e9 + 7;
 
int n, m[M], s, t, k, ans, a[N][N]; 
int siz[M];
int u[N][N], v[N][N];

int qpow(int a, int b)
{
	int res = 1;
	while(b) {
		if(b & 1) res = res * a % mod;
		a = a * a % mod;
		b >>= 1;
	}
	return res;
}

int guass2(int n)// 辗转相除 = TLE
{
	int det = 1;
	for (int i = 2; i <= n; ++ i) {
		for (int j = i + 1; 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] - t * a[j][k] + mod) % mod;
				swap(a[i], a[j]);
				det = -det;
			} 
		}
		det = (det * a[i][i]) % mod;
		if(det == 0) return 0;
	}
	return (det + mod) % mod;
}

int guass(int n)// 求逆元 = AC
{
	int det = 1;
	for (int i = 2; i <= n; ++ i) {
		for (int j = i; j <= n; ++ j) {
			if(a[i][j]) {
				swap(a[i], a[j]);
				if(i != j) 
					det = mod - det;
				break;
			}
		}

		int inv = qpow(a[i][i], mod - 2);

		for (int j = i + 1; j <= n; ++ j) {
			if(a[j][i]) {
				int tmp = a[j][i] * inv % mod;
				for (int k = i; k <= n; ++ k) {
					a[j][k] = (a[j][k] - a[i][k] * tmp % mod + mod) % mod;
				}
			}
		}
	}
	for (int i = 2; i <= n; ++ i)
		det = det * a[i][i] % mod;
	return det;
}

signed main()
{
	scanf("%lld", &n); 
	for (int i = 1; i <= n - 1; ++ i) {
		scanf("%lld", &m[i]);
		for (int j = 1; j <= m[i]; ++ j) {
			scanf("%lld%lld", &u[i][j], &v[i][j]);

			int x = u[i][j], y = v[i][j];
			a[x][x] ++ ;
			a[y][y] ++ ;
			a[x][y] = (a[x][y] + mod - 1) % mod;
			a[y][x] = (a[y][x] + mod - 1) % mod;
		}
	} 
	ans = guass(n);

	for (int i = 1; i <= (1 << (n - 1)) - 1; ++ i) {
		for (int j = 1; j <= n; ++ j) 
			for (int k = 1; k <= n; ++ k) 
				a[j][k] = 0;
			
		int cnt = 0;
		int tmp = i, j;

		for (j = 1; tmp; tmp >>= 1, ++ j) {
			if((tmp & 1) == 0) continue;
			cnt ++ ; 
			for (int k = 1; k <= m[j]; ++ k) {
				int x = u[j][k], y = v[j][k];
				a[x][x] ++ ;
				a[y][y] ++ ;
				a[x][y] = (a[x][y] + mod - 1) % mod;
				a[y][x] = (a[y][x] + mod - 1) % mod;
			}
		}
		if(cnt == n - 1) continue;
		ans = (ans - ((((n - 1) - cnt) & 1) ? 1 : -1) * guass(n) + mod) % mod;
	}
	printf("%lld\\n", ans);
	return 0;
}