TopCoder SRM 712 Div1 600 AverageVarianceSubtree

卡精度。。。。
今天才知道有__float128这种东西，问了一下noip不能用。。。

$Description$

传送门
$n$个点的树，每个点有一个权值，任选一棵非空子树，求其中所有点方差的期望。

$Solution$

裸的树形DP，感觉不值600，不过之前做过另一道树形DP也是600…
反正这种题吧，大概xy肯定会说要树转二叉树（并不用）

我们先来尻驴一下方差是个啥

$$S=\frac{(x_1-\bar{x}{)}^2+(x_1-\bar{x}{)}^2+\cdots +(x_n-\bar{x}{)}^2}{n}$$

emmm之前做的方差题里好像很多都有一句类似"可以证明$S\times n^2$是整数，然后整数比较好处理我就尻驴把它变整数。为了方便就让$s_1={\sum }_{i=1}^n{x}_i$，$s_2={\sum }_{i=1}^n{x}_i^2$，把式子一通展开就得到了：

$$S\cdot n^2=n\cdot s_2-s_1^2$$

这个式子一看就好看多了（并没有）

然后就开始树形DP的常规操作…（以下“所有子树”表示以$i$为根，子树大小为$j$的子树）

• $f[i][j]$ 表示方案数
• $h[i][j]$ 表示所有子树中$s_1$的和
• $s[i][j]$ 表示所有子树中$s_1^2$的和
• $s2[i][j]$ 表示所有子树中$s_2$的和

转移的话尽量自己思考吧，不难的（其实是不负责任的博主懒得写）。实在不懂就留言吧，我尽量回答。
其实感觉这种题很多就是暴力展开式子，然后把长得像的丢一起…好像做个几道就能get到其中的套路…

#include <bits/stdc++.h>
#define ll long long
#define fr(i,x,y) for(int i=x;i<=y;i++)
#define rf(i,x,y) for(int i=x;i>=y;i--)
//#define double long double
using namespace std;
const int N=55;
const int M=N<<1;
int n,w[N];
int cnt,head[N],Next[M],v[M];
__float128 f[N][N],s[N][N],s2[N][N],h[N][N];

template<class T> void checkmin(T &a,const T &b) { if (b<a) a=b; } 
template<class T> void checkmax(T &a,const T &b) { if (b>a) a=b; }

class AverageVarianceSubtree {
public:
    double average( vector <int> p, vector <int> weight ) ;
};

void add(int x,int y){
	Next[++cnt]=head[x];
	head[x]=cnt;
	v[cnt]=y;
}

ll sqr(ll x){
	return x*x;
}

void dfs(int x,int fa){
	f[x][1]=1;
	s[x][1]=sqr(w[x]);
	s2[x][1]=sqr(w[x]);
	h[x][1]=w[x];
	for(int i=head[x];i;i=Next[i]){
		dfs(v[i],x);
		//f[x]+=f[x]*f[v[i]];
		rf(j,n,1)
		 rf(k,n,1)
		  if (j+k<=n){
		  	f[x][j+k]+=f[x][j]*f[v[i]][k];
		  	s[x][j+k]+=f[v[i]][k]*s[x][j]+f[x][j]*s[v[i]][k]+2*h[x][j]*h[v[i]][k];
		  	s2[x][j+k]+=s2[x][j]*f[v[i]][k]+s2[v[i]][k]*f[x][j];
		  	h[x][j+k]+=h[x][j]*f[v[i]][k]+h[v[i]][k]*f[x][j];
		  }
	}
}

double AverageVarianceSubtree::average(vector <int> p, vector <int> weight) {
    n=weight.size();
    fr(i,2,n) add(p[i-2]+1,i);
    fr(i,1,n) w[i]=weight[i-1];
    dfs(1,0);
    __float128 sum=0,z=0;
    fr(i,1,n){
    	__float128 ss=0,ss2=0;
    	fr(j,1,n){
    		ss+=s[j][i];
    		ss2+=s2[j][i];
    		sum+=(s2[j][i]*i-s[j][i])/(i*i);
    		z+=f[j][i];
    	}
    }
    double ans=sum/z;
    return ans;
}