【安徽集训】数树

• 并没有题目名称，“数树”是我瞎起的

Description

　　有一棵由 \(n\) 个点 \(n-1\) 条边构成的树，点从 \(0\) 到 \(n-1\) 标号。
　　你可以将它的一条边换成另外一条边，要求这样操作之后它还是一棵树。
　　求进行不超过 \(k\) 次操作后，能构造出多少种不同的树。
　　“两棵树不同”的定义：存在一棵树中有边 \((x,y)\) 而另一棵树中没有。
　　\(n\le 50\)
　　\(0\le k\le n\)

Solution

　　首先一个图的生成树个数 可以用矩阵树定理计算。
　　那对于这道题，我们要往基尔霍夫矩阵中某些位置代入自变量 \(x\)，使得该矩阵的行列式为一个 \(x\) 的 \(n\) 次多项式。
　　\(x\) 表示什么呢？观察问题，我们想表示新树与原树重合的边的数量，故构造这 \(n\) 个点的完全图，然后把原树上的每条边由 \(1\) 条增加到 \(x\) 条。不难发现每把一条边增加到 \(x\) 条 行列式就会 \(\times\space x\)，那么我们可以推出 行列式中带 \(x^i\) 表示有多少棵生成树与原树重合 \(i\) 条边。
　　最终我们只需要求行列式的 \(n-k-1\) 到 \(n-1\) 次项的系数和。
　　问题是怎么求这个多项式？
　　考虑拉格朗日插值，将 \(1\) 到 \(n\) 这 \(n\) 个整数作为 \(x\) 代入基尔霍夫矩阵，任意去掉该矩阵的一行一列后，高斯消元求该矩阵的行列式，就会得到一个整数结果 \(y\)。
　　有了多项式在 \(n\) 个不同的 \(x\) 处的取值后，暴力多项式乘法即可求出多项式每次项的系数。当然也可以整体做一遍多项式乘法，插每个值时独立除以一个 \(1\) 次多项式，这样这部分的时间复杂度可以降至 \(O(n^2)\)。
　　但瓶颈复杂度在前面的 \(n\) 次高斯消元，时间复杂度为 \(O(n^4)\)。

#include<bits/stdc++.h>
#define ll long long
#define N 51
#define mod 998244353
using namespace std;
inline int read(){
    int x=0; bool f=1; char c=getchar();
    for(;!isdigit(c);c=getchar()) if(c=='-') f=0;
    for(; isdigit(c);c=getchar()) x=(x<<3)+(x<<1)+(c^'0');
    if(f) return x;
    return 0-x;
}
int n,k,f[N],a[N][N],e[N][N],y[N],ans;
int Pow(int x, int y){
    int ret=1;
    while(y){
        if(y&1) ret=(ll)ret*x%mod;
        x=(ll)x*x%mod;
        y>>=1;
    }
    return ret;
}
struct Poly{
	int t,coe[N];
	Poly() {}
	Poly(int _t) {t=_t; for(int i=0; i<=t; ++i) coe[i]=0;}
	Poly(int _t, int c0) {t=_t, coe[0]=c0;}
	Poly(int _t, int c1, int c0) {t=_t, coe[1]=c1, coe[0]=c0;}
	Poly operator + (const Poly &x)const{
		int _t = max(t, x.t);
		Poly ret(_t);
		for(int i=0; i<=_t; ++i)
			ret.coe[i] = (coe[i] + x.coe[i]) % mod;
		return ret;
	}
	Poly operator * (const Poly& x)const{
		Poly ret(t+x.t);
		for(int i=0; i<=t; ++i)
			for(int j=0; j<=x.t; ++j)
				(ret.coe[i+j] += (ll)coe[i] * x.coe[j] % mod) %= mod;
		return ret;
	}
}P;
int main(){
    n=read(), k=read();
    for(int i=2; i<=n; ++i) f[i]=read()+1;
    for(int x=1; x<=n; ++x){
        for(int i=1; i<=n; ++i)
            for(int j=1; j<=n; ++j)
				e[i][j] = (i==j ? n-1 : mod-1);
        for(int i=2; i<=n; ++i) e[i][f[i]]=e[f[i]][i]=mod-x, e[i][i]=(e[i][i]+x-1)%mod, e[f[i]][f[i]]=(e[f[i]][f[i]]+x-1)%mod;
        /*
        for(int i=1; i<n; ++i){
        	for(int j=1; j<n; ++j) printf("%d ",e[i][j]); putchar('\n');
        }*/
		y[x]=1;
        //gaussian elimination
        for(int i=1; i<n; ++i){
        	if(e[i][i]==0){
	            for(int j=i+1; j<n; ++j) if(e[j][i]) {swap(e[i],e[j]), y[x]=mod-1; break;}
            	if(e[i][i]==0) {y[x]=0; break;}
            }
            y[x] = (ll)y[x] * e[i][i] % mod;
            int inv = Pow(e[i][i],mod-2); 
            for(int j=i+1; j<n; ++j){
                int t = (ll)e[j][i] * inv % mod;
                for(int k=i; k<n; ++k) e[j][k] = (e[j][k] - (ll)e[i][k] * t % mod + mod) % mod;
            }
        }
        //printf("%d\n",y[x]);
    }
    for(int i=1; i<=n; ++i){
    	//printf("%d %d\n",i,y[i]);
    	Poly tmp(0,y[i]); int div=1;
    	for(int j=1; j<=n; ++j) if(i!=j){
    		Poly tmp2(1,1,mod-j);
			tmp = tmp * tmp2;
			div = (ll)div * (i-j+mod) % mod;
		}
		Poly tmp2(0, Pow(div,mod-2));
		tmp = tmp * tmp2;
		P = P + tmp;
	}
	for(int i=n-k-1; i<n; ++i) ans=(ans+P.coe[i])%mod;
    cout<<ans<<endl;
    return 0;
}
/*
6 1
0 1 2 2 0
*/