题解 P5252 【[LnOI2019]东京夏日相会】

Hope2075 · 2020-06-20 11:00:06 · 题解

upd:既然已经决定用解析几何暴力求了，就暴力到底吧QWQ

以下是废话

去年比赛时做了这道题，写模拟退火被卡精度，只有30分

当时就想试试用最小覆盖圆的思路，每一步用解析式直接求圆的圆心坐标和半径

然后我想了一下，一般情况需要求3条双曲线的交点（根据性质，应该至少有一个公共交点），也就是解二元二次方程，而且消元很麻烦，就放弃了

今年突然心血来潮，便尝试求了一下，发现并没有那么复杂，就用新的做法A了这道题，而且抢了最优解(2020-06-20 10:32:54 69ms)

下面是具体思路

基本思路

按照点的最小圆覆盖的做法求这些圆的最小圆覆盖，也就是P1742

最小圆覆盖：用的随机增量法，进行n步，每一步加入一个点。如果这个点在之前i-1个点的最小覆盖圆内部，那么就跳过。否则这个点必定在前i个点的最小覆盖圆上，枚举之前i-1个点组成的点对即可

如果把上面的“点”换成“圆”，也成立，所以可以用类似的思路做这道题

复杂度和点的最小覆盖圆相同，参考那道题的证明即可（常数大一些）

另外，在求交点时可以发现一个特殊的性质，就能简化计算

具体做法

以下设三个圆分别为

求一个圆的最小圆覆盖

很明显，取这个圆本身就行

求两个圆的最小圆覆盖

(两个圆的圆心分别是O_1和O_2)

分为两种情况

第一种：两圆半径相等

那么圆心就是两圆圆心连线的中点，半径就是圆心距离的一半加上一个圆的半径

表示出来就是：

（图略，应该很好画）

第二种：两圆半径不相等

这时候圆心到两圆的距离之差就是半径之差，而且圆心更靠近大圆

也就是这样

圆心到中点的距离就是半径差的一半

先把中点坐标求出来，再偏移一下就行

半径就用到某个圆心的距离加上那个圆的半径

表示出来就是

求三个圆的最小圆覆盖

如果某两个圆是包含关系，那么求两个圆的最小圆覆盖即可，而且算法可以保证不会出现这种情况

理论上要分三种情况，但实际有两种可以合并

第一种：三个圆半径均相等

这时候求三角形的外接圆，把半径加上某个圆的半径即可

求外接圆的办法可以参考最小圆覆盖的题解

可以发现这样的圆恰好与这三个圆均内切

也就是这样

第二种：存在两个圆半径不相等

这时候对于某两个圆而言，满足条件的圆的圆心在一条直线或双曲线上

方程就是\sqrt{(x-x_1)^2+(y-y_1)^2}+r_1=\sqrt{(x-x_2)^2+(y-y_2)^2}+r_2

如果两圆半径相等，那么这个方程表示中垂线，否则就表示双曲线（已经限制任意两圆不能是包含关系了）

画出来大概这样（第一个是存在两圆半径相等的情况）

注意图中的深蓝色直线，它有特殊含义，下面会进行叙述

首先，设A_i=(x-x_i)^2+(y-y_i)^2

那么方程可以表示为\sqrt{A_1}+r_1=\sqrt{A_2}+r_2

可以化为\sqrt{A_1}-\sqrt{A_2}=r_2-r_1

分子有理化，再取倒数，得\sqrt{A_1}+\sqrt{A_2}=\frac{A_1-A_2}{r_2-r_1}

相加，得2\sqrt{A_1}=\frac{A_1-A_2}{r_2-r_1}+r_2-r_1

先不要着急平方，再写一个式子

所以\frac{A_1-A_2}{r_2-r_1}+r_2-r_1=\frac{A_1-A_3}{r_3-r_1}+r_3-r_1

可以看出这是一次方程，也就是图中的深蓝色直线

继续化简，设直线方程为ax+by+c=0，那么

这样，所求圆的圆心也在这条直线上

接下来把原来双曲线的方程也化简一下（平方再整理就行，但整理很麻烦）

设双曲线的方程是Ax^2+Bxy+Cy^2+Dx+Ey+F=0

可以代入消元解方程组求坐标

用直线方程消元后得到一个二次方程，用x=\frac{-b\pm\sqrt{b^2-4ac}}{2a}求根，再带入直线方程求出另一个坐标就行

注意直线方程的系数可能为0，需要判断一下

这时会发现有两个解，原因就是平方化简时出现了增根

另一个解恰好是与三个圆分别外切的圆的圆心

但我们需要选出其中一个，这时计算圆心到两个点的距离，比较大小关系就能确定这个根是不是需要的点（目标点距离较大的圆更远）

这样求出来的坐标和半径理论上就没有任何误差（但会受计算机精度限制）

暴力解方程

暴力要彻底

为了防止混淆，令

过程很繁琐，只说一下简单思路和解出的结果，具体过程全部省略

