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

Description

某天WJMZBMR学习了一个神奇的算法：树的点分治！
这个算法的核心是这样的：
消耗时间=0
Solve(树 a)
 消耗时间 += a 的 大小
 如果 a 中 只有 1 个点
  退出
 否则在a中选一个点x，在a中删除点x

 那么a变成了几个小一点的树，对每个小树递归调用Solve
我们注意到的这个算法的时间复杂度跟选择的点x是密切相关的。
如果x是树的重心，那么时间复杂度就是O(nlogn)
但是由于WJMZBMR比较傻逼，他决定随机在a中选择一个点作为x！
Sevenkplus告诉他这样做的最坏复杂度是O(n^2)

但是WJMZBMR就是不信>_<。。。
于是Sevenkplus花了几分钟写了一个程序证明了这一点。。。你也试试看吧^_^
现在给你一颗树，你能告诉WJMZBMR他的傻逼算法需要的期望消耗时间吗？（消耗时间按在Solve里面的那个为标准）

Input

第一行一个整数n，表示树的大小
接下来n-1行每行两个数a,b，表示a和b之间有一条边
注意点是从0开始标号的

Output

一行一个浮点数表示答案
四舍五入到小数点后4位
如果害怕精度跪建议用long double或者extended

Sample Input

3
0 1
1 2

Sample Output

5.6667

HINT

n<=30000

分析：
期望好题

对于一棵树的期望，实际上等于ta所有子树的期望之和

考虑对于点$u$，$v$在其点分子树中的概率是多少，
$v$在$u$的点分子树中，就是说$u$到$v$的路径中，$u$是第一个被选的
那么期望即为$1 \over dis(u,v)$
那么题目所求$\sum{1 \over dis(u,v)}$

点分治
对于一棵树内的通过重心的所有链，两两都可以配对
所以我们可以用FFT解决这个类似组合数学的问题
设计生成函数$A$（如果有长度为$l$的链，那么$x_l$的系数++），$A^2$就是答案了

注意

点分有不同的写法，不过这道题要采用：
总的 - 同一子树中路径 的方法
不能挨个枚举子树与子树间转移，因为FFT的复杂度与最大深度有关，相同深度转移多次，菊花图就可以卡掉

tip

分析中所说的dis指的是经过的点数（不是边数）
所以统计时：

for (int i=1;i<fn;i++) {
    t=(ll)(A[i].x/fn+0.5);      //dis的次数
    sum+=(LD)t/(LD)(i+1);       
    }

而单点也会给答案贡献1

#include<cstdio>
#include<cstring>
#include<iostream>
#include<cmath>
#define LD double
#define ll long long 

using namespace std;

const int N=80005;
const double Pi=acos(-1.0);

struct node{
    double x,y;
    node(double xx=0,double yy=0) {
        x=xx;y=yy;
    }
};
node A[N];
node operator +(const node &A,const node &B) {return node(A.x+B.x,A.y+B.y);}
node operator -(const node &A,const node &B) {return node(A.x-B.x,A.y-B.y);}
node operator *(const node &A,const node &B) {return node(A.x*B.x-A.y*B.y,A.x*B.y+A.y*B.x);}

struct E{
    int y,nxt;
};
E way[N<<1];
int n,st[N],tot=0,sz,size[N],root,f[N];
bool vis[N];
LD ans;

void add(int u,int w) {
    tot++;way[tot].y=w;way[tot].nxt=st[u];st[u]=tot;
    tot++;way[tot].y=u;way[tot].nxt=st[w];st[w]=tot;
}

void findroot(int now,int fa) {
    size[now]=1;
    f[now]=0;
    for (int i=st[now];i;i=way[i].nxt)
        if (way[i].y!=fa&&!vis[way[i].y]) {
            findroot(way[i].y,now);
            size[now]+=size[way[i].y];
            f[now]=max(f[now],size[way[i].y]);
        }
    f[now]=max(f[now],sz-size[now]);
    if (f[now]<f[root]) root=now;
}

void FFT(int n,node *a,int opt) {
    int i,j=0,k;
    for (i=0;i<n;i++) {
        if (i>j) swap(a[i],a[j]);
        for (int l=n>>1;(j^=l)<l;l>>=1);
    }
    for (i=2;i<=n;i<<=1) {
        int m=i>>1;
        node wn(cos(2.0*opt*Pi/i),sin(2.0*opt*Pi/i));
        for (j=0;j<n;j+=i) {
            node w(1,0);
            for (k=0;k<m;k++,w=w*wn) {
                node z=a[j+k+m]*w;
                a[j+k+m]=a[j+k]-z;
                a[j+k]=a[j+k]+z;
            }
        }
    }
}

int d[N],deep[N],cnt,maxdeep;

void dfs(int now,int fa) {
    d[++cnt]=deep[now];
    maxdeep=max(maxdeep,deep[now]);
    for (int i=st[now];i;i=way[i].nxt)
        if (way[i].y!=fa&&!vis[way[i].y]) {
            deep[way[i].y]=deep[now]+1;
            dfs(way[i].y,now);
        }
}

void cal(int now,int dep,int opt) {
    deep[now]=dep;
    cnt=0; maxdeep=0;
    dfs(now,0);

    int fn=1;
    while (fn<=maxdeep*2) fn<<=1;
    for (int i=0;i<=fn;i++) A[i]=node(0.0,0.0);
    for (int i=1;i<=cnt;i++) A[d[i]].x++;
    FFT(fn,A,1);
    for (int i=0;i<fn;i++) A[i]=A[i]*A[i];
    FFT(fn,A,-1);
    LD sum=0; ll t;
    for (int i=1;i<fn;i++) {
        t=(ll)(A[i].x/fn+0.5);      //dis的次数
        sum+=(LD)t/(LD)(i+1);       // /(i+1)
    }
    ans=ans+1.0*opt*sum;
}

void solve(int now) {
    ans+=1.0;                       //单点的贡献 
    cal(now,0,1);
    vis[now]=1;
    for (int i=st[now];i;i=way[i].nxt)
        if (!vis[way[i].y]) {
            cal(way[i].y,1,-1);
            sz=size[way[i].y]; root=0;
            findroot(way[i].y,0);
            solve(root);
        }
}

int main()
{
    scanf("%d",&n); ans=0;
    for (int i=1;i<n;i++) {
        int u,w;
        scanf("%d%d",&u,&w);
        u++; w++;
        add(u,w);
    }
    sz=n; root=0; f[0]=N;
    findroot(1,0);
    solve(root);
    printf("%0.4lf\n",(double)ans);
    return 0;
}