#include <iostream>
#include <vector>
using namespace std;

long long N, A[100009], B[100009], dp[100009], fact[5009], ans[5009][5009][2], mod = 1000000007;
vector<int> X[100009]; int nr[5009][5009];

void dfs(int pos) {
	dp[pos] = 1;
	for (int i = 0; i < X[pos].size(); i++) {
		if (dp[X[pos][i]] >= 1) continue;
		dfs(X[pos][i]);
		dp[pos] += dp[X[pos][i]];
	}
}

void init() {
	fact[0] = 1; for (int i = 1; i <= N; i++) fact[i] = (1LL * i * fact[i - 1]) % mod;
	for (int i = 0; i <= N; i++) {
		for (int j = 0; j <= N; j++) {
			if (i == 0 || j == 0) nr[i][j] = 1;
			else nr[i][j] = (nr[i - 1][j] + nr[i][j - 1]) % mod;
		}
	}
}

long long ncr(long long n, long long r) {
	if (r < 0 || n < r) return 0;
	return nr[r][n - r];
}

int main() {
	cin >> N; init();
	for (int i = 1; i <= N - 1; i++) {
		cin >> A[i] >> B[i];
		X[A[i]].push_back(B[i]);
		X[B[i]].push_back(A[i]);
	}
	dfs(1); vector<int> vec;
	for (int j = 1; j <= N; j++) {
		int sum = N - 1; bool flag = true;
		for (int k = 0; k < X[j].size(); k++) {
			if (dp[j] <= dp[X[j][k]]) continue;
			if (dp[X[j][k]] > (N / 2)) flag = false;
			sum -= dp[X[j][k]];
		}
		if (sum > (N / 2)) flag = false;
		if (flag == true) vec.push_back(j);
	}

	if (vec.size() == 2) {
		cout << (fact[N / 2] * fact[N / 2]) % mod << endl;
	}
	else {
		vector<int> D; int sum = N - 1;
		for (int i = 0; i < X[vec[0]].size(); i++) {
			if (dp[X[vec[0]][i]] < dp[vec[0]]) { D.push_back(dp[X[vec[0]][i]]); sum -= D[D.size() - 1]; }
		}
		D.push_back(sum);

		ans[0][0][0] = 1;
		for (int i = 0; i < D.size(); i++) {
			for (int j = 0; j <= N; j++) {
				if (ans[i][j][0] + ans[i][j][1] == 0) continue;
				for (int k = 0; k <= D[i]; k++) {
					long long K = ncr(D[i], k); K *= K; K %= mod; K *= fact[k]; K %= mod;
					int L = 0; if (k >= 1) L = 1;
					ans[i + 1][j + k][(k + 0) & 1] += ans[i][j][0] * K; ans[i + 1][j + k][(k + 0) & 1] %= mod;
					ans[i + 1][j + k][(k + 1) & 1] += ans[i][j][1] * K; ans[i + 1][j + k][(k + 1) & 1] %= mod;
				}
			}
		}

		long long ret = 0;
		for (int i = 0; i <= N; i++) {
			ret += ans[D.size()][i][0] * fact[N - i];
			ret -= ans[D.size()][i][1] * fact[N - i];
			ret += mod * mod; ret %= mod;
		}
		cout << ret << endl;
	}
	return 0;
}
### Prompt

Generate a CPP solution to the following problem:

There is a tree with N vertices. The vertices are numbered 1 through N. The i-th edge (1 \leq i \leq N - 1) connects Vertex x_i and y_i. For vertices v and w (1 \leq v, w \leq N), we will define the distance between v and w d(v, w) as "the number of edges contained in the path v-w".

A squirrel lives in each vertex of the tree. They are planning to move, as follows. First, they will freely choose a permutation of (1, 2, ..., N), p = (p_1, p_2, ..., p_N). Then, for each 1 \leq i \leq N, the squirrel that lived in Vertex i will move to Vertex p_i.

Since they like long travels, they have decided to maximize the total distance they traveled during the process. That is, they will choose p so that d(1, p_1) + d(2, p_2) + ... + d(N, p_N) will be maximized. How many such ways are there to choose p, modulo 10^9 + 7?