而且大部分解析式我并没有化简，只是猜的结论，但应该是对的

如果手动把系数代进去，会得到一个一元二次方程

解方程即可没这么简单

以下假设消去的是y

首先化简x^2的系数

通过化简，可以得到4(r_1-r_2)^2(T^2-a^2-b^2)

其中T=\left|\begin{array}{cccc} x_1 & y_1 & 1\\ x_2 & y_2 & 1\\ x_3 & y_3 & 1 \end{array}\right|

可以想到，(r_1-r_2)^2可以消掉，否则交换两个点的顺序，结果就会不同

然后化简x的系数

明确了变形的方向就会简单一点，可以先把不含(r_1-r_2)^2的项合并，然后再和剩余部分合并，而且结果应该是交换任意两个点后保持不变的（也就是有对称性）

结果是4(r_1-r_2)^2(x_1m_1n_1+x_2m_2n_2+x_3m_3n_3)

其中m_1=[(x_2-x_3)^2+(y_2-y_3)^2-(r_2-r_3)^2]

求b^2-4ac就行（看起来复杂，实际上比化简系数再乘起来要简单）

根据前面的结果，可以求出两个解对应圆的圆心连线中点的坐标，而两个解则是这个点在直线上移动了一段距离

可以想到对于x和y的方程，后面\pm\sqrt{b^2-4ac}项分别含有因子B和A，并且均含(r_1-r_2)^2因子

按照上面的思路化简即可（最好先消掉因子A或B，这样较为简单）

可以得到b^2-4ac=16B^2m_1m_2m_3(r_1-r_2)^2（对于x）

表示出来的时候，应该是一个\pm，一个\mp

另一个问题：如何确定符号？

可以把半径也暴力求出来，求的时候就会发现根号限定了代入的\pm应该取哪一个（求出的结果应当对称）

（算到最后我没耐心化简了，就根据式子猜了几个结论，代数检验了一下）

最终结果：

要求T>0。如果T=0，则三点共线。如果T<0，那么任意交换一对点后再算即可，或者改变符号也行

可以看出，对于半径两两相等的情况也适用

这样理论上可以运行得很快（用中间变量保留相同的值，每次计算只有一次开根号），而且可以用一些特殊方法做到任意精度

但测试时好像评测机性能不稳定，有些点会慢一点，总用时大约0.2s

求最终答案

和点的最小圆覆盖基本一样，只要把求两点或三点的最小覆盖圆换成相应圆的做法就行

关于特判

我想了一下，应该需要特判三个圆的情况（有时候只和两个圆相切就能盖住三个圆），否则可以手动造几组数据卡掉，但提交的时候没有WA，应该是因为数据随机

代码

这是没卡常的代码，基本就是按照上面的式子写的，用时405ms

（我在纸上推式子时就写的x_1之类的字母，为了省事，计算时没用结构体，调用函数时先把结构体拆开）

69ms的代码就看提交记录吧

主要是IO优化，另外计算过程也优化了一下

手动解方程的方法把对应部分换掉即可，应该比这个好写

