题解 洛谷 P3340 【[ZJOI2014]星系调查】

根据题意，发现题目中的图，其实就是一颗树或者是一颗基环树，每个节点上有一个点对\((x,y)\)，每次询问为给定端点，找一条直线到端点间的所有点的距离之和最小。

设这条直线为\(y=kx+b\)，根据点到直线公式得，我们要求\(\sum\limits_{i=1}^n \frac{(kx_i-y_i+b)^2}{k^2+1}\)的最小值，其中\(n\)为路径上的点的个数。

然后开始处理这个式子：

\[\begin{aligned} &\sum\frac{(kx_i-y_i+b)^2}{k^2+1} \\ =&\sum\frac{(kx_i-y_i)^2+2(kx_i-y_i)b+b^2}{k^2+1} \\ =&\sum\frac{k^2x_i^2-2kx_iy_i+y_i^2+2kbx_i-2by_i+b^2}{k^2+1} \\ =&\frac{\sum k^2x_i^2-\sum 2kx_iy_i+\sum y_i^2+\sum 2kbx_i-\sum 2by_i+nb^2}{k^2+1} \\ \end{aligned} \]

发现对于直线\(y=kx+b\)，\(k\)和\(b\)的取值是互不影响的，所以可以先对\(b\)进行处理，使其取到能让式子达到最小值的值，然后代入原式，再求\(k\)能让式子达到最小值的值。

先对含有\(b\)的项进行分类，得：

\[\frac{nb^2+\sum (2kx_i-2y_i)b+\sum (k^2x_i^2-2kx_iy_i+y_i^2)}{k^2+1} \]

考虑\(b\)作为变量，发现其为开口向上的二次函数，当\(b= \frac{\sum (2y_i-2kx_i)}{2n}=\frac{\sum y_i}{n}-\frac{\sum kx_i}{n}\)，发现\(\frac{\sum x_i}{n}\)的意义就是对于所有\(x_i\)的平均值，所以设\(\bar x=\frac{\sum x_i}{n},\bar y=\frac{\sum y_i}{n}\)，将\(b=\bar y-k \bar x\)代入原式，得：

\[\frac{\sum k^2x_i^2-\sum 2kx_iy_i+\sum y_i^2+\sum 2k(\bar y-k \bar x)x_i-\sum 2(\bar y-k \bar x)y_i+n(\bar y-k \bar x)^2}{k^2+1} \\ \]

然后像上面一样，对含有\(k\)的项进行分类，得：

\[\frac{\sum (x_i^2-2x_i \bar x+ \bar x^2)k^2+\sum(-2x_iy_i+2x_i \bar y+2 \bar x y_i -2 \bar x \bar y)k+\sum (y_i^2-2y_i \bar y+ \bar y^2)}{k^2+1} \]

为方便接下来的处理，设出\(k\)的系数：

\[\begin{aligned} &A=\sum (x_i^2-2x_i \bar x+ \bar x^2)=\sum x_i^2-\frac{(\sum x_i)^2}{n}\\ &B=\sum(-2x_iy_i+2x_i \bar y+2 \bar x y_i -2 \bar x \bar y)=-2\sum x_iy_i+\frac{2\sum x_i\sum y_i}{n}\\ &C=\sum (y_i^2-2y_i \bar y+ \bar y^2)=\sum y_i^2-\frac{(\sum y_i)^2}{n}\\ \end{aligned} \]

设原式的值为\(S\)，得：

\[\begin{aligned} &S=\frac{Ak^2+Bk+C}{k^2+1}\\ &Sk^2+S=Ak^2+Bk+C\\ &(A-S)k^2+Bk+C-S=0\\ \end{aligned} \]

发现得到一个关于\(k\)的一元二次方程，因为\(k\)一定会有合法解，所以得：

\[\begin{aligned} \Delta&=B^2-4(A-S)(C-S) \geqslant 0\\ &=B^2-4AC+4AS+4CS-4S^2 \geqslant 0\\ &=-4S^2+4(A+C)S+B^2-4AC \geqslant 0\\ \end{aligned} \]

发现得到一个关于\(S\)的一元二次不等式，将其看作开口向下的二次函数，得其值必须大于等于\(0\)，所以\(S\)的最小取值即为其方程的较小根，得：

\[\begin{aligned} &S=\frac{A+C-\sqrt{A^2+B^2+C^2-2AC}}{2}\\ &A=\sum x_i^2-\frac{(\sum x_i)^2}{n}\\ &B=-2\sum x_iy_i+\frac{2\sum x_i\sum y_i}{n}\\ &C=\sum y_i^2-\frac{(\sum y_i)^2}{n}\\ \end{aligned} \]

然后维护每个点到根节点\(\sum x_i,\sum y_i,\sum x_i^2,\sum y_i^2,\sum x_iy_i,n\)的和，然后就可以通过树上差分求解了。对于基环树的情况，若两点在一个子树内，就直接树上差分，若两点不在一个子树内，需要经过环，则从经过环的两种方式得出的两个值取较小值作为答案。

