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是第一个被选的，
于是概率就是$\frac{1}{dis(u,v)}$。
那么原题即求$\sum\frac{1}{dis(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[N],to[2*N],nxt[2*N],num=0,size[N],ss[N],sz=0,root,cnt[MXN],mx,R[MXN],len,p,dis[MXN];
bool del[N];
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;
}