NOIp 模拟赛 T3，AT 3200 分。还是很可做的！
考虑一条边 (u,v) 对答案贡献的上界：即 2\min(siz_u,siz_v)，如果固定根的话就是 2\min(siz_u,n-siz_u)。这实际上是很多问题的经典思路。
每一条边尝试能否取到上界，实际上是和重心相关的经典思路。先把重心拉成根，那么需要每个 y\in T(x) 都是往子树外连的，即 p_y\not\in T(x)。容易发现都可以取到上界。
进一步，发现我们只需对重心的儿子（称为二级根）考虑。那么问题转换成了：有多少 \{p_i\}，使得 (i,p_i) 不在二级根的同一子树。
考虑容斥，限制为：i,p_i 在根的同一子树。设 g_{u,k} 为二级根 u 的子树内钦定有 k 个限制，f_{rt,k} 表示重心内钦定有 k 个限制。最终答案要恰好 0 个限制，那么二项式反演即可求出答案。
ans=\sum (-1)^k f_{rt,k} (n-k)!
而重心 $rt$ 的答案只需要做一次树形背包合并即可。
时空复杂度均为 $O(n^2)$。
```cpp
#include <bits/stdc++.h>
#define int long long
#define pii pair<int,int>
#define fi first
#define se second
using namespace std;
int read(){
char c=getchar();int h=0,tag=1;
while(!isdigit(c)) tag=(c=='-'?-1:1),c=getchar();
while(isdigit(c)) h=(h<<1)+(h<<3)+(c^48),c=getchar();
return h*tag;
}
void fil(){
freopen("data.in","r",stdin);
freopen("data.out","w",stdout);
}
int rt,n;
const int N=5010;
int siz[N];vector<int>s[N];
void findrt(int x,int fa) {
siz[x]=1;
int maxson=-1;
for(int y:s[x]) {
if(y==fa) continue;
findrt(y,x);
siz[x]+=siz[y];
maxson=max(maxson,siz[y]);
}
int tmp=max(n-siz[x],maxson);
if(tmp<=n/2) rt=x;
}
void calcsiz(int x,int fa) {
siz[x]=1;
for(int y:s[x]) {
if(y==fa) continue;
calcsiz(y,x);siz[x]+=siz[y];
}
}
int fac[N],inv[N];
const int mod=1e9+7;
int ksm(int a,int b) {
if(b==1) return a;
int s=ksm(a,b/2);
s=s*s%mod;
if(b%2==1) s=s*a%mod;
return s;
}
int binom(int n,int m) {
return fac[n]*inv[m]%mod*inv[n-m]%mod;
}
int g[N][N],tmp[N],dp[N];
int sgn(int x) {
if(x%2==1) return mod-1;
else return 1;
}
signed main(){
//fil();
n=read();
fac[0]=inv[0]=1;
for(int i=1;i<=n;i++) fac[i]=fac[i-1]*i%mod,inv[i]=ksm(fac[i],mod-2)%mod;
for(int i=1;i<n;i++) {
int u=read(),v=read();
s[u].push_back(v);s[v].push_back(u);
}
findrt(1,0);
calcsiz(rt,0);
vector<int>vec;
for(int y:s[rt]) {
vec.push_back(y);
}
for(int i=0;i<vec.size();i++) {
int v=vec[i];
for(int k=0;k<=siz[v];k++) g[v][k]=binom(siz[v],k)*binom(siz[v],k)%mod*fac[k]%mod;
}
int up=1;
dp[0]=1;
for(int i=0;i<vec.size();i++) {
int v=vec[i];int nup=up+siz[v];
for(int j=0;j<=up;j++) {
for(int k=0;k<=siz[v]&&j+k<=nup;k++) {
tmp[j+k]=(tmp[j+k]+g[v][k]*dp[j]%mod)%mod;
}
}
for(int j=0;j<=nup;j++) dp[j]=tmp[j],tmp[j]=0;
up=nup;
}
int ans=0;
for(int i=0;i<=n;i++) ans=(ans+sgn(i)*dp[i]%mod*fac[n-i]%mod)%mod;
cout<<ans<<'\n';
return 0;
}
```