#include <bits/stdc++.h>

#define _overload(_1,_2,_3,name,...) name
#define _rep(i,n) _range(i,0,n)
#define _range(i,a,b) for(int i=int(a);i<int(b);++i)
#define rep(...) _overload(__VA_ARGS__,_range,_rep,)(__VA_ARGS__)

#define _rrep(i,n) _rrange(i,n,0)
#define _rrange(i,a,b) for(int i=int(a)-1;i>=int(b);--i)
#define rrep(...) _overload(__VA_ARGS__,_rrange,_rrep,)(__VA_ARGS__)

#define _all(arg) begin(arg),end(arg)
#define uniq(arg) sort(_all(arg)),(arg).erase(unique(_all(arg)),end(arg))
#define getidx(ary,key) lower_bound(_all(ary),key)-begin(ary)
#define clr(a,b) memset((a),(b),sizeof(a))
#define bit(n) (1LL<<(n))
#define popcount(n) (__builtin_popcountll(n))

using namespace std;

template<class T>bool chmax(T &a, const T &b) { return (a < b) ? (a = b, 1) : 0;}
template<class T>bool chmin(T &a, const T &b) { return (b < a) ? (a = b, 1) : 0;}

using ll = long long;
using R = long double;
const R EPS = 1e-9L; // [-1000,1000]->EPS=1e-8 [-10000,10000]->EPS=1e-7
inline int sgn(const R& r) {return (r > EPS) - (r < -EPS);}
inline R sq(R x) {return sqrt(max(x, 0.0L));}

const int dx[8] = {1, 0, -1, 0, 1, -1, -1, 1};
const int dy[8] = {0, 1, 0, -1, 1, 1, -1, -1};

// Problem Specific Parameter:

const ll mod = 1000000007LL;
const ll limit = 5000;
ll F[limit + 10];
ll C[limit + 10][limit + 10];

ll dp[limit + 10][limit + 10];

void init() {
	F[0] = 1LL;
	rep(i, 1, limit + 1) F[i] = 1LL * i *  F[i - 1] % mod;

	rep(i, limit + 1) {
		rep(j, i + 1) {
			if (j == 0 or j == i)
				C[i][j] = 1LL;
			else
				C[i][j] = (C[i - 1][j - 1] + C[i - 1][j]) % mod;
		}
	}
}

vector<int> graph[limit];
int sz[limit];
vector<int> centroid;

int dfs(int v, int p, int total) {
	sz[v] = 1;
	const int thr = total / 2;
	bool ok = true;

	for (auto &v2 : graph[v]) {
		if (v2 == p) continue;
		sz[v] += dfs(v2, v, total);
		if (sz[v2] > thr) ok = false;
	}

	if (total - sz[v] > thr) ok = false;
	if (ok) centroid.push_back(v);

	return sz[v];
}

int main(void) {
	init();

	int n;
	cin >> n;
	rep(i, n - 1) {
		int a, b;
		cin >> a >> b;
		a--, b--;
		graph[a].push_back(b);
		graph[b].push_back(a);
	}

	dfs(0, -1, n);

	if (centroid.size() == 2) {
		ll ans = F[n / 2] * F[n / 2] % mod;
		cout << ans << endl;
		return 0;
	}

	const int c = centroid[0];
	dfs(c, -1, n);

	vector<int> ary;
	for (auto &v : graph[c]) ary.push_back(sz[v]);
	const int m = ary.size();
	dp[0][0] = 1LL;
	int total = 0;

	rep(i, m) {
		const int ni = ary[i];
		rep(j, total + 1) {
			if (dp[i][j] == 0) continue;
			rep(k, ni + 1) {
				const ll add = C[ni][k] * C[ni][k] % mod * F[k] % mod * dp[i][j] % mod;
				dp[i + 1][j + k] = (dp[i + 1][j + k] + add) % mod;
			}
		}
		total += ni;
	}

	rep(i, total)  dp[m][i] = dp[m][i] * F[total + 1 - i] % mod;


	ll ans = 0LL;
	rep(i, total + 1) {
		if (i % 2 == 1)
			ans = (ans + mod - dp[m][i]) % mod;
		else
			ans = (ans + dp[m][i]) % mod;
	}

	cout << ans << endl;
	return 0;
}
### Prompt

