Closest Pair （平面最近点对）

思路讲解

https://usaco.guide/plat/sweep-line?lang=cpp

主要参考了上面的题解。

这个分治法也比较玄学，就是先按照 x 排序，然后开始分治，从一对开始，逐渐变成两对，四对，八对，。。。

那么我们怎么样把两边的结果合并起来呢？这个就是分治的难点。

就是两边得到的答案，我们选一个较小值，然后把两边至少 x 满足这个值的点加入进来。

即

ret=dfs(l,mid);
ret=min(ret,dfs(mid+1,r));
Real x=A[mid+1].x;
vector<Point> wp;
ROF(i,mid,l){
    Real d=abs(A[i].x-x);
    if(d*d<ret.d){
        wp.EB(A[i]);
    }else{
        break;
    }
}
FOR(i,mid+1,r){
    Real d=abs(A[i].x-x);
    if(d*d<ret.d){
        wp.EB(A[i]);
    }else{
        break;
    }
}

然后得到了 wp 数组，将其按照 y 排序

Best mindis(vector<Point> &a,Best res){
    sort(all(a),[&](Point a,Point b){
        return a.y<b.y;
    });
    FOR(i,1,SZ(a)-1){  // 选一个点
		    // 向后扫描
        ROF(j,i-1,0){
            Point p=a[i],q=a[j];
            // 如果这个点的 y 已经不符合要求了，就直接开摆，后面的点y肯定更离谱
            if((p.y-q.y)*(p.y-q.y)>res.d){
                break;
            }
            Point v=p-q;
            res=min(res,{dot(v, v),p,q});
        }
    }
    return res;
}

struct Best{
    Real d;
    Point a,b;
    bool operator<(const Best &o)const{
        return d<o.d;
    }
};
Best mindis(vector<Point> &a,Best res){
    sort(all(a),[&](Point a,Point b){
        return a.y<b.y;
    });
    FOR(i,1,SZ(a)-1){
        ROF(j,i-1,0){
            Point p=a[i],q=a[j];
            if((p.y-q.y)*(p.y-q.y)>res.d){
                break;
            }
            Point v=p-q;
            res=min(res,{dot(v, v),p,q});
        }
    }
    return res;
}

Best dfs(int l,int r){
    if(l>=r) return {(ld)INF,A[l],A[l]}; 
    if(r-l==1){
        Point v=A[r]-A[l];
        return {dot(v,v),A[l],A[r]};
    }
    Best ret;
    int mid=r+l>>1;
    ret=dfs(l,mid);
    ret=min(ret,dfs(mid+1,r));
    Real x=A[mid+1].x;
    vector<Point> wp;
    ROF(i,mid,l){
        Real d=abs(A[i].x-x);
        if(d*d<ret.d){
            wp.EB(A[i]);
        }else{
            break;
        }
    }
    FOR(i,mid+1,r){
        Real d=abs(A[i].x-x);
        if(d*d<ret.d){
            wp.EB(A[i]);
        }else{
            break;
        }
    }
    return mindis(wp,ret);
}

AC代码

https://vjudge.net/solution/63334518

源代码

// by Gospel_rock
#include <bits/stdc++.h>
#include <vector>
#define all(vec) vec.begin(),vec.end()
#define PB push_back
#define EB emplace_back
#define EM emplace
#define FI first
#define SE second
#define pct __builtin_popcountll
#define ctz __builtin_ctzll
#define SZ(a) ((long long) a.size())
#define FOR(i, a, b) for (int i = (a); i <= (b); ++i)
#define ROF(i, a, b) for (int i = (a); i >= (b); --i)
#define debug(var) cerr << #var <<":"<<var<<"\n";
#define cend cerr<<"\n-----------\n"
#define fsp(x) fixed<<setprecision(x)
#define minn(vec) *min_element(vec.begin(),vec.end())
#define maxx(vec) *max_element(vec.begin(),vec.end())
using namespace std;
#ifdef LOCAL
template<typename T> void _print(T x) { cerr << x; }
// 可以打印多层数组，但是格式不是很好
template<typename T> void _print(vector<T> v) { cerr << "[ "; for (T i : v) _print(i), cerr << " "; cerr << "]"; }
template<typename T> void _print(set<T> s) { cerr << "{ "; for (T i : s) _print(i), cerr << " "; cerr << "}"; }
template<typename T, typename U> void _print(pair<T, U> p) { cerr << "{"; _print(p.first); cerr << ", "; _print(p.second); cerr << "}"; }
template<typename T, typename U> void _print(map<T, U> m) { cerr << "{ "; for (auto &p : m) _print(p), cerr << " "; cerr << "}"; }
// 打印两层数组，格式好一些
template<typename T> void _print(vector<vector<T>> v){cerr<<"{\n";for (auto i : v) _print(i), cerr << "\n";cerr<<" }";}
#endif

