【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;

}

巴特西

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