#include<bits/stdc++.h>
#define X first
#define Y second
#define pb emplace_back
#define FOR(i,a,b) for(int (i)=(a);i<(b);++(i))
#define EFOR(i,a,b) for(int (i)=(a);i<=(b);++(i))
#define rep(X,Y) for (int (X) = 0;(X) < (Y);++(X))
#define REP rep
#define rrep(X,Y) for (int (X) = (Y)-1;(X) >=0;--(X))
#define all(X) (X).begin(),(X).end()
#define eb emplace_back

using namespace std;
typedef long long ll;
typedef pair<int,int> pii;
typedef pair<ll,ll> pll;
typedef ll LL;
typedef pii PII;
typedef pll PLL;
const ll MOD=1e9+7;

#define rall(X) (X).rbegin(),(X).rend()
#define UNIQUE(X) (X).erase(unique(all(X)),(X).end())
#define reps(X,S,Y) for (int (X) = S;(X) < (Y);++(X))
#define rreps(X,S,Y) for (int (X) = (Y)-1;(X) >= (S);--(X))
template<class T> inline bool MX(T &l,const T &r){return l<r?l=r,1:0;}
template<class T> inline bool MN(T &l,const T &r){return l>r?l=r,1:0;}

int N;
int sz[5145];
LL fact[5145];
LL rev[5145];
vector<int> cent;
vector<int> es[5145];
vector<LL> dp[5145];

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

  bool fail = false;
  for (int u : es[v]) {
    if (u == p) continue;
    Dfs(u, v);
    sz[v] += sz[u];
    if (N/2 < sz[u]) fail = true;
  }
  if (N/2 < N-sz[v]) fail = true;

  if (!fail) cent.eb(v);
}

LL mypow(LL base, LL exp) {
  if (!exp) return 1;
  LL res = mypow(base*base%MOD, exp/2);
  if (exp%2) res = res*base%MOD;
  return res;
}

LL Comb(LL n, LL k) {
  return fact[n] * rev[k] % MOD * rev[n-k] % MOD;
}

LL Solve(vector<int> &S) {
  int m = S.size();
  int sum = 0;
  for (int s : S) sum += s;

  dp[0].assign(1, 1);
  
  int lim = 0;
  rep(i, m) {
    dp[i+1].assign(dp[i].size() + S[i], 0);
    rep(j, dp[i].size()) {
      rep(t, S[i]+1) {
        LL c = Comb(S[i], t);
        c = c*c % MOD * fact[t] % MOD;
        dp[i+1][j+t] += c * dp[i][j];
        dp[i+1][j+t] %= MOD;
      }
    }
  }

  LL ans = fact[sum];
  reps(t, 1, dp[m].size()) {
    LL v = dp[m][t] * fact[sum - t] % MOD;
    if (t%2) {
      ans += MOD - v;
    } else {
      ans += v;
    }
    if (ans >= MOD) ans -= MOD;
  }
  //cout << "ret: " << ans << endl;
  return ans;
}

signed main(){
  ios_base::sync_with_stdio(false);
  cout<<fixed<<setprecision(0);

  cin >> N;
  fact[0] = rev[0] = 1;
  reps(i, 1, N+1) {
    fact[i] = fact[i-1] * i % MOD;
    rev[i] = rev[i-1] * mypow(i, MOD-2) % MOD;
  }

  rep(i, N-1) {
    int x, y;
    cin >> x >> y;
    --x;
    --y;
    es[x].eb(y);
    es[y].eb(x);
  }

  Dfs(0, 0);
  assert(!cent.empty() && cent.size() <= 2);
  vector<int> ss;
  LL ans;
  if (cent.size() == 2) {
    assert(N%2 == 0);
    ss.eb(N/2);
    ss.eb(N/2);
    ans = Solve(ss);
  } else {
    int r = cent[0];
    //cout << "cent: " << r << endl;
    fill(sz, sz+N, 0);
    Dfs(r, r);

    for (int v : es[r]) {
      ss.eb(sz[v]);
    }
    ss.eb(1);
    ans = Solve(ss);
    ss.pop_back();
    ans += Solve(ss);
    if (ans >= MOD) ans -= MOD;
  }
  cout << ans << endl;
}

