【安徽集訓】數樹

• 並無題目名稱，「數樹」是我瞎起的

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
*/