bzoj3451 Tyvj1953 Normal（期望概率+点分治+FFT）

bzoj3451 Tyvj1953 Normal

原题地址：http://www.lydsy.com/JudgeOnline/problem.php?id=3451

题意：

某天WJMZBMR学习了一个神奇的算法：树的点分治！

这个算法的核心是这样的：

消耗时间=0

Solve(树 a)

消耗时间 += a 的 大小

如果 a 中 只有 1 个点

退出

否则在a中选一个点x，在a中删除点x

那么a变成了几个小一点的树，对每个小树递归调用Solve

我们注意到的这个算法的时间复杂度跟选择的点x是密切相关的。

如果x是树的重心，那么时间复杂度就是O(nlogn)

给定一棵树，假如每次在a中随机选一个点作为x，请求出这个算法的期望复杂度。

数据范围

树的节点数n<=30000

题解：

好题，期望妙妙。

首先，根据期望的线性性，可以算每个点的期望贡献，就是它期望的点分树子树大小。

这题是要求和，那么依旧可以分开算贡献，

考虑对于点u，v在其点分子树中的概率是多少，

v在u的点分子树中，就是说u到v的路径中，u是第一个被选的，

于是概率就是1dis(u,v)。

那么原题即求∑1dis(u,v)

要求对于每个距离，有多少条路径，

那么就点分治来找出从分治重心出发的各种路径，合成一条用FFT来转移。

注意要采用 答案=总的 - 同一子树中的 的方法，不能挨个枚举子树子树间转移，因为FFT的复杂度与最大深度有关，相同深度转移多次，菊花图就可以卡。

代码：

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<cmath>
#define LD long double
using namespace std;
const long double Pi=acos(-1);
const int N=30005;
const int MXN=65536+1000;
struct Virt
{
LD r,i;
Virt(){}
Virt(LD r,LD i):r(r),i(i){}
Virt operator+(const Virt &A){return Virt(r+A.r,i+A.i);}
Virt operator-(const Virt &A){return Virt(r-A.r,i-A.i);}
Virt operator*(const Virt &A){return Virt(r*A.r-i*A.i,r*A.i+i*A.r);}
}a[MXN];
int n,head
,to[2*N],nxt[2*N],num=0,size
,ss
,sz=0,root,cnt[MXN],mx,R[MXN],len,p,dis[MXN];
bool del
;
void build(int u,int v)
{
num++;
to[num]=v;
nxt[num]=head[u];
head[u]=num;
}
void getroot(int u,int fa)
{
size[u]=1; ss[u]=0;
for(int i=head[u];i;i=nxt[i])
{
int v=to[i]; if(v==fa||del[v]) continue;
getroot(v,u);
size[u]+=size[v];
ss[u]=max(ss[u],size[v]);
}
ss[u]=max(ss[u],sz-size[u]);
if(ss[u]<ss[root]) root=u;
}
void getdis(int u,int d,int fa)
{
cnt[d]++; size[u]=1; mx=max(mx,d);
for(int i=head[u];i;i=nxt[i])
{
int v=to[i]; if(v==fa||del[v]) continue;
size[u]+=size[v];
getdis(v,d+1,u);
}
}
void FFT(Virt *x,int opt)
{
for(int i=0;i<len;i++) if(i<R[i]) swap(x[i],x[R[i]]);
for(int m=2;m<=len;m<<=1)
{
int l=m>>1; Virt wn=Virt(cos((LD)2*Pi/m),sin((LD)2*Pi/m)); if(opt==-1) wn.i=-wn.i;
for(int j=0;j<len;j+=m)
{
Virt wi=Virt(1,0);
for(int i=0;i<l;i++)
{
Virt y=wi*x[i+j+l];
x[i+j+l]=x[i+j]-y;
x[i+j]=x[i+j]+y;
wi=wi*wn;
}
}
}
if(opt==-1) for(int i=0;i<len;i++) x[i].r=(LD)x[i].r/(LD)len;
}
void cal(int u,int d,int opt)
{
mx=0; getdis(u,d,0);
for(len=1,p=0;len<(2*mx+2);len<<=1,p++);
R[0]=0; for(int i=1;i<len;i++) R[i]=(R[i>>1]>>1)|((i&1)<<(p-1));
for(int i=0;i<len;i++) a[i]=Virt((LD)cnt[i],0),cnt[i]=0;
FFT(a,1);
for(int i=0;i<len;i++) a[i]=a[i]*a[i]; FFT(a,-1);
for(int i=0;i<len;i++) dis[i]=dis[i]+opt*(int)(a[i].r+0.5);
}
void dfs(int x)
{

getroot(x,0); int u=root;
cal(u,0,1);
del[u]=1;
for(int i=head[u];i;i=nxt[i])
{
int v=to[i]; if(del[v]) continue;
cal(v,1,-1);
root=0; sz=size[v]; ss[0]=n+1;
dfs(v);
}
}
int main()
{
scanf("%d",&n);
for(int i=1;i<n;i++)
{
int u,v; scanf("%d%d",&u,&v); u++,v++;
build(u,v); build(v,u);
}
sz=n; root=0; ss[0]=n+1;
dfs(1);
double ans=0;
for(int i=1;i<=n;i++) if(dis[i-1]) ans+=(double)dis[i-1]/(double)i;
printf("%0.4lf\n",ans);
return 0;
}