bzoj 3451: 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

分析：
答案其实就是所有节点的期望深度和。
对于一个点$x$，考虑一个点$y$对他的贡献。当前仅当在$x$到$y$路径上的点中直接选择$y$的概率，因为贡献是1，所以就是$\frac{1}{dis(x,y)}$。
所以答案就是${\sum }_{x=1}^n{\sum }_{y=1}^n\frac{1}{dis(x,y)}$
其中$dis(x,y)$是指$x$到$y$路径中点的数目。
我们考虑进行点分治，对于一个分治中心，跑出各种路径的长度，合并可以使用fft。
对于分治中心，我们不能使用子树与前面子树合并的方法，因为fft的复杂度取决于结果的多项式的长度，复杂化度可能达到$O(n^{2}logn)$级。
所以我们可以直接对所有链进行匹配，然后减去两条链都在同一棵子树的答案。
复杂度是$O(nlo{g}^{2}n)$的。

代码：

/**************************************************************
    Problem: 3451
    User: ypxrain
    Language: C++
    Result: Accepted
    Time:976 ms
    Memory:10820 kb
****************************************************************/
 
#include <iostream>
#include <cstdio>
#include <cmath>
#include <queue>
 
const int maxn=1e5+7;
const double pi=acos(-1);
 
using namespace std;
 
int n,x,y,cnt,len,sum,root;
int ls[maxn],size[maxn],height[maxn],dis[maxn],vis[maxn];
int r[maxn],f[maxn],g[maxn],h[maxn];
 
double ans;
 
struct edge{
    int y,next;
}e[maxn];
 
queue <int> q;
 
struct rec{
    double x,y;
}dft[maxn],idft[maxn],w[maxn];
 
rec operator +(rec a,rec b)
{
    return (rec){a.x+b.x,a.y+b.y};
}
 
rec operator -(rec a,rec b)
{
    return (rec){a.x-b.x,a.y-b.y};
}
 
rec operator *(rec a,rec b)
{
    return (rec){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};
}
 
rec operator !(rec a)
{
    return (rec){a.x,-a.y};
}
 
void fft(rec *a,int f)
{
    for (int i=0;i<len;i++)
    {
        if (i<r[i]) swap(a[i],a[r[i]]);
    }
    w[0]=(rec){1,0};
    for (int i=2;i<=len;i*=2)
    {
        rec wn=(rec){cos(2*pi/i),f*sin(2*pi/i)};
        for (int j=i/2;j>=0;j-=2) w[j]=w[j/2];
        for (int j=1;j<i/2;j+=2) w[j]=w[j-1]*wn;
        for (int j=0;j<len;j+=i)
        {
            for (int k=0;k<i/2;k++)
            {
                rec u=a[j+k],v=a[j+k+i/2]*w[k];
                a[j+k]=u+v;
                a[j+k+i/2]=u-v;
            }
        }
    }
}
 
void FFT(int *a,int *b,int *c,int n,int m)
{
    len=1;
    int k=0;
    while (len<(n+m)) len*=2,k++;
    for (int i=0;i<len;i++)
    {
        r[i]=(r[i>>1]>>1)|((i&1)<<(k-1));
    }
    for (int i=0;i<len;i++)
    {
        if (i<n) dft[i].x=a[i]; else dft[i].x=0;
        if (i<m) dft[i].y=b[i]; else dft[i].y=0;
    }
    fft(dft,1); 
    for (int i=0;i<len;i++)
    {
        int j=(len-1)&(len-i);
        rec da,db;
        da=(dft[i]+(!dft[j]))*(rec){0.5,0};
        db=(dft[i]-(!dft[j]))*(rec){0,-0.5};
        idft[i]=da*db;
    }
    fft(idft,-1);
    for (int i=0;i<len;i++) c[i]=trunc(idft[i].x/len+0.5);
}
 
void add(int x,int y)
{
    e[++cnt]=(edge){y,ls[x]};
    ls[x]=cnt;
}
 
void findroot(int x,int fa)
{
    size[x]=1;
    height[x]=0;
    for (int i=ls[x];i>0;i=e[i].next)
    {
        int y=e[i].y;
        if ((y==fa) || (vis[y])) continue;
        findroot(y,x);
        size[x]+=size[y];
        height[x]=max(height[x],size[y]);
    }
    height[x]=max(height[x],sum-size[x]);
    if ((!root) || (height[x]<height[root])) root=x;
}
 
void dfs(int x,int fa)
{
    size[x]=1;
    for (int i=ls[x];i>0;i=e[i].next)
    {
        int y=e[i].y;
        if ((y==fa) || (vis[y])) continue;
        dis[y]=dis[x]+1;
        g[dis[y]]++;
        f[dis[y]]++;
        dfs(y,x);
        size[x]+=size[y];
        if (x==root)
        {                   
            FFT(g,g,h,size[y]+1,size[y]+1);
            for (int i=0;i<(size[y]+1)*2;i++)
            {
                ans-=(double)h[i]/(double)(i+1);
                g[i]=h[i]=0;
            }
            q.push(y);
        }
    }
}
 
int main()
{
    scanf("%d",&n);
    for (int i=1;i<n;i++)
    {
        scanf("%d%d",&x,&y);
        add(x+1,y+1);
        add(y+1,x+1);
    }       
    q.push(1);
    size[1]=n;  
    while (!q.empty())
    {
        int x=q.front();
        q.pop();
        sum=size[x];
        root=0;
        findroot(x,0);
        dis[root]=0;
        for (int i=0;i<sum;i++) f[i]=g[i]=0;
        f[0]=1;
        dfs(root,0);
        vis[root]=1;
        FFT(f,f,h,sum,sum);
        for (int i=0;i<sum*2;i++) ans+=(double)h[i]/(double)(i+1);
    }
    printf("%.4lf",ans);
}