[BZOJ3528][Zjoi2014][数学]星系调查

上一发大神的推导

后面代入求根公式就行。
因为数据只有M=N或M=N-1两种情况，所以这张图要不是一棵树要不只有一个环，是一棵树的情况很好考虑，有一个环的话要分类讨论。

感觉ZJOI2014DAY1的题代码都好长，然而杜教AK了…….

感觉这张图上一坨水印看着不爽，不知道怎么去……

#include <cstdio>
#include <cmath>
#include <string>
#include <cstring>
#include <iostream>
#define N 100010

using namespace std;

typedef long long ll;

int n,m,q,cnt,u,v,l,Root,ccnt;
int G[N],X[N],Y[N],dpt[N],sta[N],tp,V[N],inC[N],onT[N],r[N];
int fa[N][30];
ll sgx[N],sgy[N],sgx_2[N],sgy_2[N],sgxy[N];
ll csgx,csgy,csgx_2,csgy_2,csgxy;
double A,B,C;

struct edge{
    int t,nx;
}E[N<<2];

inline char gc(){
    static char buf[100000],*p1=buf,*p2=buf;
    if(p1==p2){
        p2=(p1=buf)+fread(buf,1,100000,stdin);
        if(p1==p2) return EOF;
    }
    return *p1++;
}

inline void reaD(int &x){
    char Ch=gc();x=0;
    for(;Ch>'9'||Ch<'0';Ch=gc());
    for(;Ch>='0'&&Ch<='9';x=x*10+Ch-'0',Ch=gc());
}

inline void InserT(int x,int y){
    E[++cnt].t=y;E[cnt].nx=G[x];G[x]=cnt;
    E[++cnt].t=x;E[cnt].nx=G[y];G[y]=cnt;
}

void dfs_t(int x,int f,int flg=0){
    onT[x]=flg;
    dpt[x]=dpt[f]+1;fa[x][0]=f;
    sgxy[x]=sgxy[f]+X[x]*Y[x];
    sgx[x]=sgx[f]+X[x],sgy[x]=sgy[f]+Y[x];
    sgx_2[x]=sgx_2[f]+X[x]*X[x],sgy_2[x]=sgy_2[f]+Y[x]*Y[x];
    for(int i=1;i<=20;i++) fa[x][i]=fa[fa[x][i-1]][i-1];
    for(int i=G[x];i;i=E[i].nx)
        if(E[i].t!=f&&!inC[E[i].t]) dfs_t(E[i].t,x,flg);
}

void dfs_g(int x,int f){
    V[x]=1;
    dpt[x]=dpt[f]+1;
    sgxy[x]=sgxy[f]+X[x]*Y[x];
    sgx[x]=sgx[f]+X[x],sgy[x]=sgy[f]+Y[x];
    sgx_2[x]=sgx_2[f]+X[x]*X[x],sgy_2[x]=sgy_2[f]+Y[x]*Y[x];
    for(int i=G[x];i;i=E[i].nx)
        if(E[i].t!=f&&inC[E[i].t]&&!V[E[i].t]) dfs_g(E[i].t,x);
}

inline void swap(int &x,int &y){
    int z=x;x=y;y=z;
}

int lca(int x,int y){
    if(dpt[x]<dpt[y]) swap(x,y);
    for(int i=20;~i;i--)
        if(dpt[fa[x][i]]>=dpt[y]) x=fa[x][i];
    if(x==y) return x;
    for(int i=20;~i;i--)
        if(fa[x][i]!=fa[y][i]) x=fa[x][i],y=fa[y][i];
    return fa[x][0];
}

void Calc_t(int u,int v){
    l=lca(u,v);
    ll nn=dpt[u]+dpt[v]-2*dpt[l]+1,_sgx=sgx[u]+sgx[v]-2*sgx[l]+X[l],_sgy=sgy[u]+sgy[v]-2*sgy[l]+Y[l],
       _sgx_2=sgx_2[u]+sgx_2[v]-2*sgx_2[l]+X[l]*X[l],_sgy_2=sgy_2[u]+sgy_2[v]-2*sgy_2[l]+Y[l]*Y[l],_sgxy=sgxy[u]+sgxy[v]-2*sgxy[l]+Y[l]*X[l];
    double _x=(double)_sgy/nn,_y=(double)_sgx/nn;
    A=nn*_y*_y-2*_y*_sgx+_sgx_2;
    B=2*_x*_sgx-2*nn*_x*_y+2*_y*_sgy-2*_sgxy;
    C=nn*_x*_x-2*_x*_sgy+_sgy_2;
    printf("%.5lf\n",(A+C-sqrt((A-C)*(A-C)+B*B))/2.0+1e-7);
}

void Calc_c(int u,int v){
    if(dpt[u]<dpt[v]) swap(u,v);
    ll nn=dpt[u]-dpt[v]+1,_sgx=sgx[u]-sgx[v]+X[v],_sgy=sgy[u]-sgy[v]+Y[v],
       _sgx_2=sgx_2[u]-sgx_2[v]+X[v]*X[v],_sgy_2=sgy_2[u]-sgy_2[v]+Y[v]*Y[v],_sgxy=sgxy[u]-sgxy[v]+Y[v]*X[v];
    double _x=(double)_sgy/nn,_y=(double)_sgx/nn,Ans;
    A=nn*_y*_y-2*_y*_sgx+_sgx_2;
    B=2*_x*_sgx-2*nn*_x*_y+2*_y*_sgy-2*_sgxy;
    C=nn*_x*_x-2*_x*_sgy+_sgy_2;
    Ans=(A+C-sqrt((A-C)*(A-C)+B*B))/2.0;
    nn=ccnt-nn+2,_sgx=csgx-_sgx+X[u]+X[v],_sgy=csgy-_sgy+Y[u]+Y[v],_sgx_2=csgx_2-_sgx_2+X[u]*X[u]+X[v]*X[v],_sgy_2=csgy_2-_sgy_2+Y[u]*Y[u]+Y[v]*Y[v];_sgxy=csgxy-_sgxy+X[u]*Y[u]+X[v]*Y[v];
    _x=(double)_sgy/nn,_y=(double)_sgx/nn,Ans;
    A=nn*_y*_y-2*_y*_sgx+_sgx_2;
    B=2*_x*_sgx-2*nn*_x*_y+2*_y*_sgy-2*_sgxy;
    C=nn*_x*_x-2*_x*_sgy+_sgy_2;
    printf("%.5lf\n",min((A+C-sqrt((A-C)*(A-C)+B*B))/2.0,Ans)+1e-7);
}

