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

· · 题解

原题链接

闲话：模拟赛出的，然后考场推的三元函数，又臭又长还是错的，赛后去掉一元后十分简单，写篇题解记录下。

可能比较长，式子比较多，但还是挺好理解的。

题目分析

人话题意：给出一棵树或基环树，每个节点有一个二维平面上的店，求一条路径上点的线性回归（定一条直线 l，最小化 \sum_i\operatorname{distance}(P_i,l)^2）。

\text{Part 1}：如何求线性回归

首先抛开什么基环树，先假设你已经知道了点集 \{P_i\}。利用点到直线距离公式：对于点 P(x_0,y_0) 和直线 l:Ax+By+C=0，有

\operatorname{distance}(P,l)=\frac{|Ax_0+By_0+C|}{\sqrt{A^2+B^2}}

在这个问题中，我们不妨令直线方程为 Ax+y+B=0（一条竖直的线可以看做 A 取很大的数），所以我们写出距离的平方和：

\begin{aligned} \sum_i\operatorname{distance}(P_i,l)^2&=\sum_i\frac{(Ax_i+y_i+B)^2}{A^2+1}\\ &=\frac{1}{A^2+1}\sum_i(Ax_i+y_i+B)^2\\ &=\frac{1}{A^2+1}\sum_i(A^2x_i^2+B^2+2ABx_i+2Ax_iy_i+2By_i+y_i^2) \end{aligned}

显然这个东西化简后形如 \dfrac{uA^2+vB^2+wAB+pA+qB+r}{A^2+1}。所以我们可以把距离和看做一个关于 A,B 的二元函数 f(A,B)。

所以我们现在要求 f(A,B) 的最值，考虑二分一个值 z，现在要解决的问题就是判断 f(A,B)=z 是否有解。将式子变形：

\begin{aligned} \frac{uA^2+vB^2+wAB+pA+qB+r}{A^2+1}&=z\\ uA^2+vB^2+wAB+pA+qB+r&=zA^2+z\\ (u-z)A^2+vB^2+wAB+pA+qB+(r-z)&=0 \end{aligned}

这又是一个二元二次函数，我们叫它 g(A,B)，现在就是要判断 g(A,B)=0 是否有解。我们将从图像上考虑。

我们知道二元二次方程与圆锥曲线有着密不可分的关系，事实上，对于一个二元二次方程 Ax+By+Cxy+Dx+Ey+F=0，我们有如下的结论：

\Delta=C^2-4AB,\begin{cases} \text{图像是双曲线}&,\Delta>0\\ \text{图像是抛物线}&,\Delta=0\\ \text{图像是椭圆}&,\Delta<0 \end{cases}

如果这个方程真正有解的话。（一个巧妙的证明）

接下来，借助一些画图工具，我们可以把这个结论搬到三维空间中：

现在考虑 \Delta<0 的时候怎么算最低点。实际上因为这个函数图像的凸性，我们只需要令 x,y 的偏导都为 0 即可，也就是说：

\frac{\partial g}{\partial A}=\frac{\partial g}{\partial B}=0

即

\begin{cases} 2(u-z)A+wB+p=0\\ wA+2vB+q=0 \end{cases}

得到的是一个二元一次方程组，大家都会解。算出最低点 (A,B) 后带入 g(A,B) 即可。

\text{Part 2}：回到基环树上

看看我们最开始的 f(A,B)，发现所有系数都是一个和式，满足可加性。所以我们只要维护 x_i^2,y_i^2,x_iy_i,1 的和即可。

在树上可以通过树上前缀和实现。放到基环树上其实也一样。随便取一棵生成树，然后判断询问的路径是否经过环，如果经过环，就算一下走环的另外一边的答案。

所以说这题的基环树其实除了增加代码难度之外没什么意义。

代码实现

代码逻辑清晰，张弛有度，长而不乱，让人神清气爽、心旷神怡。

#include <bits/stdc++.h>
using namespace std;

/* 省略快读快写，即下面的 qin、qout */

int n,m,q;
vector<int> g1[100005],g2[100005];