using  ll = long long;using ull = unsigned long long;
using ld=double;using LD=long double;
using i128 = __int128;
using pdd = pair<ld,ld>;using plb = pair<ll,bool>;
using pll = pair<ll,ll>;using vll=vector<ll>;using vb=vector<bool>;
using tll = tuple<ll,ll,ll>;using vii=vector<int>;
using vpll = vector<pll>;using vtll = vector<tll>;
using vdb = vector<ld>;using vc=vector<char>;
using arr3 = array<ll,3> ;using arr2 = array<ll,2>;
using CD=complex<double>;
static constexpr ll MAXN=(ll)1e5+10,INF=(ll)1e17+3; // 1e17+3是素数 1e18+9是素数
static constexpr ll mod=(ll)1e9+7; // 998244353;
static constexpr double eps=1e-8;const double pi=acos(-1.0);
ll N,M,Q,X,K,T,lT;

inline void __();
signed main()
{
    ios::sync_with_stdio(false);
    cin.tie(0);cout.tie(0);
    
    while(cin>>N)
        __();
    return 0;
}
/*

*/
using Real = ld;
int sgn(Real x) {
    if (-eps < x && x < eps) return 0;
    return x < 0 ? -1 : 1;
}
struct Point{
    Real x,y;
    Point(){}
    Point(Real x,Real y):x(x),y(y){}
    //                      & 是指只对左值起效果，不能对Point(1,2)这样的临时值
    // friend Point& operator+=(Point &a,Point b){
    //     a.x+=b.x;
    //     a.y+=b.y;
    //     return a;
    // }
    friend Point operator +(Point a,Point b){
        return Point(a.x+b.x,a.y+b.y);
    }
    friend Point operator -(Point a,Point b){
        return Point(a.x-b.x,a.y-b.y);
    }
    friend bool operator <(Point a,Point b) {
        return sgn(a.x-b.x)?a.x<b.x:a.y<b.y; 
    }
    friend bool operator ==(Point a,Point b){
        if(sgn(a.x-b.x)) return false;
        if(sgn(a.y-b.y)) return false;
        return true;
    }
    friend Point operator *(Point a,Real k){
        return Point(a.x*k,a.y*k);
    }
    friend Point operator /(Point a,Real k){
        return Point(a.x/k,a.y/k);
    }
    friend bool operator !=(Point a,Point b){
        return !(a==b);
    }
    friend istream& operator >>(istream& is,Point &e){
        Real a,b;is>>a>>b;
        e=Point(a,b);
        return is;
    }
    friend ostream &operator<<(ostream &os, Point p) {
        return os << "(" << p.x << ", " << p.y << ")";
    }
};

Real cross(Point a,Point b){
    return a.x*b.y-a.y*b.x;
}
Real cross(Point p1, Point p2, Point p0) { // 叉乘 (p1 - p0) x (p2 - p0);
    return cross(p1 - p0, p2 - p0);
}
Real dot(Point a,Point b){
    return a.x*b.x+a.y*b.y;
}
// 采用两点式
struct Line{
    Point p,q;
    Line(){}
    Line(Point a,Point b):p(a),q(b) {}
    friend bool operator<(Line a,Line b){
        if(a.p!=b.p) return a.p<b.p;
        return a.q<b.q; 
    }
    friend bool operator==(Line a,Line b){
        if(a.p!=b.p) return false;
        if(a.q!=b.q) return false;
        return true;
    } 
    friend istream& operator>>(istream& is, Line& e) {
        Real a,b,c,d;
        is >> a >> b >> c >> d;
        e = Line(Point(a,b),Point(c,d));
        return is;
    }
    friend ostream&operator<<(ostream &os, Line l) {
        return os << "<" << l.p << ", " << l.q << ">";
    }
};

