BZOJ3451: Tyvj1953 Normal【FFT+点分治】

3451: Tyvj1953 Normal

【题目描述】

传送门

【题解】

我们考虑x对y点的贡献,当x点对y点有贡献,说明在选择y之前没有选择在x->y的路径上的点,也就是y是x点分树的祖先,所以期望就是1/(dis(x,y)+1)
总的贡献就是\large \sum_{i=1}^{n} \sum_{j=1}^{n}\frac{1}{dis(i,j)+1}
求出距离为i的点对数为S[i],答案就是\large \sum_{i=1}^{n} \frac{2 * S[i]}{i+1}
用点分治处理,所以距离为i的点对数就是\sum_{k=0}^{i}dis(k) * dis(i-k)
点分治+FFT可以解决。

我们粗略的估计一下时间复杂度,假设每次FFT的长度都是n,那么复杂度就是O(n log^2(n)),但是因为是分治,所以每次FFT的长度是log(n),所以复杂度应该是O(nlog(n)+nlog(log(n)))

代码如下

#include<cmath>
#include<cstdio>
#include<complex>
#include<algorithm>
using namespace std;
typedef complex<double> CP;
typedef long long LL;
const int MAXN=50005;
const int MOD=998244353;
const double PI=acos(-1.0);
CP a[MAXN<<2];
int n,m,Al,Now,Root,MaxSiz,rev[MAXN],Siz[MAXN],Num[MAXN],Tim[MAXN];
LL Ans[MAXN<<1];
bool vis[MAXN];
struct Edge{
    int tot,lnk[MAXN],nxt[MAXN<<1],son[MAXN<<1];
    void Add(int x,int y){son[++tot]=y;nxt[tot]=lnk[x];lnk[x]=tot;}
}E;
void DFS(int x,int fa,int Dep){
    if(Tim[Dep]!=Now) Num[Dep]=0,Tim[Dep]=Now;
    Num[Dep]++;m=max(m,Dep);Siz[x]=1;
    for(int j=E.lnk[x];j;j=E.nxt[j])
    if(E.son[j]!=fa&&!vis[E.son[j]]) DFS(E.son[j],x,Dep+1),Siz[x]+=Siz[E.son[j]];
}
int Get_Root(int x,int fa){
    int NowMax=0;
    for(int j=E.lnk[x];j;j=E.nxt[j])
    if(E.son[j]!=fa&&!vis[E.son[j]]){
        Get_Root(E.son[j],x);
        NowMax=max(NowMax,Siz[E.son[j]]);
    }
    NowMax=max(Al-Siz[x],NowMax);
    if(NowMax<MaxSiz) Root=x,MaxSiz=NowMax;
}
void FFT(CP *a,int Len,int f){
    for(int i=0;i<Len;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
    for(int Step=1;Step<Len;Step<<=1){
        CP WN=exp(CP(0,f*PI/Step));
        for(int j=0;j<Len;j+=Step<<1){
            CP WNK(1,0);
            for(int k=j;k<j+Step;k++){
                CP x=a[k],y=WNK*a[k+Step];
                a[k]=x+y,a[k+Step]=x-y,WNK*=WN;
            }
        }
    }
    if(f==-1) for(int i=0;i<Len;i++) a[i]/=Len;
}
void Cal(int opt){
    int Len=1,bit=0;
    for(int END=2*m;Len<=END;Len<<=1,bit++);
    for(int i=0;i<Len;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    for(int i=m+1;i<Len;i++) Num[i]=0;
    for(int i=0;i<Len;i++) a[i]=(double)Num[i];
    FFT(a,Len,1);
    for(int i=0;i<Len;i++) a[i]=a[i]*a[i];
    FFT(a,Len,-1);
    for(int i=1,END=m*2;i<=END;i++) Ans[i]+=opt*(int)(a[i].real()+0.1);
}
void Solve(int x){
    m=0;Now++;DFS(x,0,0);Cal(1);vis[x]=1;
    for(int j=E.lnk[x];j;j=E.nxt[j])
    if(!vis[E.son[j]]){
        m=0;Now++;Num[0]=0; 
        DFS(E.son[j],x,1),Cal(-1);
        MaxSiz=1e9,Root=0,Al=Siz[E.son[j]];
        Get_Root(E.son[j],x),Solve(Root);
    }
}
int main(){
    #ifndef ONLINE_JUDGE
    freopen("3451.in","r",stdin);
    freopen("3451.out","w",stdout);
    #endif
    scanf("%d",&n);
    for(int i=1;i<n;i++){
        int x,y;scanf("%d%d",&x,&y);x++,y++;
        E.Add(x,y);E.Add(y,x);
    }
    Root=0;Al=n;MaxSiz=1e9;Get_Root(1,0);
    Solve(Root);
    double OUT=n;
    for(int i=1;i<=n*2;i++) if(Ans[i]) OUT+=(double)1.0/(i+1)*Ans[i];
    printf("%.4lf\n",OUT);
    return 0;
}

发表评论

开始在上面输入您的搜索词,然后按回车进行搜索。按ESC取消。

返回顶部