P16629 [GKS 2017 #E] Blackhole

题面

首先三点必共平面，直接变二维做，具体来说算出三条边长，解三角形即可。

考虑二分最小半径 \(r\)。

首先是三个圆心 \(O_1, O_2, O_3\) 均包含至少一个黑洞的情况。

此时三个点分别在以三个黑洞 \(P_1, P_2, P_3\) 为圆心，半径为 \(r\) 的圆内。

若 \(O_x,O_y\) 形成单一联通区域，则其距离小于 \(2r\)，也即 \(|\vec{O_x}-\vec{O_y}|\le 2r\)。

若 \(O_1,O_2,O_3\) 联通，则至少存在一个 \(O_i\) 与另外两个 \(O_j,O_k\) 均联通。

由 \(|\vec{O_j}-\vec{P_j}| \le r\)，\(|\vec{O_i}-\vec{O_j}|\le 2r\)，\(O_iO_j\) 联通充要条件为 \(|\vec{O_i}-\vec{P_j}|\le 3r\)。

同理 \(|\vec{O_i}-\vec{P_k}|\le 3r\)

欲使两者同时成立，则 \(\max(|\vec{O_i}-\vec{P_j}|,|\vec{O_i}-\vec{P_k}|)\le 3r\)。

记前者为 \(\max(d_1,d_2)\)，我们不难证明其极小时要么 \(d_1\) 或 \(d_2\) 取得局部最小值，要么 \(d_1=d_2\)。

也即 \(O_i\) 位于 \(P_iP_j\) 或 \(P_iP_k\) 上，或者位于 \(P_jP_k\) 的中垂线上。

二分时直接算出最小值即可判定。

另外还有一种情况是，\(O_1\) 包含了 \(P_j,P_k\)，\(O_3\) 包含 \(P_i\)，\(O_2\) 位于 \(O_1O_3\) 连线上。

此时只需把上文中 \(O_3\) 的半径换为 \(5r\)，\(|\vec{O_i}-\vec{P}|\) 最大值上界换为 \(r\) 即可。

考虑 \(i\) 的任意性，对三角形三条边轮换做 \(3\) 次即可。

#include <algorithm>
#include <iostream>
#include <random>
#include <array>
#include <cmath>

typedef double Nx;

inline Nx dist(Nx a[3], Nx b[3]) {
    Nx d[3] = { a[0] - b[0], a[1] - b[1], a[2] - b[2] };
    return std::sqrt(d[0] * d[0] + d[1] * d[1] + d[2] * d[2]);
}

const Nx pi = acos(-1);

struct vec2 {
    Nx x, y;
    inline Nx len() { return hypot(x, y); }
    inline vec2 norm() {
        Nx t = len();
        return { x/t, y/t };
    }
};

inline vec2 operator+ (const vec2& x, const vec2& y) { return { x.x + y.x, x.y + y.y }; }
inline vec2 operator- (const vec2& x, const vec2& y) { return { x.x - y.x, x.y - y.y }; }
inline vec2 operator- (const vec2& x) { return { -x.x, -x.y }; }
inline vec2 operator* (const vec2& x, Nx k) { return { x.x * k, x.y * k }; }

inline void solve() {
    static Nx v[3][3];
    std::cin >> v[0][0] >> v[0][1] >> v[0][2]
             >> v[1][0] >> v[1][1] >> v[1][2]
             >> v[2][0] >> v[2][1] >> v[2][2];

    auto work = [](Nx r, vec2 p, Nx w, Nx ratio) {
        Nx ans = 1e9;
        vec2 q{w, 0};

        auto upd = [&](vec2 g) { ans = std::min(ans, std::max(g.len(), (q - g).len())); };

        upd(p + (-p).norm() * r);
        upd(p + (q-p).norm() * r);
        if(Nx d = std::abs(p.x - w / 2); d <= r)
            upd({w / 2, std::max((Nx)0, p.y - std::sqrt(r * r - d * d))});

        return ans <= ratio * r;
    };

    Nx a = dist(v[0], v[1]), b = dist(v[1], v[2]), c = dist(v[0], v[2]);

    auto gen = [=](Nx a, Nx b, Nx c) {
        Nx A = (b * b + c * c - a * a) / (2 * b * c);
        vec2 w { b * A, b * std::sqrt(1 - A * A) };
        return [=](auto&& r) -> bool {
            return work(r, w, c, 3) || work(5 * r, w, c, .2);
        };
    };

    auto c1 = gen(a, b, c), c2 = gen(a, c, b), c3 = gen(b, c, a);
    auto check = [&](Nx r) { return c1(r) || c2(r) || c3(r); };

    Nx l = 0, r = 1e3;
    while(r - l > 1e-9) {
        Nx m = (l + r) / 2;
        if(check(m)) r = m;
        else l = m;
    }

    static int cas = 0; std::cout.precision(10);
    std::cout << "Case #" << ++cas << ": " << l << "\n";
}

int main() {
    std::ios::sync_with_stdio(0), std::cin.tie(0), std::cout.tie(0);
    int t; std::cin >> t; while(t--)
    solve();
}