// p+tv=q+sw
// (p-q) x w+tv x w = 0
// t=-(p-q)x w/(vxw)
// 传入的是方向式  拿到的是两直线交点
Point getintersect(const Line &a,const Line &b){
    Point p=a.p, v=a.q-a.p,q=b.p,w=b.q-b.p;
    Point u=p-q;
    ld t=cross(w,u)/cross(v,w);
    return p+v*t;
}
bool seginterline(Line seg,Line e){
    if(sgn(cross(e.p-seg.p, e.q-seg.p))*sgn(cross(e.p-seg.q, e.q-seg.q))>0){
        return false;
    }
    return true;
}
bool onseg(Point p,Line s){
    return sgn(cross(p-s.p, p-s.q))==0 && sgn(dot(p-s.p,p-s.q))<=0;
}
// AC https://vjudge.net/problem/OpenJ_Bailian-1410
bool seginterseg(Line l1, Line l2) {
    Point s1 = l1.p, e1 = l1.q;
    Point s2 = l2.p, e2 = l2.q;
    // s2 相对于 s1->e1 其在哪一侧？
    int o1 = sgn(cross(e1, s2, s1));
    int o2 = sgn(cross(e1, e2, s1));
    int o3 = sgn(cross(e2, s1, s2));
    int o4 = sgn(cross(e2, e1, s2));
    // 防止把共线但不重叠的两线段判定为相交
    if (o1 == 0 && onseg(s2, l1)) return true;
    if (o2 == 0 && onseg(e2, l1)) return true;
    if (o3 == 0 && onseg(s1, l2)) return true;
    if (o4 == 0 && onseg(e1, l2)) return true;
    // 需要都在不同侧
    return (o1 * o2 < 0) && (o3 * o4 < 0);
}
Point A[MAXN];
struct Best{
    Real d;
    Point a,b;
    bool operator<(const Best &o)const{
        return d<o.d;
    }
};
Best mindis(vector<Point> &a,Best res){
    sort(all(a),[&](Point a,Point b){
        return a.y<b.y;
    });
    FOR(i,1,SZ(a)-1){
        ROF(j,i-1,0){
            Point p=a[i],q=a[j];
            if((p.y-q.y)*(p.y-q.y)>res.d){
                break;
            }
            Point v=p-q;
            res=min(res,{dot(v, v),p,q});
        }
    }
    return res;
}

Best dfs(int l,int r){
    if(l>=r) return {(ld)INF,A[l],A[l]}; 
    if(r-l==1){
        Point v=A[r]-A[l];
        return {dot(v,v),A[l],A[r]};
    }
    Best ret;
    int mid=r+l>>1;
    ret=dfs(l,mid);
    ret=min(ret,dfs(mid+1,r));
    Real x=A[mid+1].x;
    vector<Point> wp;
    ROF(i,mid,l){
        Real d=abs(A[i].x-x);
        if(d*d<ret.d){
            wp.EB(A[i]);
        }else{
            break;
        }
    }
    FOR(i,mid+1,r){
        Real d=abs(A[i].x-x);
        if(d*d<ret.d){
            wp.EB(A[i]);
        }else{
            break;
        }
    }
    return mindis(wp,ret);
}
inline void __(){
    if(N==0) return;
    FOR(i,1,N){
        cin>>A[i];
    }
    sort(A+1,A+1+N);
    Best ans=dfs(1, N);
    cout<<fsp(8)<<ans.a.x<<" "<<ans.a.y<<" "<<ans.b.x<<" "<<ans.b.y<<"\n";
}
/*
AC
https://open.kattis.com/submissions/17958967
*/

心路历程（WA，TLE，MLE……）