#include <bits/stdc++.h>
using namespace std;

typedef pair<int, int> pii;
typedef long long ll;
template<class T> using V = vector<T>;
template<class T> using VV = V<V<T>>;

#define pb push_back
#define eb emplace_back
#define mp make_pair
#define fi first
#define se second
#define rep(i,n) rep2(i,0,n)
#define rep2(i,m,n) for(int i=m;i<(n);i++)
#define ALL(c) (c).begin(),(c).end()
#define dump(x) cerr << #x << " = " << (x) << endl
constexpr ll TEN(int n) { return (n == 0) ? 1 : 10 * TEN(n-1); }

template<class T, class U> void chmin(T& t, const U& u) { if (t > u) t = u; }
template<class T, class U> void chmax(T& t, const U& u) { if (t < u) t = u; }

template<class T, class U>
ostream& operator<<(ostream& os, const pair<T, U>& p) {
	os<<"("<<p.first<<","<<p.second<<")";
	return os;
}

template<class T>
ostream& operator<<(ostream& os, const vector<T>& v) {
	os<<"{";
	rep(i, v.size()) {
		if (i) os<<",";
		os<<v[i];
	}
	os<<"}";
	return os;
}

const ll MOD = TEN(9) + 7;
inline void add(ll &x, ll y) { x += y; if (x >= MOD) x %= MOD; }

const int maxn = 5010;

int N;
V<int> g[maxn];
int sz[maxn], mc[maxn];
ll fac[maxn];
int par[maxn];
ll comb[maxn][maxn];

void dfs(int v, int p) {
	sz[v] = 1;
	par[v] = p;

	for (int to : g[v]) if (to != p) {
		dfs(to, v);
		sz[v] += sz[to];
		chmax(mc[v], sz[to]);
	}

	chmax(mc[v], N - sz[v]);
}

ll dp[maxn][maxn];

int main() {
	fac[0] = 1;
	rep(i, maxn - 1) {
		fac[i+1] = fac[i] * (i+1) % MOD;
	}

	rep(i, maxn) {
		comb[i][0] = comb[i][i] = 1;
		for (int j = 1; j < i; ++j) {
			comb[i][j] = (comb[i-1][j-1] + comb[i-1][j]) % MOD;
		} 
	}

	cin >> N;

	rep(i, N-1) {
		int a, b;
		cin >> a >> b;
		--a; --b;
		g[a].pb(b);
		g[b].pb(a);
	}

	dfs(0, -1);

	int r = -1;

	rep(i, N) {
		if (r == -1 || mc[r] > mc[i]) {
			r = i;
		}
		if (sz[i] * 2 == N) {
			cout << fac[N/2] * fac[N/2] % MOD << endl;
			return 0;
		}
	}

	V<int> num;

	for (int to : g[r]) {
		if (to != par[r]) {
			num.pb(sz[to]);
		} else {
			num.pb(N-sz[r]);
		}
	}

	dp[0][0] = 1;

	int sz = num.size();
	rep(i, sz) {
		int x = num[i];
		rep(j, maxn) if (dp[i][j]) {
			for (int k = 0; k <= x; ++k) {
				ll coef = comb[x][k] * comb[x][k] % MOD * fac[k] % MOD;
				ll val = dp[i][j] * coef % MOD;
				if (k % 2 == 1) {
					val = (MOD - val) % MOD;
				}
				add(dp[i+1][j+k], val);
			}
		}
	}

	ll ans = 0;
	for (int i = 0; i <= N-1; ++i) {
		ll val = dp[sz][i] * fac[N - i] % MOD;
		add(ans, val);
	}
	cout << ans << endl;
	return 0;
}
### Prompt

Construct a cpp code solution to the problem outlined:

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>
using namespace std;

typedef pair<int, int> pii;
typedef long long ll;
template<class T> using V = vector<T>;
template<class T> using VV = V<V<T>>;

#define pb push_back
#define eb emplace_back
#define mp make_pair
#define fi first
#define se second
#define rep(i,n) rep2(i,0,n)
#define rep2(i,m,n) for(int i=m;i<(n);i++)
#define ALL(c) (c).begin(),(c).end()
#define dump(x) cerr << #x << " = " << (x) << endl
constexpr ll TEN(int n) { return (n == 0) ? 1 : 10 * TEN(n-1); }

template<class T, class U> void chmin(T& t, const U& u) { if (t > u) t = u; }
template<class T, class U> void chmax(T& t, const U& u) { if (t < u) t = u; }

template<class T, class U>
ostream& operator<<(ostream& os, const pair<T, U>& p) {
	os<<"("<<p.first<<","<<p.second<<")";
	return os;
}

template<class T>
ostream& operator<<(ostream& os, const vector<T>& v) {
	os<<"{";
	rep(i, v.size()) {
		if (i) os<<",";
		os<<v[i];
	}
	os<<"}";
	return os;
}

const ll MOD = TEN(9) + 7;
inline void add(ll &x, ll y) { x += y; if (x >= MOD) x %= MOD; }

const int maxn = 5010;

int N;
V<int> g[maxn];
int sz[maxn], mc[maxn];
ll fac[maxn];
int par[maxn];
ll comb[maxn][maxn];

void dfs(int v, int p) {
	sz[v] = 1;
	par[v] = p;

	for (int to : g[v]) if (to != p) {
		dfs(to, v);
		sz[v] += sz[to];
		chmax(mc[v], sz[to]);
	}

	chmax(mc[v], N - sz[v]);
}

ll dp[maxn][maxn];

int main() {
	fac[0] = 1;
	rep(i, maxn - 1) {
		fac[i+1] = fac[i] * (i+1) % MOD;
	}

	rep(i, maxn) {
		comb[i][0] = comb[i][i] = 1;
		for (int j = 1; j < i; ++j) {
			comb[i][j] = (comb[i-1][j-1] + comb[i-1][j]) % MOD;
		} 
	}

	cin >> N;

	rep(i, N-1) {
		int a, b;
		cin >> a >> b;
		--a; --b;
		g[a].pb(b);
		g[b].pb(a);
	}

	dfs(0, -1);

	int r = -1;

	rep(i, N) {
		if (r == -1 || mc[r] > mc[i]) {
			r = i;
		}
		if (sz[i] * 2 == N) {
			cout << fac[N/2] * fac[N/2] % MOD << endl;
			return 0;
		}
	}

	V<int> num;

	for (int to : g[r]) {
		if (to != par[r]) {
			num.pb(sz[to]);
		} else {
			num.pb(N-sz[r]);
		}
	}

	dp[0][0] = 1;

	int sz = num.size();
	rep(i, sz) {
		int x = num[i];
		rep(j, maxn) if (dp[i][j]) {
			for (int k = 0; k <= x; ++k) {
				ll coef = comb[x][k] * comb[x][k] % MOD * fac[k] % MOD;
				ll val = dp[i][j] * coef % MOD;
				if (k % 2 == 1) {
					val = (MOD - val) % MOD;
				}
				add(dp[i+1][j+k], val);
			}
		}
	}

	ll ans = 0;
	for (int i = 0; i <= N-1; ++i) {
		ll val = dp[sz][i] * fac[N - i] % MOD;
		add(ans, val);
	}
	cout << ans << endl;
	return 0;
}
```