F. Grand Finale: Circles 详细题解

问题重述

给定 $n$ 个圆，第 $i$ 个圆由圆心 $(x_i, y_i)$ 和半径 $r_i$ 给出。
求一个圆 $C$（圆心 $(x, y)$，半径 $r$），使得 $C$ 被所有给定的圆包含，即对每个 $i$， $\sqrt{(x_i - x)^2 + (y_i - y)^2} + r \le r_i,$ 且 $r$ 最大。输出满足条件的 $C$，允许 $10^{-7}$ 的误差。

问题转化

设 $P = (x, y)$ 为待求圆心，$r$ 为半径。条件可改写为 $\|P - C_i\| \le r_i - r, \quad \forall i,$ 其中 $C_i = (x_i, y_i)$，$\|\cdot\|$ 为欧氏距离。

定义 $r_i' = r_i - r$，则问题变为：寻找最大的 $r$，使得存在点 $P$ 满足 $\|P - C_i\| \le r_i'$ 对所有 $i$ 成立。
当 $r$ 增大时，每个 $r_i'$ 减小，这些圆的交集逐渐缩小。最大可行的 $r$ 就是这些圆仍然有公共交点的最大值。
换句话说，我们要在所有给定圆内部找一个半径最大的内切圆。

算法思想

这是一个典型的最大空圆问题，但约束不是点而是圆。可以转化为最小外接圆的对偶问题，使用随机增量法（Welzl 算法）求解。

核心观察：
最大内含圆 $C$ 要么由单个圆决定（即 $C$ 与某个给定圆内切，且圆心为该圆圆心，半径为该圆半径），要么由两个圆决定（$C$ 同时与两个给定圆内切，圆心位于两圆圆心连线的中垂线上，半径由几何关系确定），要么由三个圆决定（$C$ 与三个给定圆内切，可通过解方程组得到）。
因此，最优解一定由至多三个给定圆“支撑”，这与最小圆覆盖（包含所有点的最小圆）的性质类似。

Welzl 算法（随机增量）可以高效地找到这样的圆。算法如下：

1. 随机打乱所有圆的顺序。
2. 递归函数 rec(S, R)：
   • S 是尚未处理的圆的下标范围（或向量），R 是当前已经选中的支撑圆集合（大小不超过 3）。
   • 初始 R 为空，S 包含所有圆。
   • 若 R 的大小为 3，则直接计算由这三个圆确定的最大内含圆（切于三者）。
   • 若 S 为空，则根据 R 的大小（0/1/2）分别返回对应的圆（无穷大半径、单个圆、切于两个圆的圆）。
   • 否则，先递归处理 S 去掉最后一个圆，得到当前圆 C。
   • 如果 C 已经包含在 S 的最后一个圆中（即 C 完全位于该圆内部），则直接返回 C。
   • 否则，将该圆加入 R 中，重新递归处理 S 去掉最后一个圆。
   • 返回结果。

由于每一步 R 的大小最多增加 1，递归深度为 $O(n)$。随机打乱保证了期望复杂度为 $O(n)$。

几何基础函数

代码中定义了若干辅助函数：

• midp(a, b)：两点中点。
• addi(a, b)，subt(a, b)，mult(a, k)：向量运算。
• abs(p)，dist(a, b)：长度和距离。
• basis0()：返回半径无穷大的圆（代表无约束）。
• basis1(a)：返回圆 a 本身。
• basis2(a, b)：返回同时内切于圆 a 和 b 的最大圆。
  设 $A, B$ 为圆心，$r_a, r_b$ 为半径。作向量 $\overrightarrow{AB}$，分别在 $A$ 和 $B$ 的方向上延长 $r_a$ 和 $r_b$ 得到点 $A', B'$，则所求圆心为 $A'B'$ 的中点，半径为 $\frac{1}{2}|A'B'|$。
• basis3(a, b, c)：返回同时内切于三个圆的最大圆。
  通过解方程组得到，代码中的公式是解析解。
• viol(c, a)：判断圆 c 是否被圆 a 包含（即是否 c.r + dist(c.p, a.p) > a.r）。如果超出，则 c 不满足被 a 包含的条件。

随机增量法实现

标程中的 solve(v) 函数：

• 随机种子固定为 1999999，但为了随机化，后续会打乱顺序。
• basis 存储当前选中的支撑圆（最多 3 个）。
• trivial() 根据 basis 的大小返回对应的圆（0 → 无穷大，1 → 那个圆，2 → 两圆内切圆，3 → 三圆内切圆）。
• rec 是递归函数，参数 n 表示处理前 n 个圆（下标 0..n-1）。
  • 递归出口：n == 0 或 basis.size() == 3，返回 trivial()。
  • 为了防止极端情况，加入计数器 counter，当递归次数超过 $30|V|$ 时返回 NaN，触发重新打乱。
  • 先递归处理前 n-1 个圆，得到 c。
  • 取第 n-1 个圆 p = v[n-1]，如果 c 不被 p 包含（viol(c, p) 为真）且 c.r 不是 NaN，则说明 c 需要更新：将 p 加入 basis，递归求解前 n-1 个圆，然后弹出 p。
  • 返回结果。
• 第一次调用 rec(rec, v.size()) 得到候选圆 c，如果 c.r 是 NaN（表示陷入死循环），则重新打乱整个序列，重复直到获得有效结果。

最后输出圆心坐标和半径。

复杂度分析

