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();
}