●Joyoi Normal

题链:

http://www.joyoi.cn/problem/tyvj-1953
题解:


定义d(u,v)这个函数,满足:
d(u,v)=1,当且仅当在点分树中,u是v的祖先
d(u,v)=0,其它情况
对于一个确定的点分树来说,
$\sum d(u,v)$就是答案。
但是现在点分树未知,我们再来考虑这个函数的取值
d(u,v)=1时必须满足u是v的祖先,
那么由于现在是随机的,所以对于u到v的这条链上的点,
u必须是第一个被随机为重心的点才能使得在点分树中,u是v祖先。
而这个概率是1/(dis(u,v)+1),
所以d(u,v)的期望取值就是1*1/(dis(u,v)+1)
即现在的问题变为求:$ANS=\sum 1/(dis(u,v)+1)$

那么必然要求出所有点对间的距离,所以采用点分治+FFT的做法。
对于分治的每一个重心,把整个子树dfs一遍,得到所有节点到重心的距离,
并得到一个数组A[i],表示子树中到重心距离为i的点数为A[i]个
然后让A数组自乘,做FFT,可以得到改子树内任意两个点经过重心的距离,
显然对于在同一个儿子子树里的两个点,也被统计了进去,所以再对重心的每个儿子做上述操作,然后减掉无效信息。
复杂度$O(Nlog^2N)$


代码:

 

#include<bits/stdc++.h>
#define MAXN 65540
using namespace std;
const double Pi=acos(-1);
double ANS,TIME;
bool vis[MAXN];
int size[MAXN];
int N,cnt;
struct Edge{
	int ent;
	int to[MAXN*2],nxt[MAXN*2],head[MAXN];
	Edge(){ent=2;}
	void Adde(int u,int v){
		to[ent]=v; nxt[ent]=head[u]; head[u]=ent++;
	}
}E;
struct Cpx{
	double R,I;
	Cpx(){}
	Cpx(double _R,double _I):R(_R),I(_I){}
	Cpx operator - () const{return Cpx(-R,-I);}
	Cpx operator + (const Cpx &rtm) const{return Cpx(R+rtm.R,I+rtm.I);}
	Cpx operator - (const Cpx &rtm) const{return *this+(-rtm);}
	Cpx operator * (const Cpx &rtm) const{return Cpx(R*rtm.R-I*rtm.I,R*rtm.I+I*rtm.R);}
	Cpx operator / (const double k) const{return Cpx(R/k,I/k);}
}A[MAXN];
namespace FFT{
	int order[MAXN];
	int Init(int _n){
		static int n,len;
		for(n=1,len=0;n<=2*_n;n<<=1) len++;
		for(int i=0;i<n;i++) order[i]=(order[i>>1]>>1)|((i&1)<<(len-1));
		return n;
	}
	void Fft(Cpx *Y,int n,int sign){
	//	static double timebefore; timebefore=1.0*clock()/CLOCKS_PER_SEC;
		for(int i=0;i<n;i++) if(i<order[i]) swap(Y[i],Y[order[i]]);
		for(int d=2;d<=n;d<<=1){
			Cpx dw(cos(2*Pi/d),sin(sign*2*Pi/d)),w,tmp;
			for(int i=0;w=Cpx(1,0),i<n;i+=d)
				for(int k=i;k<i+d/2;w=w*dw,k++)
					tmp=w*Y[k+d/2],Y[k+d/2]=Y[k]-tmp,Y[k]=Y[k]+tmp;
		}
		if(sign==-1) for(int i=0;i<n;i++) Y[i]=Y[i]/n;
//		TIME+=1.0*clock()/CLOCKS_PER_SEC-timebefore;
	}
}
void getroot(int u,int dad,int num,int &root,int &rootmax){
	int umax=0; size[u]=0; cnt++;
	for(int e=E.head[u];e;e=E.nxt[e]){
		int v=E.to[e]; 
		if(v==dad||vis[v]) continue;
		getroot(v,u,num,root,rootmax);
		size[u]+=size[v];
		umax=max(umax,size[v]);
	}
	size[u]++; umax=max(umax,num-size[u]);
	if(rootmax>umax) root=u,rootmax=umax; 
}
void calc(int s,int dep,int sign){
	static queue<int> Q;
	static int dis[MAXN],reach[MAXN],tim,n;
	Q.push(s); reach[s]=++tim; n=0;
	dis[s]=dep; A[dis[s]].R+=1; n=dis[s];
	while(!Q.empty()){
		cnt++;
		int u=Q.front(); Q.pop();
		for(int e=E.head[u];e;e=E.nxt[e]){
			int v=E.to[e]; 
			if(reach[v]==tim||vis[v]) continue;
			reach[v]=tim; dis[v]=dis[u]+1; 
			n=max(n,dis[v]); A[dis[v]].R++;
			Q.push(v);
		}
	}
	n=FFT::Init(n);
	FFT::Fft(A,n,1);
	for(int i=0;i<n;i++) A[i]=A[i]*A[i];
	FFT::Fft(A,n,-1);
	for(int i=1;i<n;i++) 
		ANS+=sign*A[i].R*1.0/(i+1),A[i]=Cpx(0,0);
	ANS+=sign*A[0].R; A[0]=Cpx(0,0);
}
void solve(int u){
	if(size[u]!=N) calc(u,1,-1);
	int root=u,rootmax=size[u];
	getroot(u,0,size[u],root,rootmax);
	//cout<<rootmax<<endl;
	vis[root]=1;
	calc(root,0,1);
	for(int e=E.head[root];e;e=E.nxt[e]){
		int v=E.to[e]; if(vis[v]) continue;
		if(size[v]>size[root]) size[v]=size[u]-size[root];
		solve(v);
	}
}
int main(){
	ios::sync_with_stdio(0);
	cin>>N;
	for(int i=1,u,v;i<N;i++){
		cin>>u>>v; u++; v++;
		E.Adde(u,v); E.Adde(v,u);
	}
	size[1]=N; solve(1);
	cout<<fixed<<setprecision(4)<<ANS<<endl;
	//cout<<TIME<<" "<<cnt<<endl;	
	return 0;
}

  

 



Do not go gentle into that good night.
Rage, rage against the dying of the light.
————Dylan Thomas
posted @ 2018-03-11 08:58  *ZJ  阅读(246)  评论(0编辑  收藏  举报