bzoj3451/Tyvj1953:Normal(点分治+FFT)

来源:互联网 发布:淘宝卖家怎么催快递 编辑:程序博客网 时间:2024/05/17 22:13

题面
题意:给你一棵树,XJB选点来做点分治,问期望的复杂度。

根据我对期望的粗鄙理解,期望就是个积分,所以它满足积分加减等线性运算。故这里可以考虑对每个点求期望复杂度。

对于点x,它产生的期望就是x在点分树中的期望深度。

再次运用期望的线性性质,我们可以对于每个y,求出x在y子树中的概率,由此推出期望。

经过小小的分析,发现x在y的子树中,仅当在x-y的路径上,y被第一个选中来分治,则概率为1dis(x,y),也就是y能为x带来的期望。

然后原问题的答案为

i=1nj=1n1dis(i,j)

若dis(i,j)不在分母里,就是一个果的点分治了。

我们换下枚举的顺序,设f[x]为距离为x的点对数,变成nx=1f[x],就可以算答案了。

同样考虑点分治,找到重心x,考虑过x的路径。
对于每个连通块。求出它们每个点到x的距离,设g[i]为到x距离为i的点数。我们发现两个g能对答案的贡献是个卷积的形式,用fft优化即可。

期望复杂度O(nlog2n),但理论上说,不加优化的话,是会被菊花图卡的。
比如说这样
这里写图片描述

用一条长链把fft的次数界撑大了,然后用大大的次数界不断地做fft。
应该先做深度小的连通块,再做深度大的。

#include <iostream>#include <fstream>#include <algorithm>#include <cmath>#include <ctime>#include <cstdio>#include <cstdlib>#include <cstring>using namespace std;#define mmst(a, b) memset(a, b, sizeof(a))#define mmcp(a, b) memcpy(a, b, sizeof(b))typedef long long LL;const int N=60060,oo=1e9+7;const LL g=3,p=1004535809;int n,rev[N];int nn;int to[N],nex[N],head[N],cnt;int root,sum,siz[N],d[N],dep[N],len,slen;LL a[N],b[N],c[N],ans[N];double res=0.0;bool vis[N];void read(int &hy){    hy=0;    char cc=getchar();    while(cc<'0'||cc>'9')    cc=getchar();    while(cc>='0'&&cc<='9')    {        hy=(hy<<3)+(hy<<1)+cc-'0';        cc=getchar();    }}LL cheng(LL a,LL b){    LL res=1ll;    for(;b;b>>=1,a=a*a%p)    if(b&1)    res=res*a%p;    return res;}void init(int lim){    int k=-1;    n=1;    while(n<lim)    n<<=1,k++;    for(int i=0;i<n;i++)    rev[i]=(rev[i>>1] >> 1) | ((i&1)<<k);}void ntt(LL *a,int ops){    for(int i=0;i<n;i++)    if(i<rev[i])    swap(a[i],a[rev[i]]);    for(int l=2;l<=n;l<<=1)    {        int m=l>>1;        LL wn;        if(ops)        wn=cheng(g,(p-1)/l);        else        wn=cheng(g,p-1-(p-1)/l);        for(int i=0;i<n;i+=l)        {            LL w=1ll;            for(int k=0;k<m;k++)            {                LL t=a[i+k+m]*w%p;                a[i+k+m]=(a[i+k]-t+p)%p;                a[i+k]=(a[i+k]+t)%p;                w=w*wn%p;            }        }    }    if(!ops)    {        LL Inv=cheng(n,p-2);        for(int i=0;i<n;i++)        a[i]=a[i]*Inv%p;    }}void add(int u,int v){    to[++cnt]=v;    nex[cnt]=head[u];    head[u]=cnt;}void dfsRoot(int x,int fa){    d[x]=0;    siz[x]=1;    for(int h=head[x];h;h=nex[h])    if(!vis[to[h]]&&to[h]!=fa)    {        dfsRoot(to[h],x);        siz[x]+=siz[to[h]];        d[x]=max(d[x],siz[to[h]]);    }    d[x]=max(d[x],sum-siz[x]);    if(d[x]<d[root])    root=x;}void dfsLen(int x,int fa){    len=max(len,dep[x]);    a[dep[x]]++;    for(int h=head[x];h;h=nex[h])    if(!vis[to[h]]&&to[h]!=fa)    {        dep[to[h]]=dep[x]+1;            dfsLen(to[h],x);    }}void dfsSol(int x){    vis[x]=1;    slen=1;    b[0]=1;    for(int h=head[x];h;h=nex[h])    if(!vis[to[h]])    {        len=0;        dep[to[h]]=1;        dfsLen(to[h],0);        init(len+slen+1);        for(int i=0;i<n;i++)        c[i]=b[i];        for(int i=0;i<=len;i++)        b[i]+=a[i];        ntt(c,1);        ntt(a,1);        for(int i=0;i<n;i++)        c[i]=c[i]*a[i]%p;        ntt(c,0);        for(int i=0;i<n;i++)        ans[i]+=c[i];        for(int i=0;i<n;i++)        a[i]=0;        slen=max(slen,len);    }    for(int i=0;i<=slen;i++)    b[i]=0;    for(int h=head[x];h;h=nex[h])    if(!vis[to[h]])    {        root=0;        sum=siz[to[h]];        dfsRoot(to[h],0);        dfsSol(root);    }}int main(){    cin>>nn;    for(int i=1;i<nn;i++)    {        int u,v;        read(u);        read(v);        u++;        v++;        add(u,v);        add(v,u);    }    sum=nn;    root=0;    d[0]=oo;    dfsRoot(1,0);    dfsSol(root);    for(int i=1;i<nn;i++)    res+=2.0*ans[i]/(i+1);    res+=nn;    printf("%.4lf\n",res);}

这里写图片描述

原创粉丝点击