[bzoj3451]Tyvj1953 Normal——点分治+fft

题目大意：

求随机点分治的期望复杂度，每次对一颗大小为(n)的子树需要(O(n))的复杂度。

思路：

考虑计算每个点期望下被算的次数，根据期望的线性性，最后将每个点的答案加起来就可以了。
计算点u的计算次数可以考虑v对点u的贡献，即在v作为分治重心的时候u在v所在的子树里面。
不难发现如果v对u产生了贡献，那么从u到v的路径上，v必定是第一个选的，路径外的点怎么选没有影响，于是期望贡献为(frac{1}{dis(u,v)+1})。
答案即(sum_{i=1}^{n}sum_{j=1}^{n}frac{1}{dis(i,j)+1})，又转化成了树上路径问题，考虑点分治，计算出每颗子树内的所有路径长度的出现次数，直接FFT优化即可。

/*=======================================
 * Author : ylsoi
 * Time : 2019.2.13
 * Problem : bzoj3451
 * E-mail : ylsoi@foxmail.com
 * ====================================*/
#include<bits/stdc++.h>

#define REP(i,a,b) for(int i=a,i##_end_=b;i<=i##_end_;++i)
#define DREP(i,a,b) for(int i=a,i##_end_=b;i>=i##_end_;--i)
#define debug(x) cout<<#x<<"="<<x<<" "
#define fi first
#define se second
#define mk make_pair
#define pb push_back
typedef long long ll;

using namespace std;

void File(){
	freopen("bzoj3451.in","r",stdin);
	freopen("bzoj3451.out","w",stdout);
}

template<typename T>void read(T &_){
	_=0; T f=1; char c=getchar();
	for(;!isdigit(c);c=getchar())if(c=='-')f=-1;
	for(;isdigit(c);c=getchar())_=(_<<1)+(_<<3)+(c^'0');
	_*=f;
}

const int maxn=3e4+10;
const int inf=0x3f3f3f3f;
const double pi=acos(-1);
int n;
int beg[maxn],las[maxn<<1],to[maxn<<1],cnte=1;
double ans;

void add(int u,int v){
	las[++cnte]=beg[u],beg[u]=cnte,to[cnte]=v;
	las[++cnte]=beg[v],beg[v]=cnte,to[cnte]=u;
}

struct cp{
	double x,y;
	cp(double xx=0,double yy=0){
		x=xx,y=yy;
	}
};
cp operator + (cp a,cp b){return cp(a.x+b.x,a.y+b.y);}
cp operator - (cp a,cp b){return cp(a.x-b.x,a.y-b.y);}
cp operator * (cp a,cp b){return cp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
cp operator / (cp a,double b){return cp(a.x/b,a.y/b);}

int lim,cnt,dn[maxn<<2];
cp g[maxn<<2],ig[maxn<<2];

void fft(cp *A,int ty){
	REP(i,0,lim-1)if(i<dn[i])swap(A[i],A[dn[i]]);
	for(int len=1;len<lim;len<<=1){
		cp w= ty==1 ? g[len<<1] : ig[len<<1];
		for(int L=0;L<lim;L+=len<<1){
			cp wk=cp(1,0);
			REP(i,L,L+len-1){
				cp u=A[i],v=A[i+len]*wk;
				A[i]=u+v;
				A[i+len]=u-v;
				wk=wk*w;
			}
		}
	}
	if(ty==-1)
		REP(i,0,lim-1)A[i]=A[i]/lim;
}

int sz[maxn],tot_sz,Min_sz,rt;
bool vis[maxn];

void findrt(int u,int fh){
	int Max_sz=0;
	sz[u]=1;
	for(int i=beg[u];i;i=las[i]){
		int v=to[i];
		if(v==fh || vis[v])continue;
		findrt(v,u);
		sz[u]+=sz[v];
		Max_sz=max(Max_sz,sz[v]);
	}
	Max_sz=max(Max_sz,tot_sz-sz[u]);
	if(Max_sz<Min_sz){
		Min_sz=Max_sz;
		rt=u;
	}
}

int dis[maxn],cnt_dis;

void get_dis(int u,int fh,int d){
	dis[++cnt_dis]=d;
	for(int i=beg[u];i;i=las[i]){
		int v=to[i];
		if(v==fh || vis[v])continue;
		get_dis(v,u,d+1);
	}
}

cp a[maxn<<2];

void solve(int u,int s,int ty){
	cnt_dis=0;
	get_dis(u,0,0);
	int max_dis=0;
	REP(i,1,cnt_dis)max_dis=max(max_dis,dis[i]);
	lim=1,cnt=0;
	while(lim<=max_dis*2)lim<<=1,++cnt;
	if(!cnt)cnt=1;
	REP(i,0,lim-1){
		dn[i]=dn[i>>1]>>1|((i&1)<<(cnt-1));
		a[i]=cp(0,0);
	}
	REP(i,1,cnt_dis)a[dis[i]].x+=1;
	fft(a,1);
	REP(i,0,lim-1)a[i]=a[i]*a[i];
	fft(a,-1);
	REP(i,0,lim-1)if(i+s>0)
		ans=ans+a[i].x/(i+s)*ty;
}

void divide(int u){
	vis[u]=1; findrt(u,0);
	solve(u,1,1);
	for(int i=beg[u];i;i=las[i]){
		int v=to[i];
		if(vis[v])continue;
		solve(v,3,-1);
		tot_sz=sz[v],Min_sz=inf;
		findrt(v,0);
		divide(rt);
	}
}

int main(){
	File();
	read(n);
	int u,v;
	REP(i,1,n-1){
		read(u),read(v);
		add(u+1,v+1);
	}

	lim=1;
	while(lim<=n+n)lim<<=1;
	g[lim]=cp(cos(pi*2.0/lim),sin(pi*2.0/lim));
	ig[lim]=cp(cos(pi*2.0/lim),-sin(pi*2.0/lim));
	for(int i=lim>>1;i;i>>=1){
		g[i]=g[i<<1]*g[i<<1];
		ig[i]=ig[i<<1]*ig[i<<1];
	}

	tot_sz=n,Min_sz=inf;
	findrt(1,0);
	divide(rt);

	printf("%.4lf
",round(ans*1e4)/1e4);

	return 0;
}