### 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 X first
#define Y second
#define pb emplace_back
#define FOR(i,a,b) for(int (i)=(a);i<(b);++(i))
#define EFOR(i,a,b) for(int (i)=(a);i<=(b);++(i))
#define rep(X,Y) for (int (X) = 0;(X) < (Y);++(X))
#define REP rep
#define rrep(X,Y) for (int (X) = (Y)-1;(X) >=0;--(X))
#define all(X) (X).begin(),(X).end()
#define eb emplace_back

using namespace std;
typedef long long ll;
typedef pair<int,int> pii;
typedef pair<ll,ll> pll;
typedef ll LL;
typedef pii PII;
typedef pll PLL;
const ll MOD=1e9+7;

#define rall(X) (X).rbegin(),(X).rend()
#define UNIQUE(X) (X).erase(unique(all(X)),(X).end())
#define reps(X,S,Y) for (int (X) = S;(X) < (Y);++(X))
#define rreps(X,S,Y) for (int (X) = (Y)-1;(X) >= (S);--(X))
template<class T> inline bool MX(T &l,const T &r){return l<r?l=r,1:0;}
template<class T> inline bool MN(T &l,const T &r){return l>r?l=r,1:0;}

int N;
int sz[5145];
LL fact[5145];
LL rev[5145];
vector<int> cent;
vector<int> es[5145];
vector<LL> dp[5145];

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

  bool fail = false;
  for (int u : es[v]) {
    if (u == p) continue;
    Dfs(u, v);
    sz[v] += sz[u];
    if (N/2 < sz[u]) fail = true;
  }
  if (N/2 < N-sz[v]) fail = true;

  if (!fail) cent.eb(v);
}

LL mypow(LL base, LL exp) {
  if (!exp) return 1;
  LL res = mypow(base*base%MOD, exp/2);
  if (exp%2) res = res*base%MOD;
  return res;
}

LL Comb(LL n, LL k) {
  return fact[n] * rev[k] % MOD * rev[n-k] % MOD;
}

LL Solve(vector<int> &S) {
  int m = S.size();
  int sum = 0;
  for (int s : S) sum += s;

  dp[0].assign(1, 1);
  
  int lim = 0;
  rep(i, m) {
    dp[i+1].assign(dp[i].size() + S[i], 0);
    rep(j, dp[i].size()) {
      rep(t, S[i]+1) {
        LL c = Comb(S[i], t);
        c = c*c % MOD * fact[t] % MOD;
        dp[i+1][j+t] += c * dp[i][j];
        dp[i+1][j+t] %= MOD;
      }
    }
  }

  LL ans = fact[sum];
  reps(t, 1, dp[m].size()) {
    LL v = dp[m][t] * fact[sum - t] % MOD;
    if (t%2) {
      ans += MOD - v;
    } else {
      ans += v;
    }
    if (ans >= MOD) ans -= MOD;
  }
  //cout << "ret: " << ans << endl;
  return ans;
}

signed main(){
  ios_base::sync_with_stdio(false);
  cout<<fixed<<setprecision(0);

  cin >> N;
  fact[0] = rev[0] = 1;
  reps(i, 1, N+1) {
    fact[i] = fact[i-1] * i % MOD;
    rev[i] = rev[i-1] * mypow(i, MOD-2) % MOD;
  }

  rep(i, N-1) {
    int x, y;
    cin >> x >> y;
    --x;
    --y;
    es[x].eb(y);
    es[y].eb(x);
  }

  Dfs(0, 0);
  assert(!cent.empty() && cent.size() <= 2);
  vector<int> ss;
  LL ans;
  if (cent.size() == 2) {
    assert(N%2 == 0);
    ss.eb(N/2);
    ss.eb(N/2);
    ans = Solve(ss);
  } else {
    int r = cent[0];
    //cout << "cent: " << r << endl;
    fill(sz, sz+N, 0);
    Dfs(r, r);

    for (int v : es[r]) {
      ss.eb(sz[v]);
    }
    ss.eb(1);
    ans = Solve(ss);
    ss.pop_back();
    ans += Solve(ss);
    if (ans >= MOD) ans -= MOD;
  }
  cout << ans << endl;
}

```