【BZOJ3451】Normal （點分治）

題面

BZOJ

題解

顯然考慮每個點的貢獻。但是發現似乎怎麽算都不好計算其在點分樹上的深度。
那麽考慮一下這個點在點分樹中每一次被計算的情況，顯然就是其在某個點的點分樹內時才會被計算答案。
那麽設\(p[i][j]\)表示\(i\)在\(j\)的點分樹裏面的概率。
那麽答案就變成了\(\sum_i\sum_j p[i][j]\)
那麽\(i\)在\(j\)的點分樹的概率顯然就是兩點之間路徑不被斷開的概率，即\(\frac{1}{dis(i,j)+1}\)。
那麽點分治+\(FFT\)統計即可。

#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
#define MAX 70000
const double Pi=acos(-1);
inline int read()
{
    int x=0;bool t=false;char ch=getchar();
    while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
    if(ch=='-')t=true,ch=getchar();
    while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
    return t?-x:x;
}
struct Complex{double a,b;}A[MAX],B[MAX],W[MAX];
Complex operator+(Complex a,Complex b){return (Complex){a.a+b.a,a.b+b.b};}
Complex operator-(Complex a,Complex b){return (Complex){a.a-b.a,a.b-b.b};}
Complex operator*(Complex a,Complex b){return (Complex){a.a*b.a-a.b*b.b,a.b*b.a+a.a*b.b};}
int r[MAX];
void FFT(Complex *P,int N,int opt)
{
    for(int i=0;i<N;++i)if(i<r[i])swap(P[i],P[r[i]]);
    for(int i=1;i<N;i<<=1)
        for(int p=i<<1,j=0;j<N;j+=p)
            for(int k=0;k<i;++k)
            {
                Complex w=W[N/i*k];w.b*=opt;
                Complex X=P[j+k],Y=P[i+j+k]*w;
                P[j+k]=X+Y;P[i+j+k]=X-Y;
            }
    if(opt==-1)for(int i=0;i<N;++i)P[i].a/=N;
}
struct Line{int v,next;}e[MAX<<1];
int h[MAX],cnt=1;
inline void Add(int u,int v){e[cnt]=(Line){v,h[u]};h[u]=cnt++;}
int size[MAX],rt,mx,Size;
bool vis[MAX];
void Getroot(int u,int ff)
{
    size[u]=1;int ret=0;
    for(int i=h[u];i;i=e[i].next)
    {
        int v=e[i].v;if(v==ff||vis[v])continue;
        Getroot(v,u);size[u]+=size[v];
        ret=max(ret,size[v]);
    }
    ret=max(ret,Size-size[u]);
    if(mx>ret)mx=ret,rt=u;
}
int Ans[MAX];
int t[MAX],s[MAX],ln[MAX],md;
void dfs(int u,int ff,int dep)
{
    t[dep]+=1;md=max(md,dep);
    for(int i=h[u];i;i=e[i].next)
    {
        int v=e[i].v;if(v==ff||vis[v])continue;
        dfs(v,u,dep+1);
    }
}
void Multi(int *a,int *b,int n,int m,int *c)
{
    int N,l=0;for(N=1;N<=n+m;N<<=1)++l;
    for(int i=0;i<N;++i)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    for(int i=1;i<N;i<<=1)
        for(int k=0;k<i;++k)
            W[N/i*k]=(Complex){cos(Pi/i*k),sin(Pi/i*k)};
    for(int i=0;i<=n;++i)A[i].a=a[i];
    for(int i=0;i<=m;++i)B[i].a=b[i];
    FFT(A,N,1);FFT(B,N,1);
    for(int i=0;i<N;++i)A[i]=A[i]*B[i];
    FFT(A,N,-1);
    for(int i=0;i<=n+m;++i)c[i]=A[i].a+0.5;
    for(int i=0;i<N;++i)A[i]=B[i]=(Complex){0,0};
}
void Divide(int u)
{
    vis[u]=true;int mxd=0;
    s[0]=1;
    for(int i=h[u];i;i=e[i].next)
    {
        int v=e[i].v;if(vis[v])continue;
        md=0;dfs(e[i].v,u,1);
        Multi(s,t,mxd,md,ln);
        for(int i=0;i<=mxd+md;++i)Ans[i]+=ln[i],ln[i]=0;
        for(int i=0;i<=md;++i)s[i]+=t[i],t[i]=0;
        mxd=max(mxd,md);
    }
    for(int i=0;i<=mxd;++i)s[i]=0;
    for(int i=h[u];i;i=e[i].next)
    {
        int v=e[i].v;if(vis[v])continue;
        Size=mx=size[v];Getroot(v,u);
        Divide(rt);
    }
}
int n;double ans;
int main()
{
    n=read();
    for(int i=1;i<n;++i)
    {
        int u=read()+1,v=read()+1;
        Add(u,v);Add(v,u);
    }
    int N;for(N=1;N<=n;N<<=1);
    for(int i=1;i<N;i<<=1)
        for(int k=0;k<i;++k)
            W[N/i*k]=(Complex){cos(Pi/i*k),sin(Pi/i*k)};
    Size=mx=n;Getroot(1,0);Divide(rt);
    for(int i=1;i<=n;++i)ans+=1.0*Ans[i]/(i+1);
    ans*=2;ans+=n;
    printf("%.4lf\n",ans);
    return 0;
}
 

【BZOJ3451】Normal （點分治）