Please create a solution in Cpp 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 <bits/stdc++.h>

#define _overload(_1,_2,_3,name,...) name
#define _rep(i,n) _range(i,0,n)
#define _range(i,a,b) for(int i=int(a);i<int(b);++i)
#define rep(...) _overload(__VA_ARGS__,_range,_rep,)(__VA_ARGS__)

#define _rrep(i,n) _rrange(i,n,0)
#define _rrange(i,a,b) for(int i=int(a)-1;i>=int(b);--i)
#define rrep(...) _overload(__VA_ARGS__,_rrange,_rrep,)(__VA_ARGS__)

#define _all(arg) begin(arg),end(arg)
#define uniq(arg) sort(_all(arg)),(arg).erase(unique(_all(arg)),end(arg))
#define getidx(ary,key) lower_bound(_all(ary),key)-begin(ary)
#define clr(a,b) memset((a),(b),sizeof(a))
#define bit(n) (1LL<<(n))
#define popcount(n) (__builtin_popcountll(n))

using namespace std;

template<class T>bool chmax(T &a, const T &b) { return (a < b) ? (a = b, 1) : 0;}
template<class T>bool chmin(T &a, const T &b) { return (b < a) ? (a = b, 1) : 0;}

using ll = long long;
using R = long double;
const R EPS = 1e-9L; // [-1000,1000]->EPS=1e-8 [-10000,10000]->EPS=1e-7
inline int sgn(const R& r) {return (r > EPS) - (r < -EPS);}
inline R sq(R x) {return sqrt(max(x, 0.0L));}

const int dx[8] = {1, 0, -1, 0, 1, -1, -1, 1};
const int dy[8] = {0, 1, 0, -1, 1, 1, -1, -1};

// Problem Specific Parameter:

const ll mod = 1000000007LL;
const ll limit = 5000;
ll F[limit + 10];
ll C[limit + 10][limit + 10];

ll dp[limit + 10][limit + 10];

void init() {
	F[0] = 1LL;
	rep(i, 1, limit + 1) F[i] = 1LL * i *  F[i - 1] % mod;

	rep(i, limit + 1) {
		rep(j, i + 1) {
			if (j == 0 or j == i)
				C[i][j] = 1LL;
			else
				C[i][j] = (C[i - 1][j - 1] + C[i - 1][j]) % mod;
		}
	}
}

vector<int> graph[limit];
int sz[limit];
vector<int> centroid;

int dfs(int v, int p, int total) {
	sz[v] = 1;
	const int thr = total / 2;
	bool ok = true;

	for (auto &v2 : graph[v]) {
		if (v2 == p) continue;
		sz[v] += dfs(v2, v, total);
		if (sz[v2] > thr) ok = false;
	}

	if (total - sz[v] > thr) ok = false;
	if (ok) centroid.push_back(v);

	return sz[v];
}

int main(void) {
	init();

	int n;
	cin >> n;
	rep(i, n - 1) {
		int a, b;
		cin >> a >> b;
		a--, b--;
		graph[a].push_back(b);
		graph[b].push_back(a);
	}

	dfs(0, -1, n);

	if (centroid.size() == 2) {
		ll ans = F[n / 2] * F[n / 2] % mod;
		cout << ans << endl;
		return 0;
	}

	const int c = centroid[0];
	dfs(c, -1, n);

	vector<int> ary;
	for (auto &v : graph[c]) ary.push_back(sz[v]);
	const int m = ary.size();
	dp[0][0] = 1LL;
	int total = 0;

	rep(i, m) {
		const int ni = ary[i];
		rep(j, total + 1) {
			if (dp[i][j] == 0) continue;
			rep(k, ni + 1) {
				const ll add = C[ni][k] * C[ni][k] % mod * F[k] % mod * dp[i][j] % mod;
				dp[i + 1][j + k] = (dp[i + 1][j + k] + add) % mod;
			}
		}
		total += ni;
	}

	rep(i, total)  dp[m][i] = dp[m][i] * F[total + 1 - i] % mod;


	ll ans = 0LL;
	rep(i, total + 1) {
		if (i % 2 == 1)
			ans = (ans + mod - dp[m][i]) % mod;
		else
			ans = (ans + dp[m][i]) % mod;
	}

	cout << ans << endl;
	return 0;
}
```