library/geometry/2dPointSet.hpp

• View this file on GitHub
• Last update: 2024-03-15 16:42:38+06:00
• Include: #include "library/geometry/2dPointSet.hpp"

Depends on

• library/geometry/2dGeometryTemplate.hpp

Verified with

• library/test/aoj/CGL_5_A-Closest_Pair.test.cpp

Code

#pragma once
#include "2dGeometryTemplate.hpp"

Real closest_pair(Polygon ps) {
    sort(ps.begin(), ps.end());
    Polygon a(ps.size());

    function<Real(int, int)> rec = [&](int left, int right) -> Real {
        if (right - left <= 1) return 1e18;
        int mid = (left + right) / 2;
        Real x = ps[mid].real();
        Real ret = min(rec(left, mid), rec(mid, right));
        inplace_merge(ps.begin() + left, ps.begin() + mid, ps.begin() + right,
                      [&](const Point& a, const Point& b) { return a.imag() < b.imag(); });
        int pos = 0;
        for (int i = left; i < right; i++) {
            if (fabs((ps[i].real()) - x) >= ret) continue;
            for (int j = 0; j < pos; j++) {
                auto tar = ps[i] - a[pos - j - 1];
                if (tar.imag() >= ret) break;
                ret = min(ret, abs(tar));
            }
            a[pos++] = ps[i];
        }
        return ret;
    };
    return rec(0, (int)ps.size());
}

#line 2 "library/geometry/2dGeometryTemplate.hpp"

using Real = double;
using Point = complex<Real>;
using Polygon = vector<Point>;
const Real EPS = 1e-8, PI = acos(-1);

Point operator*(const Point& p, const Real& d) {
    return Point(p.real() * d, p.imag() * d);
}

Point operator/(const Point& p, const Real& d) {
    return Point(p.real() / d, p.imag() / d);
}

istream& operator>>(istream& is, Point& p) {
    Real a, b;
    is >> a >> b;
    p = Point(a, b);
    return is;
}

ostream& operator<<(ostream& os, const Point& p) {
    return os << fixed << setprecision(20) << p.real() << " " << p.imag();
}

int sign(const Real& r) {
    if (r <= -EPS) return -1;
    if (r >= +EPS) return +1;
    return 0;
}

bool equals(const Real& a, const Real& b) {
    return sign(a - b) == 0;
}

namespace std {
bool operator<(const Point& a, const Point& b) {
    if (equals(a.real(), b.real())) return a.imag() < b.imag();
    return a.real() < b.real();
}
}  // namespace std

Real dot(const Point& a, const Point& b) {
    return (conj(a) * b).real();
}

Real cross(const Point& a, const Point& b) {
    return (conj(a) * b).imag();
}

struct Line {
    Point a, b;
    Line() = default;
    Line(Point a, Point b) : a(a), b(b) {}
};
using Segment = Line;
#line 3 "library/geometry/2dPointSet.hpp"

Real closest_pair(Polygon ps) {
    sort(ps.begin(), ps.end());
    Polygon a(ps.size());

    function<Real(int, int)> rec = [&](int left, int right) -> Real {
        if (right - left <= 1) return 1e18;
        int mid = (left + right) / 2;
        Real x = ps[mid].real();
        Real ret = min(rec(left, mid), rec(mid, right));
        inplace_merge(ps.begin() + left, ps.begin() + mid, ps.begin() + right,
                      [&](const Point& a, const Point& b) { return a.imag() < b.imag(); });
        int pos = 0;
        for (int i = left; i < right; i++) {
            if (fabs((ps[i].real()) - x) >= ret) continue;
            for (int j = 0; j < pos; j++) {
                auto tar = ps[i] - a[pos - j - 1];
                if (tar.imag() >= ret) break;
                ret = min(ret, abs(tar));
            }
            a[pos++] = ps[i];
        }
        return ret;
    };
    return rec(0, (int)ps.size());
}

Back to top page