BZOJ3451 Normal 点分治+FFT

题意：陈老师在点分治时随机选择重心然后分治，每次代价为树的大小，求期望代价

n<=30000

Sol:

这只是Normal 啊！！Lunatic要难成什么样啊TAT

由期望的线性性质，我们可以算出每个点的期望代价求和即为答案。

现在选定一个点，考虑其他点对她的贡献，当这个点是她到选定点路径上第一个被钦定的点时会产生1的贡献，所以贡献为1/dis(x,y) dis(i,j)表示i到j路径上的点数(包括端点)

然后可以通过点分治求出从一个重心向下的所有边的长度，看成多项式自己卷积自己一下再减去重复的后，对于现在多项式的第i项，对答案的贡献是a[i]/(i+1)

注意点分治的写法，用"总体减去各个子树"式肯定没问题，但如果要用"枚举每棵子树法"可能需要对深度排序来防止复杂度退化，太麻烦了我还是写第一种吧

Code:

#include<bits/stdc++.h>
#define debug(x) cout<<#x<<"="<<x<<endl
using namespace std;
const int N = 30009 , M = (1<<16)+9;
const double pi = acos(-1);

int first[N];
struct edg
{
	int next;
	int to;
}e[N<<1];
int e_sum;
int siz[N],W[N],rev[M],all[N],sub[N],b[M];
int n,sum,rt,len,bit,tot;
bool vis[N];
struct Complex
{
	double x,y;
	Complex(double _x=0.0,double _y=0.0):x(_x),y(_y){}
	Complex operator + (const Complex &X){return Complex(x+X.x,y+X.y);}
	Complex operator - (const Complex &X){return Complex(x-X.x,y-X.y);}
	Complex operator * (const Complex &X){return Complex(x*X.x-y*X.y,x*X.y+y*X.x);}
};
Complex a[M];
double ans;

inline int read()
{
	int x=0,f=1;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
	return x*f;
}
inline void add_edg(int x,int y)
{
	e_sum++;
	e[e_sum].next=first[x];
	first[x]=e_sum;
	e[e_sum].to=y;
}

void FFT(Complex *x,int type)
{
	for(int i=0;i<tot;i++) if(i<rev[i]) swap(x[i],x[rev[i]]);
	for(int i=1;i<tot;i<<=1)
	{
		Complex wn(cos(pi/i),1.0*type*sin(pi/i));
		for(int p=(i<<1),j=0;j<tot;j+=p)
		{
			Complex tmp(1,0);
			for(int k=0;k<i;k++,tmp=tmp*wn)
			{
				Complex L=x[j+k],R=tmp*x[j+k+i];
				x[j+k]=L+R;x[j+k+i]=L-R;
			}
		}
	}
}
void calc(int *arr,int opt)
{
	if(arr[0]==0) return ;
	len=0;
	for(int i=1;i<=arr[0];i++)
	{
		a[arr[i]].x+=1.0;
		len=max(len,arr[i]);
	}
	for(len<<=1,bit=0,tot=1;tot<=len;tot<<=1) bit++;
	for(int i=0;i<tot;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
	FFT(a,1);
	for(int i=0;i<tot;i++) a[i]=a[i]*a[i];
	FFT(a,-1);
	for(int i=1;i<=tot;i++) b[i]=a[i-1].x/tot+0.5;
	for(int i=1;i<=arr[0];i++) b[arr[i]*2+1]--;
	for(int i=1;i<=tot;i++) ans+=1.0*opt*b[i]/i;
	for(int i=0;i<tot;i++) a[i]=Complex(0,0);
}

void dfs(int x,int fa,int dis)
{
	ans+=2.0/(dis+1);
	all[++all[0]]=dis;sub[++sub[0]]=dis;
	for(int i=first[x];i;i=e[i].next)
	{
		int w=e[i].to;
		if(w==fa||vis[w]) continue;
		dfs(w,x,dis+1);
	}
}
void find(int x,int fa)
{
	siz[x]=1;W[x]=0;
	for(int i=first[x];i;i=e[i].next)
	{
		int w=e[i].to;
		if(w==fa||vis[w]) continue;
		find(w,x);
		siz[x]+=siz[w];
		W[x]=max(W[x],siz[w]);
	}W[x]=max(W[x],sum-siz[x]);
	if(W[x]<W[rt]) rt=x;
}
void solve(int x)
{
	vis[x]=1;all[0]=0;ans+=1.0;
	for(int i=first[x];i;i=e[i].next)
	{
		int w=e[i].to;
		if(vis[w]) continue;
		sub[0]=0;dfs(w,x,1);
		calc(sub,-1);
	}
	calc(all,1);
	for(int i=first[x];i;i=e[i].next)
	{
		int w=e[i].to;
		if(vis[w]) continue;
		sum=siz[w];rt=0;
		find(w,0);solve(rt);
	}
}

int main()
{
	n=read();
	for(int i=1;i<n;i++)
	{
		int x=read()+1,y=read()+1;
		add_edg(x,y);add_edg(y,x);
	}
	rt=0;W[0]=n+1;sum=n;
	find(1,0);solve(rt);
	printf("%.4lf\n",ans);
	return 0;
}