[noi31]MST

定义dp[i]表示当前连通块状态为i的方案数（状态记录该状态每一个连通块的大小），那么从小到大枚举每条边，考虑这条边在不在最小生成树上：

1. 如果不在最小生成树上，那么这条边有$\sum_{i=1}^{scc}\sum_{j=i+1}^{scc}Si\cdot Sj$（Si表示第i个连通块的点数）种位置，即每一个状态的方案都乘上这个数；

2. 如果在最小生成树上，那么他一定会把某两个连通块连起来，枚举这两个连通块并递推到新的状态。

那么这样的时间复杂度是多少呢？大约是$o(Pn\cdot n^{3})$（Pn表示将n划分成一些数的和的方案数，$P_{40}\approx 40000$，$n^{2}$是枚举连通块，另一个n是hash的时间），显然这样是会炸掉的……

似乎这个状态的记录方式可以改变，可改为每一个大小的连通块出现次数，由于最多只有$\sqrt{n}$种方案，枚举两个连通块也仅有n的时间复杂度，总时间复杂度降为$o(Pn\cdot n^{2})$

#include<bits/stdc++.h>
using namespace std;
#define mod 1000000007
map<int,int>vis;
int n,m,k,di[40001][41],a[41],mi[41],b[1001],f[40001];
vector<int>vec[40001];
void dfs(int k,int s,int p){
    if (!s){
        memcpy(di[++m],a,sizeof(a));
        di[m][0]=k;
    }
    if (s<1)return;
    k++;
    for(int i=1;i<=p;i++){
        a[i]++;
        dfs(k,s-i,i);
        a[i]--;
    }
}
int ha(int k){
    int ans=0;
    for(int i=1;i<=n;i++)ans=(ans+1LL*mi[i]*di[k][i])%mod;
    return ans;
}
int main(){
    mi[0]=1;
    for(int i=1;i<=40;i++)mi[i]=mi[i-1]*41LL%mod; 
    scanf("%d",&n);
    dfs(0,n,n);
    for(int i=1;i<=m;i++)vec[di[i][0]].push_back(i);
    for(int i=1;i<n;i++){
        scanf("%d",&k);
        b[k]=1;
    }
    f[1]=1;
    m=n*(n-1)/2;
    k=n;
    for(int i=1;i<=m;i++)
        if (!b[i])
            for(int j=0;j<vec[k].size();j++){
                int v=vec[k][j],s=m-i+1;
                for(int x=1;x<=n;x++)
                    for(int y=x+1;y<=n;y++)s-=di[v][x]*di[v][y]*x*y;
                for(int x=1;x<=n;x++)s-=(di[v][x]-1)*di[v][x]*x*x/2;
                f[vec[k][j]]=f[vec[k][j]]*1LL*s%mod; 
            }
        else{
            for(int j=0;j<vec[k-1].size();j++)vis[ha(vec[k-1][j])]=vec[k-1][j];
            for(int j=0;j<vec[k].size();j++){
                int v=vec[k][j];
                for(int x=1;x<=n;x++)
                    for(int y=x+1;y<=n-x;y++){
                        if ((!di[v][x])||(!di[v][y]))continue;
                        di[v][x]--;
                        di[v][y]--;
                        di[v][x+y]++;
                        f[vis[ha(v)]]=(f[vis[ha(v)]]+f[v]*(di[v][x]+1LL)*(di[v][y]+1)*x*y)%mod;
                        di[v][x]++;
                        di[v][y]++;
                        di[v][x+y]--;
                    }
                for(int x=1;x<=n/2;x++)
                    if (di[v][x]>1){
                        di[v][x]-=2;
                        di[v][x+x]++;
                        f[vis[ha(v)]]=(f[vis[ha(v)]]+f[v]*(di[v][x]+2LL)*(di[v][x]+1)*x*x/2)%mod;
                        di[v][x]+=2;
                        di[v][x+x]--;
                    }
            }
            k--;
            for(int j=0;j<vec[k].size();j++)vis[ha(vec[k][j])]=0;
        }
    printf("%d",f[vec[1][0]]);
}