三維最近點對?

//   __________________
//  | ________________ |
//  ||          ____  ||
//  ||   /\    |      ||
//  ||  /__\   |      ||
//  || /    \  |____  ||
//  ||________________||
//  |__________________|
//  \###################\
//   \###################\
//    \        ____       \
//     \_______\___\_______\
// An AC a day keeps the doctor away.

#pragma g++ optimize("Ofast")
#pragma loop_opt(on)
#include <bits/extc++.h>
#ifdef local
#define debug(x) (cerr<<#x<<" = "<<(x)<<'\n')
#else
#include <bits/stdc++.h>
#define debug(x) ((void)0)
#endif // local
#define all(v) begin(v),end(v)
#define siz(v) (ll(v.size()))
#define get_pos(v,x) (lower_bound(all(v),x)-begin(v))
#define sort_uni(v) sort(begin(v),end(v)),v.erase(unique(begin(v),end(v)),end(v))
#define pb emplace_back
#define ff first
#define ss second
#define mem(v,x) memset(v,x,sizeof v)

using namespace std;
using namespace __gnu_pbds;
typedef int64_t ll;
typedef long double ld;
typedef pair<ll,ll> pll;
typedef pair<ld,ld> pld;
template <typename T> using max_heap = std::priority_queue<T,vector<T>,less<T> >;
template <typename T> using min_heap = std::priority_queue<T,vector<T>,greater<T> >;
template <typename T> using rbt = tree<T,null_type,less<T>,rb_tree_tag,tree_order_statistics_node_update>;
constexpr ld PI = acos(-1), eps = 1e-5;
constexpr ll N = 125, INF = 1e18, MOD = 1000000007, K = 26, inf = 1e9;
constexpr inline ll cdiv(ll x, ll m) { return x/m + ((x<0 ^ m>0) && (x%m)); } // ceiling divide
constexpr inline ll modpow(ll e,ll p,ll m=MOD) {ll r=1; for(e%=m;p;p>>=1,e=e*e%m) if(p&1) r=r*e%m; return r;}

struct pt {
    int x,y,z,i;
    friend bool operator<(const pt &a, const pt &b) {
        return a.z!=b.z ? a.z<b.z : a.i<b.i;
    }
};
signed main() {
    ios_base::sync_with_stdio(0), cin.tie(0);
    int n;
    cin >> n;
    vector<pt> p(n);
    int cnt=0;
    for(auto &[x, y, z, i]: p) cin >> x >> y >> z, i=++cnt;
    sort(p.begin(), p.end(), [](pt &a, pt &b){
        return a.x!=b.x ? a.x<b.x : a.i < b.i;
    });
    auto mag2 = [](pt a, pt b) -> long long {
        auto [xa, ya, za, ia] = a;
        auto [xb, yb, zb, ib] = b;
        return 1LL*(xa-xb)*(xa-xb) + 1LL*(ya-yb)*(ya-yb) + 1LL*(za-zb)*(za-zb);
    };
    pt x=p[0], y=p[1];
    function<void(int,int)> solve = [&](int l, int r) {
        if(l+1 == r) return;
        int m = l+(r-l)/2;
        solve(l, m), solve(m, r);
        vector<pt> tmp(r-l);
        merge(p.begin()+l, p.begin()+m, p.begin()+m, p.begin()+r, tmp.begin(),
            [](pt &a, pt &b){
                return a.y!=b.y ? a.y<b.y : a.i<b.i;
            });
        set<pt> s;
        long double best = sqrt(mag2(x,y));
        for(int i = 0, j = 0; i < r-l; i++) {
            //if(tmp[i].z)
            auto [it, succ] = s.insert(tmp[i]);
            while(tmp[j].y < tmp[i].y - best) s.erase(tmp[j++]);
            for(auto lit = it; ++lit != s.end(); ) {
                if(lit->z > tmp[i].z + best) break;
                if(mag2(*lit, tmp[i]) < mag2(x,y))
                    x=*lit, y=tmp[i];
            }
            for(auto rit = it; rit-- != s.begin(); ) {
                if(rit->z < tmp[i].z - best) break;
                if(mag2(*rit, tmp[i]) < mag2(x,y))
                    x=*rit, y=tmp[i];
            }
        }
        for(int i = l; i < r; i++) p[i] = tmp[i-l];
    };
    solve(0, n);
    int a = x.i, b = y.i;
    if(a > b) swap(a,b);
    cout << "WARNING: galaxy" << a << " and galaxy" << b << " in ";
    cout << fixed << setprecision(10) << sqrt(mag2(x,y)) << " Uu\n";
}