void Calc_t_c(int u,int v){
    ll n=0,sgx1=0,sgy1=0,sgx21=0,sgy21=0,sgxy1=0;
    if(onT[u]) n=dpt[u],sgx1=sgx[u],sgy1=sgy[u],sgx21=sgx_2[u],sgy21=sgy_2[u],sgxy1=sgxy[u],u=r[onT[u]];
    if(onT[v]) n+=dpt[v],sgx1+=sgx[v],sgy1+=sgy[v],sgx21+=sgx_2[v],sgy21+=sgy_2[v],sgxy1+=sgxy[v],v=r[onT[v]];

    if(dpt[u]<dpt[v]) swap(u,v);
    ll nn=dpt[u]-dpt[v]+1,_sgx=sgx[u]-sgx[v]+X[v],_sgy=sgy[u]-sgy[v]+Y[v],
       _sgx_2=sgx_2[u]-sgx_2[v]+X[v]*X[v],_sgy_2=sgy_2[u]-sgy_2[v]+Y[v]*Y[v],_sgxy=sgxy[u]-sgxy[v]+Y[v]*X[v];
    double _x=(double)(_sgy+sgy1)/(nn+n),_y=(double)(_sgx+sgx1)/(nn+n),Ans;
    A=(nn+n)*_y*_y-2*_y*(_sgx+sgx1)+(_sgx_2+sgx21);
    B=2*_x*(_sgx+sgx1)-2*(nn+n)*_x*_y+2*_y*(_sgy+sgy1)-2*(_sgxy+sgxy1);
    C=(nn+n)*_x*_x-2*_x*(_sgy+sgy1)+(_sgy_2+sgy21);
    Ans=(A+C-sqrt((A-C)*(A-C)+B*B))/2.0;
    nn=ccnt-nn+2,_sgx=csgx-_sgx+X[u]+X[v],_sgy=csgy-_sgy+Y[u]+Y[v],_sgx_2=csgx_2-_sgx_2+X[u]*X[u]+X[v]*X[v],_sgy_2=csgy_2-_sgy_2+Y[u]*Y[u]+Y[v]*Y[v];_sgxy=csgxy-_sgxy+X[u]*Y[u]+X[v]*Y[v];
    _x=(double)(_sgy+sgy1)/(nn+n),_y=(double)(_sgx+sgx1)/(nn+n);
    A=(nn+n)*_y*_y-2*_y*(_sgx+sgx1)+(_sgx_2+sgx21);
    B=2*_x*(_sgx+sgx1)-2*(nn+n)*_x*_y+2*_y*(_sgy+sgy1)-2*(_sgxy+sgxy1);
    C=(nn+n)*_x*_x-2*_x*(_sgy+sgy1)+(_sgy_2+sgy21);
    printf("%.5lf\n",min((A+C-sqrt((A-C)*(A-C)+B*B))/2.0,Ans)+1e-7);
}

void WorkOnTree(){
    dfs_t(1,0);
    for(reaD(q);q;q--)
        reaD(u),reaD(v),Calc_t(u,v);
}

bool Getcir(int x,int f){
    if(V[x]){
        while(tp){
            int k=sta[tp--];inC[k]=1;ccnt++;
            csgx+=X[k],csgy+=Y[k],csgx_2+=X[k]*X[k];
            csgy_2+=Y[k]*Y[k],csgxy+=X[k]*Y[k];
            if(k==x) break;
        }
        return true;
    }
    V[x]=1;
    sta[++tp]=x;
    for(int i=G[x];i;i=E[i].nx)
        if(E[i].t!=f){
            if(Getcir(E[i].t,x))return true;
        }else f=0;
    tp--;
    return false;
}

void WorkOnGraph(){
    Getcir(1,0);cnt=0;
    for(int k=1;k<=n;k++)
        if(inC[k])for(int i=G[k];i;i=E[i].nx)
            if(!inC[E[i].t]) dfs_t(E[i].t,0,E[i].t),r[E[i].t]=k;
    memset(V,0,sizeof(V));
    for(int k=1;k<=n;k++)
        if(inC[k]){dfs_g(k,0);Root=k;break;}
    for(reaD(q);q;q--){
        reaD(v);reaD(u);
        if(onT[v]==onT[u]&&!inC[v]&&!inC[u]) Calc_t(u,v);
        else if(inC[v]&&inC[u]) Calc_c(u,v);
        else Calc_t_c(u,v);
    }
}

int main(){
    reaD(n);reaD(m);
    for(int i=1;i<=n;i++) reaD(X[i]),reaD(Y[i]);
    for(int i=1;i<=m;i++) reaD(u),reaD(v),InserT(u,v);
    if(m==n-1) WorkOnTree();
    else WorkOnGraph();
}