struct Data
{
    int A2,B2,AB,A,B,C;
    Data()=default;
    Data(int _a2,int _b2,int _ab,int _a,int _b,int _c):
        A2(_a2),B2(_b2),AB(_ab),A(_a),B(_b),C(_c){}
    inline Data &operator +=(const Data &rhs)
        { return A2+=rhs.A2,B2+=rhs.B2,AB+=rhs.AB,A+=rhs.A,B+=rhs.B,C+=rhs.C,*this; }
    inline Data &operator -=(const Data &rhs)
        { return A2-=rhs.A2,B2-=rhs.B2,AB-=rhs.AB,A-=rhs.A,B-=rhs.B,C-=rhs.C,*this; }
    inline Data operator +(const Data &rhs)const
        { return Data(*this)+=rhs; }
    inline Data operator -(const Data &rhs)const
        { return Data(*this)-=rhs; }
}a[100005],sum[100005],cir;

bool vis[100005],isc[100005];
int top[100005];
void buildTree(int u)
{
    vis[u]=true;
    for(auto &v: g1[u])
    {
        if(vis[v]) continue;
        g2[u].push_back(v);
        g2[v].push_back(u);
        buildTree(v);
    }
}
bool flag;
stack<int> st;
void findCirc(int u,int fa=0)
{
    if(vis[u])
    {
        while(st.top()!=u)
            isc[st.top()]=true,st.pop();
        isc[u]=flag=true;
        return;
    }
    vis[u]=true,st.push(u);
    for(auto &v: g1[u])
    {
        if(v==fa) continue;
        findCirc(v,u);
        if(flag) return;
    }
    st.pop();
}
void getTop(int u,int fa=0,int tp=0)
{
    top[u]=(isc[u]?tp=u:tp);
    for(auto &v: g1[u])
        if(v!=fa && !isc[v]) getTop(v,u,tp);
}

int f[100005][20],dep[100005];
void dfs(int u,int fa=0)
{
    sum[u]=sum[fa]+a[u],f[u][0]=fa,dep[u]=dep[fa]+1;
    for(int i=1;(1<<i)<=dep[u];i++)
        f[u][i]=f[f[u][i-1]][i-1];
    for(auto &v: g2[u])
        if(v!=fa) dfs(v,u);
}
inline int lca(int x,int y)
{
    if(dep[x]<dep[y]) swap(x,y);
    for(int i=17;i>=0;i--)
        if(dep[f[x][i]]>=dep[y])
            x=f[x][i];
    if(x==y) return x;
    for(int i=17;i>=0;i--)
        if(f[x][i]!=f[y][i])
            x=f[x][i],y=f[y][i];
    return f[x][0];
}
inline Data getsum(int u,int v)
{
    int l=lca(u,v);
    return sum[u]+sum[v]-sum[l]-sum[f[l][0]];
}

inline bool check(const Data &c,LD x)
{
    LD A2=c.A2,B2=c.B2,AB=c.AB,A=c.A,B=c.B,C=c.C;
    A2-=x,C-=x;
    if(AB*AB-4*A2*B2>=0) return true;
    LD a=(-2*B2*A+AB*B)/(4*A2*B2-AB*AB),
       b=(-2*A2*B+AB*A)/(4*A2*B2-AB*AB);
    return A2*a*a+B2*b*b+AB*a*b+A*a+B*b+C<=0;
}
inline LD solve(Data c)
{
    LD l=0,r=min(c.A,c.C),mid;
    while(r-l>1e-4)
    {
        mid=(l+r)/2;
        if(check(c,mid)) r=mid;
        else l=mid;
    }
    return l;
}
int main()
{
    int x,y;
    qin>>n>>m;
    for(int i=1;i<=n;i++)
    {
        qin>>x>>y,x+=2,y+=2;
        a[i].A2=x*x,a[i].B2=1,a[i].AB=2*x;
        a[i].A=2*x*y,a[i].B=2*y,a[i].C=y*y;
    }
    for(int i=1;i<=m;i++)
    {
        qin>>x>>y;
        g1[x].push_back(y);
        g1[y].push_back(x);
    }
    buildTree(1),dfs(1);
    if(m==n)
    {
        memset(vis,0,sizeof(vis)),findCirc(1);
        for(int i=1;i<=n;i++)
            if(isc[i]) cir+=a[i],getTop(i);
    }
    qin>>q;
    while(q--)
    {
        qin>>x>>y;
        LD ans=INFINITY;
        Data zt=getsum(x,y);
        ans=min(ans,solve(zt));
        if(n==m && top[x]!=top[y])
        {
            int tx=top[x],ty=top[y];
            Data yt=getsum(tx,ty);
            zt+=cir-yt-yt+a[tx]+a[ty];
            ans=min(ans,solve(zt));
        }
        printf("%.6Lf\n",ans);
    }
    return 0;
}