[BZOJ3451] Tyvj1953 Normal 点分治+FFT

题目链接

显然按照点的贡献来考虑,答案就是每个点在点分树的期望深度之和。

深度一般可以转化为它是几个点的子树。

考虑\(y\)\(x\)的点分树的子树中。那么需要保证在原树\(x\)\(y\)的路径上的点里面,必须是\(x\)最先被随机到。因此这种概率为\(\frac 1{dis(x,y)}\)

所以统计一下原树上每种长度的路径的条数,这个东西就是点分治+FFT的套路了。

#include<bits/stdc++.h>
using namespace std;
#define fec(i,x,y) (int i=head[x],y=g[i].to;i;i=g[i].ne,y=g[i].to)
#define dbg(...) fprintf(stderr,__VA_ARGS__)
#define File(x) freopen(#x".in","r",stdin),freopen(#x".out","w",stdout)
#define isin(x,S) (((S)>>((x)-1))&1)
#define fi first
#define se second
#define pb push_back
template<typename I>inline void read(I&x){int f=0,c;while(!isdigit(c=getchar()))c=='-'?f=1:0;x=c&15;while(isdigit(c=getchar()))x=(x<<1)+(x<<3)+(c&15);f?x=-x:0;}
template<typename A,typename B>inline char SMAX(A&a,const B&b){return a<b?a=b,1:0;}
template<typename A,typename B>inline char SMIN(A&a,const B&b){return a>b?a=b,1:0;}
typedef long long ll;typedef unsigned long long ull;typedef pair<int,int>pii;

const int N=30000+7;const double PI=acos(-1);
int n,x,y,ans[N];double ans1;
struct Edge{int to,ne;}g[N<<1];int head[N],tot;
inline void Addedge(int x,int y){g[++tot].to=y;g[tot].ne=head[x];head[x]=tot;}
inline void Adde(int x,int y){Addedge(x,y),Addedge(y,x);}

struct Complex{
	double x,y;
	inline Complex(const double&x=0,const double&y=0):x(x),y(y){}
	inline Complex operator+(const Complex&a){return Complex(x+a.x,y+a.y);}
	inline Complex operator-(const Complex&a){return Complex(x-a.x,y-a.y);}
	inline Complex operator*(const Complex&a){return Complex(x*a.x-y*a.y,x*a.y+y*a.x);}
}A[N<<2],B[N<<2];

inline void FFT(Complex*a,int n,int f){
	for(int i=0,j=0;i<n;++i){
		if(i>j)swap(a[i],a[j]);
		for(int l=(n>>1);(j^=l)<l;l>>=1);
	}
	for(int i=1;i<n;i<<=1){
		Complex w(cos(PI/i),f*sin(PI/i));
		for(int j=0;j<n;j+=(i<<1)){
			Complex e(1,0);
			for(int k=0;k<i;++k,e=e*w){
				Complex x=a[j+k],y=e*a[i+j+k];
				a[j+k]=x+y,a[i+j+k]=x-y;
			}
		}
	}
	if(f>0)return;
	for(int i=0;i<n;++i)a[i].x/=n;
}
inline void Mul(int*a,int*b,int n,int m){
	int l=1;while(l<=n+m)l<<=1;
	for(int i=0;i<=n;++i)A[i]=Complex(a[i],0);
	for(int i=0;i<=m;++i)B[i]=Complex(b[i],0);
	FFT(A,l,1);FFT(B,l,1);
	for(int i=0;i<l;++i)A[i]=A[i]*B[i];
	FFT(A,l,-1);
	for(int i=0;i<=min(n+m,::n);++i)ans[i]+=(int)(A[i].x+0.5);
	for(int i=0;i<l;++i)A[i]=B[i]=Complex();
}

int rt,mima,sum,num[N],vis[N];
inline void Grt(int x,int fa=0){
	num[x]=1;int f=0;
	for fec(i,x,y)if(y!=fa&&!vis[y])Grt(y,x),num[x]+=num[y],SMAX(f,num[y]);
	SMAX(f,sum-num[x]);if(SMIN(mima,f))rt=x;
}

int f[N],s[N],fa[N],q[N],dis[N],hd,tl;
inline void BFS(int x,int f){
	q[hd=0,tl=1]=x;fa[x]=f;
	while(hd<tl){
		int x=q[++hd];dis[x]=dis[fa[x]]+1;++::f[dis[x]-1];
		for fec(i,x,y)if(y!=fa[x]&&!vis[y])q[++tl]=y,fa[y]=x;
	}
}

inline void Calc(int x){
	s[1]=1;dis[x]=1;int mxd=1;ans[s[1]]++;
	for fec(i,x,y)if(!vis[y]){
		f[0]=0;BFS(y,x);num[y]=tl;
		Mul(s,f,mxd,dis[q[tl]]-1);
		SMAX(mxd,dis[q[tl]]);
		for(int i=1;i<=dis[q[tl]];++i)s[i]+=f[i-1],f[i-1]=0;
	}
	for(int i=1;i<=mxd;++i)s[i]=0;
}

inline void Solve(int x){
	vis[x]=1;Calc(x);
	for fec(i,x,y)if(!vis[y])mima=sum=num[y],Grt(y),Solve(rt);
}

int main(){
	#ifdef hzhkk
	freopen("hkk.in","r",stdin);
	#endif
	read(n);
	for(int i=1;i<n;++i)read(x),read(y),Adde(x+1,y+1);
	sum=mima=n;Grt(1);Solve(rt);
	ans1=ans[1];
	for(int i=2;i<=n;++i)ans1+=(double)ans[i]*2/i;
	printf("%.4lf\n",ans1);
}
posted @ 2019-01-11 20:59  hankeke303  阅读(123)  评论(0编辑  收藏  举报