Constraints

* 2 \leq N \leq 5,000
* 1 \leq x_i, y_i \leq N
* The input graph is a tree.

Input

Input is given from Standard Input in the following format:


N
x_1 y_1
x_2 y_2
:
x_{N - 1} y_{N - 1}


Output

Print the number of the ways to choose p so that the condition is satisfied, modulo 10^9 + 7.

Examples

Input

3
1 2
2 3


Output

3


Input

4
1 2
1 3
1 4


Output

11


Input

6
1 2
1 3
1 4
2 5
2 6


Output

36


Input

7
1 2
6 3
4 5
1 7
1 5
2 3


Output

396

### Response

```cpp
#include <iostream>
#include <vector>
using namespace std;

long long N, A[100009], B[100009], dp[100009], fact[5009], ans[5009][5009][2], mod = 1000000007;
vector<int> X[100009]; int nr[5009][5009];

void dfs(int pos) {
	dp[pos] = 1;
	for (int i = 0; i < X[pos].size(); i++) {
		if (dp[X[pos][i]] >= 1) continue;
		dfs(X[pos][i]);
		dp[pos] += dp[X[pos][i]];
	}
}

void init() {
	fact[0] = 1; for (int i = 1; i <= N; i++) fact[i] = (1LL * i * fact[i - 1]) % mod;
	for (int i = 0; i <= N; i++) {
		for (int j = 0; j <= N; j++) {
			if (i == 0 || j == 0) nr[i][j] = 1;
			else nr[i][j] = (nr[i - 1][j] + nr[i][j - 1]) % mod;
		}
	}
}

long long ncr(long long n, long long r) {
	if (r < 0 || n < r) return 0;
	return nr[r][n - r];
}

int main() {
	cin >> N; init();
	for (int i = 1; i <= N - 1; i++) {
		cin >> A[i] >> B[i];
		X[A[i]].push_back(B[i]);
		X[B[i]].push_back(A[i]);
	}
	dfs(1); vector<int> vec;
	for (int j = 1; j <= N; j++) {
		int sum = N - 1; bool flag = true;
		for (int k = 0; k < X[j].size(); k++) {
			if (dp[j] <= dp[X[j][k]]) continue;
			if (dp[X[j][k]] > (N / 2)) flag = false;
			sum -= dp[X[j][k]];
		}
		if (sum > (N / 2)) flag = false;
		if (flag == true) vec.push_back(j);
	}

	if (vec.size() == 2) {
		cout << (fact[N / 2] * fact[N / 2]) % mod << endl;
	}
	else {
		vector<int> D; int sum = N - 1;
		for (int i = 0; i < X[vec[0]].size(); i++) {
			if (dp[X[vec[0]][i]] < dp[vec[0]]) { D.push_back(dp[X[vec[0]][i]]); sum -= D[D.size() - 1]; }
		}
		D.push_back(sum);

		ans[0][0][0] = 1;
		for (int i = 0; i < D.size(); i++) {
			for (int j = 0; j <= N; j++) {
				if (ans[i][j][0] + ans[i][j][1] == 0) continue;
				for (int k = 0; k <= D[i]; k++) {
					long long K = ncr(D[i], k); K *= K; K %= mod; K *= fact[k]; K %= mod;
					int L = 0; if (k >= 1) L = 1;
					ans[i + 1][j + k][(k + 0) & 1] += ans[i][j][0] * K; ans[i + 1][j + k][(k + 0) & 1] %= mod;
					ans[i + 1][j + k][(k + 1) & 1] += ans[i][j][1] * K; ans[i + 1][j + k][(k + 1) & 1] %= mod;
				}
			}
		}

		long long ret = 0;
		for (int i = 0; i <= N; i++) {
			ret += ans[D.size()][i][0] * fact[N - i];
			ret -= ans[D.size()][i][1] * fact[N - i];
			ret += mod * mod; ret %= mod;
		}
		cout << ret << endl;
	}
	return 0;
}
```