\(code:\)

#include<bits/stdc++.h>
#define maxn 400010
using namespace std;
template<typename T> inline void read(T &x)
{
    x=0;char c=getchar();bool flag=false;
    while(!isdigit(c)){if(c=='-')flag=true;c=getchar();}
    while(isdigit(c)){x=(x<<1)+(x<<3)+(c^48);c=getchar();}
    if(flag)x=-x;
}
int n,m,q,tot;
int de[maxn],f[maxn][25],fa[maxn],a[maxn],b[maxn];
int cir[maxn],bel[maxn],pos[maxn];
bool flag;
bool tag[maxn];
struct edge
{
    int to,nxt;
}e[maxn];
int head[maxn],edge_cnt;
void add(int from,int to)
{
    e[++edge_cnt]=(edge){to,head[from]};
    head[from]=edge_cnt;
}
struct node
{
    int x,y,xd,yd,mul,num;
    void insert(int a,int b)
    {
        x+=a,y+=b,xd+=a*a,yd+=b*b,mul+=a*b,num++;
    }
    double calc()
    {
        double A=xd-(double)x*x/num,B=-2*mul+2*(double)x*y/num,C=yd-(double)y*y/num;
        return (A+C-sqrt(A*A+B*B+C*C-2*A*C))/2;
    }
}s[maxn],sum[maxn];
node operator +(const node &a,const node &b)
{
    return (node){a.x+b.x,a.y+b.y,a.xd+b.xd,a.yd+b.yd,a.mul+b.mul,a.num+b.num};
}
node operator -(const node &a,const node &b)
{
    return (node){a.x-b.x,a.y-b.y,a.xd-b.xd,a.yd-b.yd,a.mul-b.mul,a.num-b.num};
}
int find(int x)
{
    return fa[x]==x?x:fa[x]=find(fa[x]);
}
void dfs(int x,int fath,int col)
{
    f[x][0]=fath,de[x]=de[fath]+1,bel[x]=col,s[x]=s[fath],s[x].insert(a[x],b[x]);
    for(int i=1;i<=18;++i) f[x][i]=f[f[x][i-1]][i-1];
    for(int i=head[x];i;i=e[i].nxt)
    {
        int y=e[i].to;
        if(y==fath||tag[y]) continue;
        dfs(y,x,col);
    }
}
int lca(int x,int y)
{
    if(de[x]<de[y]) swap(x,y);
    for(int i=18;i>=0;--i)
        if(de[f[x][i]]>=de[y])
            x=f[x][i];
    if(x==y) return x;
    for(int i=18;i>=0;--i)
        if(f[x][i]!=f[y][i])
            x=f[x][i],y=f[y][i];
    return f[x][0];
}
void dfs_cir(int x,int fath)
{
    if(tag[x]&&x!=cir[1])
    {
        flag=true;
        return;
    }
    for(int i=head[x];i;i=e[i].nxt)
    {
        int y=e[i].to;
        if(y==fath) continue;
        if(tag[x]&&tag[y]) continue;
        cir[++tot]=y,dfs_cir(y,x);
        if(flag) return;
        tot--;
    }
}
double solve(int x,int y)
{
    int anc=lca(x,y);
    return (s[x]+s[y]-s[anc]-s[f[anc][0]]).calc();
}
int main()
{
    read(n),read(m);
    for(int i=1;i<=n;++i) read(a[i]),read(b[i]),fa[i]=i;
    for(int i=1;i<=m;++i)
    {
        int x,y;
        read(x),read(y),add(x,y),add(y,x);
        if(find(x)==find(y)) cir[++tot]=x,tag[x]=tag[y]=true;
        else fa[find(x)]=find(y);
    }
    if(m==n-1) dfs(1,0,0);
    else
    {
        dfs_cir(cir[1],0);
        for(int i=1;i<=tot;++i)
        {
            int x=cir[i];
            pos[x]=i,tag[x]=true;
            sum[i]=sum[i-1],sum[i].insert(a[x],b[x]);
        }
        for(int i=1;i<=tot;++i) dfs(cir[i],0,cir[i]);
    }
    read(q);
    while(q--)
    {
        int x,y;
        read(x),read(y);
        if(m==n-1) printf("%.6lf\n",solve(x,y));
        else
        {
            if(bel[x]==bel[y]) printf("%.6lf\n",solve(x,y));
            else
            {
                int px=pos[bel[x]],py=pos[bel[y]];
                node val=s[x]+s[y]-s[bel[x]]-s[bel[y]];
                if(px>py) swap(px,py);
                printf("%.6lf\n",min((val+sum[py]-sum[px-1]).calc(),(val+sum[tot]-sum[py-1]+sum[px]).calc()));
            }
        }
    }
    return 0;
}