[期望+点分治+FFT]bzoj3451: Tyvj1953 Normal

题目描述

某天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里面的那个为标准)

输入

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

输出

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

样例输入

3
0 1
1 2

样例输出

5.6667

提示

n<=30000

题解

期望+FFT+点分治;

根据期望的线性性,Ans=\sum_{i=1}^n\sum_{j=1}^nE(ji的点分树子树中)

ji的点分树子树中,即i\to j这条路径上,i第一个被选;因为i\to j上每个点被选择的概率是一样的,所以P(ji的点分树子树中)=\large{\frac{1}{dis(i,j)}},这里的dis为两点之间的点数,又E(ji的点分树子树中)=\large{\frac{1}{P}},所以我们只需要求\sum_{i=1}^n\sum_{j=1}^ndis(i,j)即可;

对于每个点i,令F(x)dis=x的路径条数,选择任意两颗子树,dis(x)=\sum_{i=1}^xdis(i)[son1]*dis(x-i)[son2],这是个标准的卷积形式,用FFT快速求解;外层再套一个点分即可;

时间复杂度:FFTlogN次,长度逐渐增加,合起来为O(logN),即为O(NloglogN),点分治O(NlogN),最后总和为O(NloglogN+NlogN)

#include<cmath>
#include<cstdio>
#include<algorithm>
using namespace std;
typedef long long LL;
const double pi=acos(-1.0);
struct Cp{
    double _r,_i;
    Cp(double x=0,double y=0):_r(x),_i(y){}
    Cp friend operator + (Cp A,Cp B){return Cp(A._r+B._r,A._i+B._i);}
    Cp friend operator - (Cp A,Cp B){return Cp(A._r-B._r,A._i-B._i);}
    Cp friend operator * (Cp A,Cp B){return Cp(A._r*B._r-A._i*B._i,A._r*B._i+A._i*B._r);}
}A[262145];
int Rev[262145];
inline void getRev(int len){for (int i=1;i<(1<<len);i++) Rev[i]=((Rev[i>>1]>>1)|((i&1)<<(len-1)));}
void FFT(Cp *F,int len,int dft){
    for (int i=0;i<len;i++) if (i<Rev[i]) swap(F[i],F[Rev[i]]);
    for (int i=1;i<len;i<<=1){
        Cp Wn(cos(pi/i),dft*sin(pi/i));
        for (int j=0;j<len;j+=i<<1){
            Cp W(1,0);
            for (int k=j;k<i+j;k++){
                Cp x=F[k],y=W*F[i+k];
                F[k]=x+y;F[i+k]=x-y;W=W*Wn;
            }
        }
    }
    if (dft==-1) for (int i=0;i<len;i++) F[i]._r/=(double)len;
}
int N;
int tot,lnk[50005],son[100005],nxt[100005];
LL Ans[100005];
long double Res;
inline void add(int x,int y){son[++tot]=y;nxt[tot]=lnk[x];lnk[x]=tot;}
int Rt,Sum,Size[50005],F[50005],dep[50005],val[50005];
bool vis[50005];
void calc(int Miu){
    int Mx=0,len=0,S=1;
    for (int i=1;i<=tot;i++) Mx=max(Mx,val[i]);
    while (S<=Mx<<1) S<<=1,len++;getRev(len);
    for (int i=0;i<S;i++) A[i]=(0,0);
    for (int i=1;i<=tot;i++) A[val[i]]._r++;
    FFT(A,S,1);
    for (int i=0;i<S;i++) A[i]=A[i]*A[i];
    FFT(A,S,-1);
    for (int i=0;i<=Mx<<1;i++) Ans[i]+=Miu*(int)(A[i]._r+0.5);
}
void getRt(int x,int fa){
    Size[x]=1;F[x]=0;
    for (int i=lnk[x];i;i=nxt[i]){
        if (vis[son[i]]||son[i]==fa) continue;
        getRt(son[i],x);
        Size[x]+=Size[son[i]];
        F[x]=max(F[x],Size[son[i]]);
    }
    F[x]=max(F[x],Sum-Size[x]);
    if (F[x]<F[Rt]) Rt=x;
}
void dfs(int x,int fa){
    val[++tot]=dep[x];Size[x]=1;
    for (int i=lnk[x];i;i=nxt[i]) if (!vis[son[i]]&&son[i]!=fa) dep[son[i]]=dep[x]+1,dfs(son[i],x),Size[x]+=Size[son[i]];
}
void Solve(int x){
    dep[x]=tot=0;dfs(x,0);calc(1);vis[x]=1;
    for (int i=lnk[x];i;i=nxt[i]){
        if (vis[son[i]]) continue;
        Rt=0;Sum=Size[son[i]];
        dep[son[i]]=1;tot=0;dfs(son[i],0);calc(-1);
        Sum=Size[son[i]];Rt=0;getRt(son[i],0);Solve(Rt);
    }
}
inline int read(){
    int Res=0;char ch=getchar();
    while (ch>'9'||ch<'0') ch=getchar();
    while (ch>='0'&&ch<='9') Res=Res*10+ch-'0',ch=getchar();
    return Res;
}
int main()
{
    N=read();
    for (int i=1;i<N;i++){int x=read(),y=read();add(x+1,y+1);add(y+1,x+1);}
    dep[0]=-1;F[0]=2e9;Sum=N;getRt(1,0);Solve(Rt);
    for (int i=0;i<N;i++) Res+=(long double)Ans[i]/(long double)(i+1);
    printf("%.4Lf\n",Res);
    return 0;
}
0 Comments
Leave a Reply