【BZOJ3451】Tyvj1953 Normal 点分治+FFT+期望

【BZOJ3451】Tyvj1953 Normal

Description

某天WJMZBMR学习了一个神奇的算法：树的点分治！
这个算法的核心是这样的：
消耗时间=0
Solve(树 a)
 消耗时间 += a 的 大小
 如果 a 中 只有 1 个点
  退出
 否则在a中选一个点x，在a中删除点x
 那么a变成了几个小一点的树，对每个小树递归调用Solve
我们注意到的这个算法的时间复杂度跟选择的点x是密切相关的。
如果x是树的重心，那么时间复杂度就是O(nlogn)
但是由于WJMZBMR比较傻逼，他决定随机在a中选择一个点作为x！
Sevenkplus告诉他这样做的最坏复杂度是O(n^2)
但是WJMZBMR就是不信>_<。。。
于是Sevenkplus花了几分钟写了一个程序证明了这一点。。。你也试试看吧^_^
现在给你一颗树，你能告诉WJMZBMR他的傻逼算法需要的期望消耗时间吗？（消耗时间按在Solve里面的那个为标准）

Input

第一行一个整数n，表示树的大小
接下来n-1行每行两个数a,b，表示a和b之间有一条边
注意点是从0开始标号的

Output

一行一个浮点数表示答案
四舍五入到小数点后4位
如果害怕精度跪建议用long double或者extended

Sample Input

3
0 1
1 2

Sample Output

5.6667

HINT

n<=30000

题解：由于期望永远是可加的，所以我们可以讨论每个点对答案的贡献（即每个点在点分树上的深度）。对于x,y，我们统计y对x的贡献，即y成为x在点分树上的祖先的概率。y是x在点分树上的祖先当且仅当y是x-y路径上的第一个被选中的点。由于路径上每个点第一次被选中的概率都是相同的，所以概率就是1/dis(x,y)。具体地，我们的答案=$\sum\limits_{x=1}^n\sum\limits_{y=1}^n {1\over dis(x,y)}$。

所以我们希望对于任意的dis，统计出有多少点对之间的距离=dis，这个点分治+FFT即可。不过这里的点分治最好采用容斥的写法，即当以x为分治中心时，先统计出x子树中任意两点间的答案，再将两点再同一个儿子中的情况减去。

#include <cstdio>
#include <cstring>
#include <iostream>
#include <cmath>
#define pi acos(-1.0)
using namespace std;
const int maxn=30010;
struct cp
{
	double x,y;
	cp () {}
	cp (double a,double b){x=a,y=b;}
	cp operator + (const cp &a) const {return cp(x+a.x,y+a.y);}
	cp operator - (const cp &a) const {return cp(x-a.x,y-a.y);}
	cp operator * (const cp &a)	const {return cp(x*a.x-y*a.y,x*a.y+y*a.x);}
}A[maxn<<2];
long double Ans;
int n,rt,cnt,mx,tot,md;
int ans[maxn<<2];
int to[maxn<<1],next[maxn<<1],head[maxn],vis[maxn],siz[maxn],dep[maxn];
inline int rd()
{
	int ret=0,f=1;	char gc=getchar();
	while(gc<'0'||gc>'9')	{if(gc=='-')f=-f;	gc=getchar();}
	while(gc>='0'&&gc<='9')	ret=ret*10+gc-'0',gc=getchar();
	return ret*f;
}
void add(int a,int b)
{
	to[cnt]=b,next[cnt]=head[a],head[a]=cnt++;
}
void FFT(cp *a,int len,int f)
{
	int i,j,k,h;
	cp t;
	for(i=k=0;i<len;i++)
	{
		if(i>k)	swap(a[i],a[k]);
		for(j=len>>1;(k^=j)<j;j>>=1);
	}
	for(h=2;h<=len;h<<=1)
	{
		cp wn(cos(2*pi*f/h),sin(2*pi*f/h));
		for(j=0;j<len;j+=h)
		{
			cp w(1,0);
			for(k=j;k<j+h/2;k++)	t=a[k+h/2]*w,a[k+h/2]=a[k]-t,a[k]=a[k]+t,w=w*wn;
		}
	}
}
void getr(int x,int fa)
{
	siz[x]=1;
	int tmp=0;
	for(int i=head[x];i!=-1;i=next[i])
		if(to[i]!=fa&&!vis[to[i]])	getr(to[i],x),siz[x]+=siz[to[i]],tmp=max(tmp,siz[to[i]]);
	tmp=max(tmp,tot-siz[x]);
	if(tmp<mx)	mx=tmp,rt=x;
}
void getd(int x,int fa,int dep)
{
	A[dep].x+=1,md=max(md,dep);
	for(int i=head[x];i!=-1;i=next[i])	if(to[i]!=fa&&!vis[to[i]])	getd(to[i],x,dep+1);
}
void calc(int x,int f)
{
	md=0,getd(x,0,0);
	int i,len;
	for(len=1;len<=md*2;len<<=1);
	FFT(A,len,1);
	for(i=0;i<len;i++)	A[i]=A[i]*A[i];
	FFT(A,len,-1);
	if(f==1)	for(i=0;i<len;i++)	ans[i+1]+=int(A[i].x/len+0.1);
	else	for(i=0;i<len;i++)	ans[i+3]-=int(A[i].x/len+0.1);
	memset(A,0,sizeof(A[0])*len);
}
void dfs(int x)
{
	vis[x]=1;
	calc(x,1);
	for(int i=head[x];i!=-1;i=next[i])	if(!vis[to[i]])	calc(to[i],0),tot=siz[to[i]],mx=1<<30,getr(to[i],x),dfs(rt);
}
int main()
{
	n=rd();
	int i,a,b;
	memset(head,-1,sizeof(head));
	for(i=1;i<n;i++)	a=rd()+1,b=rd()+1,add(a,b),add(b,a);
	tot=n,mx=1<<30,getr(1,0),dfs(rt);
	for(i=1;i<=n;i++)
	{
		Ans+=(long double)ans[i]/i;
	}
	printf("%.4lf",(double)Ans);
	return 0;
}