• 每次递归调用可能增加 basis 的大小，但最多增加到 3。
• 随机化保证了每次圆被加入 basis 的概率很低，期望总递归次数为 $O(n)$。
• 实际实现中加上计数器限制，避免最坏情况。
• 总体复杂度 $O(n)$，在 $n \le 10^5$ 时运行很快。

关于精度

算法使用 long double（80 位扩展精度），能够满足 $10^{-7}$ 的误差要求。
输出时使用 setprecision(16) 输出足够多的小数位。

标程代码

#include<bits/stdc++.h>
using namespace std;
using ll=long long;
 
using lf=long double;
using pt=pair<lf,lf>;
lf& real(pt& p){return p.first;}
lf& imag(pt& p){return p.second;}
pt midp(pt a,pt b){return pt{(real(a)+real(b))/2,(imag(a)+imag(b))/2};}
pt addi(pt a,pt b){return pt{(real(a)+real(b)),(imag(a)+imag(b))};}
pt subt(pt a,pt b){return pt{(real(a)-real(b)),(imag(a)-imag(b))};}
pt mult(pt a,lf b){return pt{real(a)*b,imag(a)*b};}
lf abs(pt a){return sqrtl(powl(real(a),2)+powl(imag(a),2));}
lf dist(pt a,pt b){return sqrtl(powl(real(a)-real(b),2)+powl(imag(a)-imag(b),2));}
 
struct circ{pt p;lf r;};
const lf inf=1e18;
 
circ basis0(){return {pt{0,0},inf};}
 
circ basis1(circ a){return a;}
 
circ basis2(circ a,circ b)
{
    pt aa=a.p,bb=b.p;
    pt ab=subt(bb,aa),ba=subt(aa,bb);
    lf aab=abs(ab),aba=abs(ba);
    real(ab)/=aab;
    imag(ab)/=aab;
    real(ba)/=aba;
    imag(ba)/=aba;
    pt ar=addi(aa,mult(ab,a.r)),br=addi(bb,mult(ba,b.r));
    
    return {midp(ar,br),dist(ar,br)/2.0L};
}
 
circ basis3(circ a,circ b,circ c)
{
    lf x1=real(a.p),y1=imag(a.p),r1=a.r;
    lf x2=real(b.p),y2=imag(b.p),r2=b.r;
    lf x3=real(c.p),y3=imag(c.p),r3=c.r;
    lf a2=x1-x2,a3=x1-x3,b2=y1-y2,b3=y1-y3,c2=r2-r1,c3=r3-r1;
    lf d1=x1*x1+y1*y1-r1*r1,d2=d1-x2*x2-y2*y2+r2*r2,d3=d1-x3*x3-y3*y3+r3*r3;
    lf ab=a3*b2-a2*b3;
    lf xa=(b2*d3-b3*d2)/(ab*2)-x1;
    lf xb=(b3*c2-b2*c3)/ab;
    lf ya=(a3*d2-a2*d3)/(ab*2)-y1;
    lf yb=(a2*c3-a3*c2)/ab;
    lf A=xb*xb+yb*yb-1;
    lf B=2*(r1+xa*xb+ya*yb);
    lf C=xa*xa+ya*ya-r1*r1;
    lf r=-(A?(B-sqrtl(B*B-4*A*C))/(2*A):C/B);
    return {pt{x1+xa+xb*r,y1+ya+yb*r},r};
}
 
bool viol(circ p,circ a)
{
    return a.r<p.r+dist(a.p,p.p);
}
 
circ solve(vector<circ>&v)
{
    lf _nan=nan("aaa");
    mt19937_64 mt(1999999);
    //shuffle(begin(v),end(v),mt);
    vector<circ>basis;
    auto trivial=[&]()
    {
        if(basis.size()==0)return basis0();
        if(basis.size()==1)return basis[0];
        if(basis.size()==2)return basis2(basis[0],basis[1]);
        return basis3(basis[0],basis[1],basis[2]);
    };
 
    int counter=0;
    
    auto rec=[&](auto rec,int n)->circ
    {
        if(n==0||basis.size()==3)return trivial();
        counter++;
        if(counter>30*(int)v.size())return {pt{0,0},_nan};
        auto c=rec(rec,n-1);
        auto p=v[n-1];
        if(!viol(c,p)||isnan(c.r))return c;
        basis.push_back(p);
        c=rec(rec,n-1);
        basis.pop_back();
        return c;
    };
    auto c=rec(rec,(int)v.size());
    while(isnan(c.r))
    {
        counter=0;
        shuffle(begin(v),end(v),mt);
        c=rec(rec,(int)v.size());
    }
    return c;
}
 
int main()
{
    cin.tie(0)->sync_with_stdio(0);
    int n;cin>>n;
    vector<circ>v(n);
    for(auto&[p,r]:v)
    {
        ll x,y,rr;cin>>x>>y>>rr;
        p={x,y};r=rr;
    }
    auto ans=solve(v);
    cout<<setprecision(16)<<fixed<<real(ans.p)<<" "<<imag(ans.p)<<" "<<ans.r<<"\n";
}

总结

本题要求在所有给定圆的内部找一个最大内切圆，可以转化为求一个点 $P$ 使得它到各圆心距离加上自身半径不超过给定半径。随机增量法（Welzl 算法）通过维护至多三个支撑圆，递归求解，具有期望线性时间复杂度，实现简洁，精度满足要求。

代码中 basis2 和 basis3 的几何推导是核心，理解了两圆/三圆内切圆的构造方法，就能理解整个算法。