【BZOJ3451】Normal

【BZOJ3451】Normal

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

题目大意就是给定一棵树，问随机进行点分治（不一定找重心）每个节点期望访问次数之和。

神题啊。

我们考虑计算点对\(P_{x,y}\)表示在\(x\)最为点分重心的时候访问了\(y\)的概率，由期望的线性性得出\(ans=\sum_{x=1}^n\sum_{y=1}^n P_{x,y}\)。

然后我们考虑怎么求这个\(P_{x,y}\)。假设\(x\)到\(y\)最短路上有\(ver_{x,y}\)个点，那么概率就是\(\frac{1}{ver_{x,y}}\)。

为什么呢？我们可以考虑将点分看成一个删点的操作。我们必须保证在删除\(x\)之前\(x\)到\(y\)的路径都是未被删除的。我们假设可以重复删除以删除的点，则：
\[ \displaystyle \begin{align} P_{x,y}&=\sum_{i=0}^{\infty}(\frac{n-ver_{x,y}}{n})^i\frac{1}{n}\\ & =\frac{1}{1-\frac{n-ver_{x,y}}{n}}\frac{1}{n}\\ &=\frac{1}{ver_{x,y}} \end{align} \]
这个公式的意义就是考虑在删除\(x\)之前删除了多少次其他点。

所以我们要求的就是\(ans=ans=\sum_{x=1}^n\sum_{y=1}^n \frac{1}{ver_{x,y}}\)

因为这个公式是非线性的，所以我们对于每个\(k\)要求出\([ver_{x,y}==k]\)的数量。这个我们就可以考虑点分治了。

#include<bits/stdc++.h>
#define ll long long
#define N 30005

using namespace std;
inline int Get() {int x=0,f=1;char ch=getchar();while(ch<'0'||ch>'9') {if(ch=='-') f=-1;ch=getchar();}while('0'<=ch&&ch<='9') {x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}return x*f;}

int n;
struct Com {
    long double r,v;
    Com() {r=v=0;}
    Com(double a,double b) {r=a,v=b;}
};
Com operator +(const Com &a,const Com &b) {return Com(a.r+b.r,a.v+b.v);}
Com operator -(const Com &a,const Com &b) {return Com(a.r-b.r,a.v-b.v);}
Com operator *(const Com &a,const Com &b) {return Com(a.r*b.r-a.v*b.v,a.r*b.v+a.v*b.r);}
Com operator /(const Com &a,const long double b) {return Com(a.r/b,a.v/b);}
const double pi=acos(-1);
void FFT(Com *a,int d,int flag) {
    int n=1<<d;
    static int rev[N<<2];
    for(int i=0;i<n;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<d-1);
    for(int i=0;i<n;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
    for(int s=1;s<=d;s++) {
        int len=1<<s,mid=len>>1;
        Com w(cos(2*flag*pi/len),sin(2*flag*pi/len));
        for(int i=0;i<n;i+=len) {
            Com t(1,0);
            for(int j=0;j<mid;j++,t=t*w) {
                Com u=a[i+j],v=a[i+j+mid]*t;
                a[i+j]=u+v;
                a[i+j+mid]=u-v;
            }
        }
    }
    if(flag==-1) for(int i=0;i<n;i++) a[i]=a[i]/n;
}

struct road {
    int to,nxt;
}s[N<<1];
int h[N],cnt;
void add(int i,int j) {s[++cnt]=(road) {j,h[i]};h[i]=cnt;}
int size[N],mx[N],sum,rt;
int fr[N];
bool vis[N];

void Get_root(int v,int fr) {
    mx[v]=size[v]=1;
    ::fr[v]=fr;
    for(int i=h[v];i;i=s[i].nxt) {
        int to=s[i].to;
        if(to==fr||vis[to]) continue ;
        Get_root(to,v);
        size[v]+=size[to];
        mx[v]=max(mx[v],size[to]);
    }
    mx[v]=max(mx[v],sum-size[v]);
    if(mx[rt]>mx[v]) rt=v;
}

int ans[N];
int tem[N];
int mxdep;

void statis(int v,int fr,int dep) {
    tem[dep]++;
    mxdep=max(mxdep,dep);
    for(int i=h[v];i;i=s[i].nxt) {
        int to=s[i].to;
        if(to==fr||vis[to]) continue ;
        statis(to,v,dep+1);
    }
}

Com A[N<<2];
void cal(int flag) {
    int d=ceil(log2(mxdep<<1|1));
    for(int i=0;i<1<<d;i++) A[i]=Com(tem[i],0);
    FFT(A,d,1);
    for(int i=0;i<1<<d;i++) A[i]=A[i]*A[i];
    FFT(A,d,-1);
    for(int i=0;i<1<<d;i++) {
        ans[i+1]+=flag*(int(A[i].r+0.5));
    }
}

void solve(int v) {
    vis[v]=1;
    if(fr[v]) size[fr[v]]=sum-size[v];
    mxdep=0;
    for(int i=h[v];i;i=s[i].nxt) {
        int to=s[i].to;
        if(vis[to]) continue ;
        statis(to,v,1);
    }
    ans[1]++;
    for(int i=1;i<=mxdep;i++) ans[i+1]+=tem[i]*2;
    cal(1);
    for(int i=1;i<=mxdep;i++) tem[i]=0;
    for(int i=h[v];i;i=s[i].nxt) {
        int to=s[i].to;
        if(vis[to]) continue ;
        mxdep=0;
        statis(to,v,1);
        cal(-1);
        for(int j=1;j<=mxdep;j++) tem[j]=0;
        sum=size[to];
        rt=0;
        Get_root(to,v);
        solve(rt);
    }
}

int main() {
    mx[0]=1e9;
    n=Get();
    int a,b;
    for(int i=1;i<n;i++) {
        a=Get()+1,b=Get()+1;
        add(a,b),add(b,a);
    }
    sum=n;
    Get_root(1,0);
    solve(rt);
    long double Ans=0;
    for(int i=1;i<=n;i++) {
        Ans+=1.0*ans[i]/(1.0*i);
    }
    cout<<fixed<<setprecision(4)<<Ans;
    return 0;
}

转载于:https://www.cnblogs.com/hchhch233/p/10585680.html