[矩阵树定理 拉格朗日插值] TCO14 Round 3B TreeDistance

考虑两棵树，最少需要多少步把其中一棵变成另一棵？
可以发现答案就是不同的边数。即存在于A中而不存在于B中的边
数。
相当于要求：有多少带标号无根树，只有≤ $k$ 条原树中没有的边？
把原树中的边设成$1$，不在原树中的边设成$x$，跑矩阵树定理，得出
来是一个$n$阶多项式。
$x_i(i ≤ k)$前的系数和就是答案

// BEGIN CUT HERE  
#include<conio.h>
#include<sstream>
// END CUT HERE  
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<iostream>
#include<algorithm>
#include<vector>
#include<map>
#include<string>
#include<set>
#define cl(x) memset(x,0,sizeof(x))
using namespace std;
typedef long long ll;


const int P=1e9+7;
const int N=85;

inline ll Pow(ll a,int b){
  ll ret=1; for (;b;b>>=1,a=a*a%P) if(b&1) ret=ret*a%P; return ret;
}
inline ll Inv(ll a){
  return Pow(a,P-2);
}

int n,K;
int fat[N];
int y[N];
int a[N][N];

inline int det(int n){
  int t=0;
  for (int i=1;i<=n;i++){
    int k=0;
    for (int j=i;j<=n;j++) if (a[j][i]) { k=j; break; }
    if (i^k) { for (int j=i;j<=n;j++) swap(a[i][j],a[k][j]); t^=1; }
    for (int j=i+1;j<=n;j++){
      ll t=(ll)Inv(a[i][i])*a[j][i]%P;
      for (int k=i;k<=n;k++)
    (a[j][k]+=P-(ll)t*a[i][k]%P)%=P;
    }
  }
  ll ret=1;
  for (int i=1;i<=n;i++) ret=ret*a[i][i]%P;
  return t?(P-ret)%P:ret;
}

inline int calc(int x){
  for (int i=1;i<=n;i++) for (int j=1;j<=n;j++) a[i][j]=0;
  for (int i=1;i<=n;i++)
    for (int j=1;j<i;j++){
      if (i==fat[j] || j==fat[i])
    a[i][j]=a[j][i]=P-1,a[i][i]+=1,a[j][j]+=1;
      else
    a[i][j]=a[j][i]=P-x,a[i][i]+=x,a[j][j]+=x;
    }
  return det(n-1);
}

int f[N],g[N];

class TreeDistance{

public:
  int countTrees(vector <int> p, int _K){
    n=p.size()+1; K=_K;
    for (int i=2;i<=n;i++) fat[i]=p[i-2]+1;
    for (int i=1;i<=n;i++) y[i]=calc(i);
    cl(f); cl(g);
    for (int i=1;i<=n;i++){
      ll t=y[i];
      for (int k=0;k<n;k++) g[k]=0; g[0]=1;
      for (int j=1;j<=n;j++){
    if (i==j) continue; t=t*Inv((i+P-j)%P)%P;
    for (int k=n-1;~k;k--)
      g[k]=((ll)g[k]*(P-j)%P+(k?g[k-1]:0))%P;
      }
      for (int k=0;k<n;k++) (f[k]+=(ll)g[k]*t%P)%=P;
    }
    ll ans=0;
    for (int k=0;k<=K;k++) ans+=f[k];
    return (int)(ans%P);
  }


  // BEGIN CUT HERE
public:
  void run_test(int Case) { if ((Case == -1) || (Case == 0)) test_case_0(); if ((Case == -1) || (Case == 1)) test_case_1(); if ((Case == -1) || (Case == 2)) test_case_2(); if ((Case == -1) || (Case == 3)) test_case_3(); if ((Case == -1) || (Case == 4)) test_case_4(); if ((Case == -1) || (Case == 5)) test_case_5(); if ((Case == -1) || (Case == 6)) test_case_6(); }
private:
  template <typename T> string print_array(const vector<T> &V) { ostringstream os; os << "{ "; for (typename vector<T>::const_iterator iter = V.begin(); iter != V.end(); ++iter) os << '\"' << *iter << "\","; os << " }"; return os.str(); }
  void verify_case(int Case, const int &Expected, const int &Received) { cerr << "Test Case #" << Case << "..."; if (Expected == Received) cerr << "PASSED" << endl; else { cerr << "FAILED" << endl; cerr << "\tExpected: \"" << Expected << '\"' << endl; cerr << "\tReceived: \"" << Received << '\"' << endl; } }
  void test_case_0() { int Arr0[] = {0, 0}; vector <int> Arg0(Arr0, Arr0 + (sizeof(Arr0) / sizeof(Arr0[0]))); int Arg1 = 1; int Arg2 = 3; verify_case(0, Arg2, countTrees(Arg0, Arg1)); }
  void test_case_1() { int Arr0[] = {0, 1, 2, 2, 0}; vector <int> Arg0(Arr0, Arr0 + (sizeof(Arr0) / sizeof(Arr0[0]))); int Arg1 = 1; int Arg2 = 28; verify_case(1, Arg2, countTrees(Arg0, Arg1)); }
  void test_case_2() { int Arr0[] = {0, 1, 2, 2, 0}; vector <int> Arg0(Arr0, Arr0 + (sizeof(Arr0) / sizeof(Arr0[0]))); int Arg1 = 2; int Arg2 = 222; verify_case(2, Arg2, countTrees(Arg0, Arg1)); }
  void test_case_3() { int Arr0[] = {0, 1, 2, 2, 0}; vector <int> Arg0(Arr0, Arr0 + (sizeof(Arr0) / sizeof(Arr0[0]))); int Arg1 = 50; int Arg2 = 1296; verify_case(3, Arg2, countTrees(Arg0, Arg1)); }
  void test_case_4() { int Arr0[] = {0, 1, 2, 2, 0}; vector <int> Arg0(Arr0, Arr0 + (sizeof(Arr0) / sizeof(Arr0[0]))); int Arg1 = 0; int Arg2 = 1; verify_case(4, Arg2, countTrees(Arg0, Arg1)); }
  void test_case_5() { int Arr0[] = {0, 1, 0, 3, 3, 4, 4, 5, 6, 8, 3, 1, 12, 12, 13, 10, 4, 8, 13, 17, 2, 10, 12, 20, 2, 14, 17, 19, 15, 0, 22, 15, 3, 8, 3, 17, 27, 2, 12, 38, 37, 4, 40, 29, 9, 22, 43, 32, 37}; vector <int> Arg0(Arr0, Arr0 + (sizeof(Arr0) / sizeof(Arr0[0]))); int Arg1 = 1; int Arg2 = 7124; verify_case(5, Arg2, countTrees(Arg0, Arg1)); }
  void test_case_6() { int Arr0[] = {0, 0, 0, 0, 2, 3, 1, 2, 3, 7, 3, 10, 8, 8, 9, 1, 2, 0, 7, 17, 19, 2, 17, 2, 0,
                     6, 4, 9, 12, 14, 8, 12, 10, 30, 20, 30, 8, 36, 28, 22, 8, 2, 2, 13, 26, 14, 46, 6, 25}; vector <int> Arg0(Arr0, Arr0 + (sizeof(Arr0) / sizeof(Arr0[0]))); int Arg1 = 10; int Arg2 = 310259667; verify_case(6, Arg2, countTrees(Arg0, Arg1)); }

  // END CUT HERE

};

// BEGIN CUT HERE
int main(){
  TreeDistance ___test;
  ___test.run_test(-1);
  getch() ;
  return 0;
}
// END CUT HERE