#include <iostream>
#include <vector>
#define llint long long
#define mod 1000000007

using namespace std;

llint n;
vector<int> G[5005];
int size[5005];
bool used[5005];
llint dp[5005][5005];

const int FACT_MAX = 5005;
llint fact[FACT_MAX], fact_inv[FACT_MAX];

llint modpow(llint a, llint n)
{
	if(n == 0) return 1;
	if(n % 2){
		return ((a%mod) * (modpow(a, n-1)%mod)) % mod;
	}
	else{
		return modpow((a*a)%mod, n/2) % mod;
	}
}

void make_fact()
{
	llint val = 1;
	fact[0] = 1;
	for(int i = 1; i < FACT_MAX; i++){
		val *= i;
		val %= mod;
		fact[i] = val;
	}
	fact_inv[FACT_MAX-1] = modpow(fact[FACT_MAX-1], mod-2);
	for(int i = FACT_MAX-2; i >= 0; i--){
		fact_inv[i] = fact_inv[i+1] * (i+1) % mod;
	}
}

llint comb(llint n, llint k)
{
	llint ret = 1;
	ret *= fact[n];
	ret *= fact_inv[k], ret %= mod;
	ret *= fact_inv[n-k], ret %= mod;
	return ret;
}

int sizedfs(int v, int pre)
{
	int ret = 1;
	for(int i = 0; i < G[v].size(); i++){
		if(G[v][i] == pre) continue;
		if(used[G[v][i]]) continue;
		ret += sizedfs(G[v][i], v);
	}
	return size[v] = ret;
}

int centdfs(int v, int pre, int sz)
{
	for(int i = 0; i < G[v].size(); i++){
		if(G[v][i] == pre) continue;
		if(used[G[v][i]]) continue;
		if(size[G[v][i]] > sz/2) return centdfs(G[v][i], v, sz);
	}
	return v;
}

int main(void)
{
	cin >> n;
	int u, v;
	for(int i = 0; i < n-1; i++){
		cin >> u >> v;
		G[u].push_back(v);
		G[v].push_back(u);
	}
	make_fact();
	
	sizedfs(1, -1);
	int cent = centdfs(1, -1, n);
	if(n%2 == 0){
		bool flag = false;
		for(int i = 0; i < G[cent].size(); i++) if(size[G[cent][i]] == n/2) flag = true;
		if(flag){
			cout << fact[n/2] * fact[n/2] % mod << endl;
			return 0;
		}
	}
	
	sizedfs(cent, -1);
	vector<llint> vec;
	for(int i = 0; i < G[cent].size(); i++) vec.push_back(size[G[cent][i]]);
	
	int m = vec.size();
	dp[0][0] = 1;
	for(int i = 0; i < m; i++){
		for(int j = 0; j < n; j++){
			for(int k = 0; k <= vec[i]; k++){
				if(j+k<n){
					dp[i+1][j+k] += dp[i][j] * comb(vec[i], k) % mod * comb(vec[i], k) % mod * fact[k] % mod;
					dp[i+1][j+k] %= mod;
				}
			}
		}
	}
	llint ans = fact[n-1], tmp = 0;
	for(int i = 1; i < n; i++){
		if(i % 2) tmp += fact[n-1-i] * dp[m][i] % mod;
		else tmp += mod - fact[n-1-i] * dp[m][i] % mod;
		tmp %= mod;
	}
	ans += mod - tmp, ans %= mod;
	
	vec.push_back(1);
	m = vec.size();
	for(int i = 0; i <= m; i++){
		for(int j = 0; j <= n; j++){
			dp[i][j] = 0;
		}
	}
	dp[0][0] = 1;
	for(int i = 0; i < m; i++){
		for(int j = 0; j <= n; j++){
			for(int k = 0; k <= vec[i]; k++){
				if(j+k<=n){
					dp[i+1][j+k] += dp[i][j] * comb(vec[i], k) % mod * comb(vec[i], k) % mod * fact[k] % mod;
					dp[i+1][j+k] %= mod;
				}
			}
		}
	}
	llint ans2 = fact[n]; tmp = 0;
	for(int i = 1; i <= n; i++){
		if(i % 2) tmp += fact[n-i] * dp[m][i] % mod;
		else tmp += mod - fact[n-i] * dp[m][i] % mod;
		tmp %= mod;
	}
	ans2 += mod - tmp, ans %= mod;
	
	cout << (ans+ans2)%mod << endl;
	return 0;
}
### Prompt