#include<iostream>
#include<cmath>
struct C{
    double x,y,r;
    C(){x=0;y=0;r=0;}
    C(double x0,double y0,double r0){x=x0;y=y0;r=r0;}
};
unsigned long long seed=666;
int rand(int l){
    return (seed=(seed*998244353+65537))%l;
}
void swap(double &a,double &b){
    double t=a;a=b;b=t;
}
void swap(C &a,C &b){
    C t=a;a=b;b=t;
}
double dis(double x1,double y1,double x2,double y2){
    return sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
}
bool inside(C a,C b){//a in b
    if(a.r>=b.r)return false;
    if(dis(a.x,a.y,b.x,b.y)<=fabs(a.r-b.r))return true;
    return false;
}
void cal(double x1,double y1,double r1,double x2,double y2,double r2,double &px,double &py,double &r){
    if(r1==r2){
        double rx=(x1+x2)/2;
        double ry=(y1+y2)/2;
        double rr=dis(x1,y1,x2,y2)/2+r1;
        px=rx;py=ry;r=rr;
    }else{
        double mx=(x1+x2)/2;
        double my=(y1+y2)/2;
        double dr=r2-r1;
        double d=dis(x1,y1,x2,y2);
        double rx=mx+(x2-x1)*dr/d/2;
        double ry=my+(y2-y1)*dr/d/2;
        double rr=dis(x1,y1,rx,ry)+r1;
        px=rx;py=ry;r=rr;
    }

}
C cal(C c1,C c2){
    C cr;
    cal(c1.x,c1.y,c1.r,c2.x,c2.y,c2.r,cr.x,cr.y,cr.r);
    return cr;
}
void cal(double x1,double y1,double r1,double x2,double y2,double r2,double x3,double y3,double r3,double &px,double &py,double &r){
    if(r1>r2){//冒泡排序，保证取两个半径不同的圆
        swap(x1,x2);
        swap(y1,y2);
        swap(r1,r2);
    }
    if(r2>r3){
        swap(x2,x3);
        swap(y2,y3);
        swap(r2,r3);
    }
    if(r1>r2){
        swap(x1,x2);
        swap(y1,y2);
        swap(r1,r2);
    }
    if(r1==r3){//如果三个圆的半径相同，那么最大值等于最小值
        double t=2*(x2-x1)*(y3-y1)-2*(x3-x1)*(y2-y1);
        double rx=((x2*x2+y2*y2-x1*x1-y1*y1)*(y3-y1)-(x3*x3+y3*y3-x1*x1-y1*y1)*(y2-y1))/t;
        double ry=-((x2*x2+y2*y2-x1*x1-y1*y1)*(x3-x1)-(x3*x3+y3*y3-x1*x1-y1*y1)*(x2-x1))/t;
        double rr=dis(x1,y1,rx,ry)+r1;
        px=rx;py=ry;r=rr;
    }else{//否则最大的圆和最小的圆半径不同
        double a=x1*(r2-r3)+x2*(r3-r1)+x3*(r1-r2);
        double b=y1*(r2-r3)+y2*(r3-r1)+y3*(r1-r2);
        double c=-0.5*((x1*x1+y1*y1-r1*r1)*(r2-r3)+(x2*x2+y2*y2-r2*r2)*(r3-r1)+(x3*x3+y3*y3-r3*r3)*(r1-r2));

        double A=4*((x1-x3)*(x1-x3)-(r1-r3)*(r1-r3));
        double B=8*(x1-x3)*(y1-y3);
        double C=4*((y1-y3)*(y1-y3)-(r1-r3)*(r1-r3));
        double D=4*((x1+x3)*(r1-r3)*(r1-r3)-(x1-x3)*(x1*x1-x3*x3+y1*y1-y3*y3));
        double E=4*((y1+y3)*(r1-r3)*(r1-r3)-(y1-y3)*(x1*x1-x3*x3+y1*y1-y3*y3));
        double F=(x1*x1-x3*x3+y1*y1-y3*y3)*(x1*x1-x3*x3+y1*y1-y3*y3)-2*(x1*x1+x3*x3+y1*y1+y3*y3)*(r3-r1)*(r3-r1)+(r3-r1)*(r3-r1)*(r3-r1)*(r3-r1);

        double rx1,rx2,ry1,ry2;
        if(fabs(b)>fabs(a)){//取系数绝对值较大的，避免除以0
            double xa=b*b*A-a*b*B+a*a*C;
            double xb=-b*c*B+2*a*c*C+b*b*D-a*b*E;
            double xc=c*c*C-b*c*E+b*b*F;

            double dx=xb*xb-4*xa*xc;
            rx1=(-xb+sqrt(dx))/(2*xa);
            rx2=(-xb-sqrt(dx))/(2*xa);

            ry1=-(c+a*rx1)/b;
            ry2=-(c+a*rx2)/b;

        }else{
            double ya=a*a*C-b*a*B+b*b*A;
            double yb=-a*c*B+2*b*c*A+a*a*E-b*a*D;
            double yc=c*c*A-a*c*D+a*a*F;

            double dy=yb*yb-4*ya*yc;
            ry1=(-yb+sqrt(dy))/(2*ya);
            ry2=(-yb-sqrt(dy))/(2*ya);

            rx1=-(c+b*ry1)/a;
            rx2=-(c+b*ry2)/a;

        }

        double d11=dis(rx1,ry1,x1,y1);
        double d13=dis(rx1,ry1,x3,y3);

        double d21=dis(rx2,ry2,x1,y1);

        if(d11>d13){
            double rr1=d11+r1;
            px=rx1;py=ry1;r=rr1;
        }else{
            double rr2=d21+r1;
            px=rx2;py=ry2;r=rr2;
        }
    }
}
C cal(C c1,C c2,C c3){
    C cr;
    cal(c1.x,c1.y,c1.r,c2.x,c2.y,c2.r,c3.x,c3.y,c3.r,cr.x,cr.y,cr.r);
    return cr;
}
int n;
C list[50007];
C cur;
int main(){
    std::ios::sync_with_stdio(false);
    std::cin>>n;
    for(int i=1;i<=n;i++){
        double x,y,r;
        std::cin>>x>>y>>r;
        list[i]=C(x,y,r);
    }
    for(int i=n;i>1;i--){
    //需要打乱顺序，这样才能保证复杂度
        int p=rand(i);
        swap(list[i],list[p+1]);
    }
    cur=list[1];
    for(int i=2;i<=n;i++){
        if(!inside(list[i],cur)){
            cur=list[i];
            for(int j=1;j<i;j++){
                if(!inside(list[j],cur)){
                    cur=cal(list[i],list[j]);
                    for(int k=1;k<j;k++){
                        if(!inside(list[k],cur)){
                            cur=cal(list[i],list[j],list[k]);
                        }
                    }
                }
            }
        }
    }
    std::cout<<cur.x<<" "<<cur.y<<" "<<cur.r<<std::endl;
}