题解：P16549 [ICPC 2026 LAC] Crop Circles

· · 题解

题意简述

给定 n 个圆盘，求圆周完全包含在这些圆盘并集中的圆的最大半径。只要求覆盖圆周，圆内部可以存在未被覆盖的区域。

其中 n\le 40，圆心与半径均为整数，圆心两两不同。

解题思路

记第 i 个圆盘的圆心为 c_i，半径为 R_i，全部圆盘的并集为 U。直接选择任意一个原圆的圆周都是合法方案，因此答案至少为 R_{\max}=\max R_i。

下面只考虑半径 r>R_{\max} 的候选圆。关键是证明：如果最优半径大于 R_{\max}，就一定可以用三个原圆之间的交点确定这个最优圆。

设某个可行圆的圆心为 p，圆周为 C。将 C 与 U 的边界的交点称为接触点。其余圆周点都严格位于 U 内部，具有一定的移动余量。

接触点不可能位于某个原圆盘的内部，否则它也位于 U 内部。因此，覆盖接触点 z 的圆盘都满足 |z-c_i|=R_i。

考虑候选圆在 z 附近的两个方向。设 e=(z-p)/r 为候选圆在 z 处的单位外法向量，u 为单位切向量，用有正负的弧长 s 参数化附近的圆周。由圆的局部展开可得：

\begin{aligned} \gamma(s) & =z+su-\frac{s^2}{2r}e+O(s^3) \\ |\gamma(s)-c_i|^2-R_i^2 & =2(z-c_i)\cdot u\,s+\left(1-\frac{(z-c_i)\cdot e}{r}\right)s^2+O(s^3) \end{aligned}

由于 |z-c_i|=R_i<r，二次项系数严格大于 0。

若一次项系数为 0，则在足够靠近 z 的两侧，候选圆周都位于这个圆盘外部。几何上，这是因为候选圆的半径更大，在切点附近弯曲得更缓，不能沿着半径更小的圆盘内部延伸。

若一次项系数不为 0，则这个圆盘恰好覆盖 z 附近的一侧。为了同时覆盖两侧，必然需要至少两个圆盘，并且它们的一次项系数异号。

所以，每个接触点都是两个原圆的交点，并且不在任何原圆盘内部。称这样的交点为外露交点。这里的外露同时包含并集外边界和孔洞边界，不能只保留最外层轮廓。

上述异号条件还说明：在保持经过 z 的前提下，略微改变候选圆，两个圆盘仍分别覆盖 z 附近的两侧。因此，经过这个接触点的局部覆盖关系是稳定的。

现在证明最优圆至少有三个不同的接触点。由于所有可行圆都处于 U 的包围盒内，圆心和半径有界；覆盖条件又是闭合的，所以最大半径能够取到。

当 r>R_{\max} 时，候选圆不可能与任何原圆重合，接触点数量有限。若接触点不足三个，可以分别讨论：

后两种移动不会破坏接触点附近的覆盖。去掉这些点的小邻域后，剩余部分是严格位于 U 内部的紧集，仍有统一的正距离余量。因此，只要移动足够小，剩余圆周也继续被覆盖。

三种情况都能得到半径更大的可行圆，与最优性矛盾。所以，最优半径若大于 R_{\max}，最优圆必然经过至少三个不同的外露交点。圆上的三个不同点不共线，可以唯一确定这个圆。

据此，先枚举所有原圆对的交点，删除严格位于其他圆盘内部的点，并合并重复点。之后枚举三个不共线的保留点，求它们的外接圆，只验证半径大于当前答案的候选圆。初始答案取 R_{\max}，便不会遗漏最优解。

还需要准确判断一个候选圆的整条圆周是否被覆盖，不能只检查三个确定点，也不能检查固定数量的采样点。

对圆心为 p、半径为 r 的候选圆，考虑一个圆心为 c_i、半径为 R_i 的原圆盘。设 d=|c_i-p|，从 p 指向 c_i 的方向角为 \alpha。候选圆周上方向角为 \theta 的点被这个圆盘覆盖，当且仅当：

d^2+r^2-2dr\cos(\theta-\alpha)\le R_i^2

若 d+r\le R_i，整个候选圆周都被这个圆盘覆盖，可以直接返回合法。

若两个圆外离，或者候选圆包住原圆盘但两条圆周不相交，则这个圆盘不能覆盖正长度的候选圆弧。相切产生的孤立点也不能填补任何正长度的空缺，可以忽略。

其余情况下，两圆相交，覆盖的是以 \alpha 为中心、半角为 \beta 的闭区间：

\beta=\arccos\frac{d^2+r^2-R_i^2}{2dr}

将跨过 0 的角度区间拆成两段，统一放入 [0,2\pi]。按左端点排序，扫描并维护已经连续覆盖到的最右端点；如果出现空隙则不合法，否则检查是否覆盖到 2\pi。

代码在区间判定前做了两个必要条件检查：候选圆不能超出全部圆盘的包围盒，上、下、左、右四个极值点也必须被覆盖。它们只用于提前排除明显不合法的候选圆，最终是否覆盖整条圆周仍由角度区间判定。

计算使用 long double。对反余弦的参数截断到 [-1,1]，避免浮点舍入导致越界；区间拼接允许很小的舍入误差。两圆交点的计算使用圆心连线上的投影及垂直方向的偏移，相切时重复产生的交点会在加入时去重。

设外露交点数为 m。枚举圆对并检查交点是否被其他圆盘遮住，需要 O(n^3)；枚举三点并检查角度覆盖，需要 O(m^3n\log n)。

