【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就是不信>_<。。。

Input

第一行一個整數n，表示樹的大小
接下來n-1行每行兩個數a,b，表示a和b之間有一條邊
註意點是從0開始標號的

Output

一行一個浮點數表示答案
四舍五入到小數點後4位
如果害怕精度跪建議用long double或者extended

Sample Input

Sample Output

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+期望