【GDOI2020模拟02.05】生成树 (矩阵树扩展)

Description：

给定一张 N 个点，M 条边的无向图，边有红、绿、蓝三种颜色，分别用 1,2,3 表示。
求这张图有多少生成树，满足绿色边数量不超过 x，蓝色边数量不超过 y，答案对10^9 + 7 取模。

n<=40

题解：

考虑正常的用矩阵树求生成树个数的方法。

基尔霍夫矩阵里的每个位置只是一个数，事实上可以把它扩展成一个多项式。

在这题中，每个位置就是一个(系数*x^{...}*y^{...})。

然后求出这个矩阵的余子式的行列式就好了。

当然这个矩阵是不好高斯消元去求行列式的，但是我们可以拉格朗日插值。

二维拉格朗日插值：

(sum_{i,j} a_{x_i,y_j}prod_{k ot= i{x - xkover xi-xk}}prod_{k ot= j{y - ykover yj-yk}})

和一维没有什么区别，

我们需要暴力展开这个式子，直接搞的复杂度是(O(n^4))，加上高斯消元的(O(n^5))能够通过。

这题比较特殊，可以使(y=x^n)去做一维插值，我比赛时写的一维插值：
(sum_{i,j}a_{xi}prod_{k ot= i{x - xkover xi-xk}})

这格东西有(O(n^2))项，如果对后面的暴力展开，总共是(O(n^6))，需要卡常才能通过。

可以用除法优化，注意除一次也就(O(n^2))，总共(O(n^4))。

Code：

#include<bits/stdc++.h>
#define fo(i, x, y) for(int i = x, _b = y; i <= _b; i ++)
#define ff(i, x, y) for(int i = x, _b = y; i <  _b; i ++)
#define fd(i, x, y) for(int i = x, _b = y; i >= _b; i --)
#define ll long long
#define pp printf
#define hh pp("
")
using namespace std;

const int mo = 1e9 + 7;

ll ksm(ll x, ll y) {
	ll s = 1;
	for(; y; y /= 2, x = x * x % mo)
		if(y & 1) s = s * x % mo;
	return s;
}

const int N = 41;

int n;

ll a[N][N];

ll solve() {
	fo(i, 1, n) a[i][n] = a[n][i] = 0;
	n --;
	ll ans = 1;
	fo(i, 1, n) {
		int u = -1;
		fo(j, i, n) if(a[j][i])
			{ u = j; break;}
		if(u == -1) continue;
		if(u > i) {
			ans = -ans;
			fo(j, 1, n) swap(a[i][j], a[u][j]);
		}
		ll ni = ksm(a[i][i], mo - 2);
		fo(j, i + 1, n) if(a[j][i]) {
			ll v = -ni * a[j][i] % mo;
			fo(k, 1, n) a[j][k] = (a[j][k] + a[i][k] * v) % mo;
		}
	}
	fo(i, 1, n) ans = ans * a[i][i] % mo;
	ans = (ans % mo + mo) % mo;
	n ++;
	return ans;
}

int m, x, y, z;
int c1, c2;

ll p[3][N], q[3][N][N];

void Init() {
	scanf("%d %d %d %d", &n, &m, &c1, &c2);
	fo(i, 1, m) {
		scanf("%d %d %d", &x, &y, &z);
		if(x == y) continue;
		z --;
		p[z][x] ++; p[z][y] ++;
		q[z][x][y] ++; q[z][y][x] ++;
	}
}


ll cz(ll x) {
	memset(a, 0, sizeof a);
	ll x1 = x, x2 = ksm(x, n);
	fo(i, 1, n) {
		a[i][i] += p[0][i];
		a[i][i] += p[1][i] * x1 % mo;
		a[i][i] += p[2][i] * x2 % mo;
	}
	fo(i, 1, n) fo(j, 1, n) {
		a[i][j] -= q[0][i][j];
		a[i][j] -= q[1][i][j] * x1 % mo;
		a[i][j] -= q[2][i][j] * x2 % mo;
		a[i][j] %= mo;
	}
	return solve();
}

int t;

ll f[N * N], g[N * N], h[N * N];

ll ni[N * N];

void work() {
	t = n * (n - 1);
	fo(i, 0, t) ni[i] = ksm(i, mo - 2);
	
	g[0] = 1;
	fo(i, 0, t) {
		fd(j, i, 0) {
			g[j + 1] = (g[j + 1] + g[j]) % mo;
			g[j] = g[j] * -i % mo;
		}
	}
	
	fo(i, 0, t) {
		ll xs = cz(i);
		fo(j, 0, t) if(i != j)
			xs = xs * (i < j ? -ni[j - i] : ni[i - j]) % mo;
		
		fo(j, 0, t + 1) h[j] = g[j];
		fd(j, t + 1, 1) h[j - 1] = (h[j - 1] + h[j] * i) % mo;
		fo(j, 0, t) f[j] = (f[j] + xs * h[j + 1]) % mo;
	}
	
	ll ans = 0;
	fo(k, 0, t) {
		if(k % n <= c1 && k / n <= c2) {
			ans = (ans + f[k]) % mo;
		}
	}
	
	ans = (ans % mo + mo) % mo;
	pp("%lld
", ans);
}

int main() {
	freopen("tree.in", "r", stdin);
	freopen("tree.out", "w", stdout);
	Init();
	work();
}