这里 m 实际上为 O(n)，而不是任意排列下两两交点的 O(n^2)。可以用幂图解释这个界：定义点 x 对第 i 个圆的幂为 h_i(x)=|x-c_i|^2-R_i^2，则并集边界满足 \min_i h_i(x)=0。外露交点处至少有两个圆同时取得最小值。

从各个 h_i 中减去公共的 |x|^2 后，比较的是 n 个仿射函数，它们的最小值区域组成平面上的幂图。每个非空区域都是凸的，这个平面细分总共只有 O(n) 条边和顶点。外露交点若位于某条边内部，就同时位于这条直线和对应原圆上，每条边至多产生两个；位于幂图顶点的外露交点也只有 O(n) 个。因此，不同外露交点总数为 O(n)。

总时间复杂度为 O(n^4\log n)。有效保留的数据仅需 O(n+m) 空间，但代码中的交点数组按所有圆对最多产生两个交点的上界预留容量，避免依赖线性界的具体常数，因此当前实现的空间复杂度为 O(n^2)。

参考代码

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

using ld=long double;
const int N=45;
const int M=1565;
const ld eps=1e-12;
const ld pi=acosl(-1);
const int dx[4]={1,0,-1,0};
const int dy[4]={0,1,0,-1};
struct point
{
    ld x,y;
    point operator+(point b){return {x+b.x,y+b.y};}
    point operator-(point b){return {x-b.x,y-b.y};}
    point operator*(ld k){return {x*k,y*k};}
}v[M];
struct circle
{
    point p;
    ld r;
}a[N];
int n,cnt;
ld ans,lx,rx,ly,ry;
ld norm(point p){return p.x*p.x+p.y*p.y;}
void add(point p)
{
    for(int i=1;i<=n;i++)if(norm(p-a[i].p)<(a[i].r-eps)*(a[i].r-eps))return;
    for(int i=1;i<=cnt;i++)if(norm(p-v[i])<eps*eps)return;
    cnt++;
    v[cnt]=p;
}
bool check(point p,ld r)
{
    if(p.x-r<lx-eps||p.x+r>rx+eps||p.y-r<ly-eps||p.y+r>ry+eps)return 0;
    for(int i=0;i<4;i++)
    {
        point q={p.x+r*dx[i],p.y+r*dy[i]};
        bool flag=0;
        for(int j=1;j<=n;j++)if(norm(q-a[j].p)<=(a[j].r+eps)*(a[j].r+eps)){flag=1;break;}
        if(!flag)return 0;
    }
    pair<ld,ld> seg[2*N];
    int cnt=0;
    for(int i=1;i<=n;i++)
    {
        point q=a[i].p-p;
        ld d=sqrt(norm(q));
        if(d+r<=a[i].r)return 1;
        if(d>=r+a[i].r||d<=abs(r-a[i].r))continue;
        ld ang=atan2(q.y,q.x);
        if(ang<0)ang+=2*pi;
        ld len=acos(clamp((d*d+r*r-a[i].r*a[i].r)/(2*d*r),(ld)-1,(ld)1));
        ld l=ang-len,r=ang+len;
        if(l<0){seg[cnt++]={l+2*pi,2*pi};l=0;}
        if(r>2*pi){seg[cnt++]={0,r-2*pi};r=2*pi;}
        seg[cnt++]={l,r};
    }
    sort(seg,seg+cnt);
    ld cur=0;
    for(int i=0;i<cnt;i++)
    {
        if(seg[i].first>cur+eps)return 0;
        cur=max(cur,seg[i].second);
    }
    return cur>=2*pi-eps;
}
int main()
{
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    cin>>n;
    lx=ly=1e9;
    rx=ry=-1e9;
    for(int i=1;i<=n;i++)
    {
        cin>>a[i].p.x>>a[i].p.y>>a[i].r;
        ans=max(ans,a[i].r);
        lx=min(lx,a[i].p.x-a[i].r);
        rx=max(rx,a[i].p.x+a[i].r);
        ly=min(ly,a[i].p.y-a[i].r);
        ry=max(ry,a[i].p.y+a[i].r);
    }
    for(int i=1;i<=n;i++)
    {
        for(int j=i+1;j<=n;j++)
        {
            point q=a[j].p-a[i].p;
            ld d=sqrt(norm(q));
            if(d>a[i].r+a[j].r||d<abs(a[i].r-a[j].r))continue;
            ld x=(a[i].r*a[i].r-a[j].r*a[j].r+d*d)/(2*d);
            ld y=sqrt(max((ld)0,a[i].r*a[i].r-x*x));
            point p=a[i].p+q*(x/d),o={-q.y*y/d,q.x*y/d};
            add(p+o);
            add(p-o);
        }
    }
    for(int i=1;i<=cnt;i++)
    {
        for(int j=i+1;j<=cnt;j++)
        {
            for(int k=j+1;k<=cnt;k++)
            {
                point q=v[j]-v[i],o=v[k]-v[i];
                ld d=2*(q.x*o.y-q.y*o.x);
                if(d==0)continue;
                point p=v[i]+point{(norm(q)*o.y-norm(o)*q.y)/d,(q.x*norm(o)-o.x*norm(q))/d};
                ld r=sqrt(norm(p-v[i]));
                if(r>ans&&check(p,r))ans=r;
            }
        }
    }
    cout<<fixed<<setprecision(10)<<ans<<'\n';
    return 0;
}