Please formulate 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>
#define llint long long
#define mod 1000000007

using namespace std;

llint n;
vector<int> G[5005];
int size[5005];
bool used[5005];
llint dp[5005][5005];

const int FACT_MAX = 5005;
llint fact[FACT_MAX], fact_inv[FACT_MAX];

llint modpow(llint a, llint n)
{
	if(n == 0) return 1;
	if(n % 2){
		return ((a%mod) * (modpow(a, n-1)%mod)) % mod;
	}
	else{
		return modpow((a*a)%mod, n/2) % mod;
	}
}

void make_fact()
{
	llint val = 1;
	fact[0] = 1;
	for(int i = 1; i < FACT_MAX; i++){
		val *= i;
		val %= mod;
		fact[i] = val;
	}
	fact_inv[FACT_MAX-1] = modpow(fact[FACT_MAX-1], mod-2);
	for(int i = FACT_MAX-2; i >= 0; i--){
		fact_inv[i] = fact_inv[i+1] * (i+1) % mod;
	}
}

llint comb(llint n, llint k)
{
	llint ret = 1;
	ret *= fact[n];
	ret *= fact_inv[k], ret %= mod;
	ret *= fact_inv[n-k], ret %= mod;
	return ret;
}

int sizedfs(int v, int pre)
{
	int ret = 1;
	for(int i = 0; i < G[v].size(); i++){
		if(G[v][i] == pre) continue;
		if(used[G[v][i]]) continue;
		ret += sizedfs(G[v][i], v);
	}
	return size[v] = ret;
}

int centdfs(int v, int pre, int sz)
{
	for(int i = 0; i < G[v].size(); i++){
		if(G[v][i] == pre) continue;
		if(used[G[v][i]]) continue;
		if(size[G[v][i]] > sz/2) return centdfs(G[v][i], v, sz);
	}
	return v;
}

int main(void)
{
	cin >> n;
	int u, v;
	for(int i = 0; i < n-1; i++){
		cin >> u >> v;
		G[u].push_back(v);
		G[v].push_back(u);
	}
	make_fact();
	
	sizedfs(1, -1);
	int cent = centdfs(1, -1, n);
	if(n%2 == 0){
		bool flag = false;
		for(int i = 0; i < G[cent].size(); i++) if(size[G[cent][i]] == n/2) flag = true;
		if(flag){
			cout << fact[n/2] * fact[n/2] % mod << endl;
			return 0;
		}
	}
	
	sizedfs(cent, -1);
	vector<llint> vec;
	for(int i = 0; i < G[cent].size(); i++) vec.push_back(size[G[cent][i]]);
	
	int m = vec.size();
	dp[0][0] = 1;
	for(int i = 0; i < m; i++){
		for(int j = 0; j < n; j++){
			for(int k = 0; k <= vec[i]; k++){
				if(j+k<n){
					dp[i+1][j+k] += dp[i][j] * comb(vec[i], k) % mod * comb(vec[i], k) % mod * fact[k] % mod;
					dp[i+1][j+k] %= mod;
				}
			}
		}
	}
	llint ans = fact[n-1], tmp = 0;
	for(int i = 1; i < n; i++){
		if(i % 2) tmp += fact[n-1-i] * dp[m][i] % mod;
		else tmp += mod - fact[n-1-i] * dp[m][i] % mod;
		tmp %= mod;
	}
	ans += mod - tmp, ans %= mod;
	
	vec.push_back(1);
	m = vec.size();
	for(int i = 0; i <= m; i++){
		for(int j = 0; j <= n; j++){
			dp[i][j] = 0;
		}
	}
	dp[0][0] = 1;
	for(int i = 0; i < m; i++){
		for(int j = 0; j <= n; j++){
			for(int k = 0; k <= vec[i]; k++){
				if(j+k<=n){
					dp[i+1][j+k] += dp[i][j] * comb(vec[i], k) % mod * comb(vec[i], k) % mod * fact[k] % mod;
					dp[i+1][j+k] %= mod;
				}
			}
		}
	}
	llint ans2 = fact[n]; tmp = 0;
	for(int i = 1; i <= n; i++){
		if(i % 2) tmp += fact[n-i] * dp[m][i] % mod;
		else tmp += mod - fact[n-i] * dp[m][i] % mod;
		tmp %= mod;
	}
	ans2 += mod - tmp, ans %= mod;
	
	cout << (ans+ans2)%mod << endl;
	return 0;
}
```