cp-library

This documentation is automatically generated by online-judge-tools/verification-helper

View the Project on GitHub rniya/cp-library

:heavy_check_mark: Furthest Pair
(src/geometry/furthest_pair.hpp)

Depends on

Verified with

Code

#pragma once
#include "convex_diameter.hpp"
#include "convex_hull.hpp"

namespace geometry {

template <typename T> std::pair<int, int> furthest_pair(const std::vector<Point<T>>& P) {
    int n = P.size();
    assert(n >= 2);
    auto convex = convex_hull(P);
    n = convex.size();
    if (n == 1) return {0, 1};
    if (n == 2) return {convex[0], convex[1]};
    Polygon<T> Q;
    for (int i : convex) Q.emplace_back(P[i]);
    auto [i, j] = convex_diameter(Q);
    return {convex[i], convex[j]};
}

}  // namespace geometry
#line 2 "src/geometry/convex_diameter.hpp"
#include <algorithm>
#include <tuple>
#include <utility>
#line 2 "src/geometry/Polygon.hpp"
#include <vector>
#line 2 "src/geometry/Point.hpp"
#include <cassert>
#include <cmath>
#include <iostream>
#include <type_traits>

namespace geometry {

template <typename T> struct Point {
    static T EPS;

    static void set_eps(T eps) { EPS = eps; }

    T x, y;

    Point() {}

    Point(T x, T y) : x(x), y(y) {}

    Point operator+(const Point& p) const { return Point(x + p.x, y + p.y); }

    Point operator-(const Point& p) const { return Point(x - p.x, y - p.y); }

    Point operator*(T t) const { return Point(x * t, y * t); }

    Point operator/(T t) const { return Point(x / t, y / t); }

    bool operator==(const Point& p) const { return x == p.x and y == p.y; }

    bool operator!=(const Point& p) const { return not((*this) == p); }

    bool operator<(const Point& p) const { return x != p.x ? x < p.x : y < p.y; }

    friend std::istream& operator>>(std::istream& is, Point& p) { return is >> p.x >> p.y; }

    friend std::ostream& operator<<(std::ostream& os, const Point& p) { return os << p.x << ' ' << p.y; }

    T norm() { return std::sqrt(x * x + y * y); }

    T norm2() { return x * x + y * y; }

    T arg() { return std::atan2(y, x); }

    T dot(const Point& p) { return x * p.x + y * p.y; }

    T det(const Point& p) { return x * p.y - y * p.x; }

    Point perp() { return Point(-y, x); }

    Point unit() { return *this / norm(); }

    Point normal() { return perp().unit(); }

    Point rotate(T rad) { return Point(std::cos(rad) * x - std::sin(rad) * y, std::sin(rad) * x + std::cos(rad) * y); }
};

template <> double Point<double>::EPS = 1e-9;
template <> long double Point<long double>::EPS = 1e-12;
template <> int Point<int>::EPS = 0;
template <> long long Point<long long>::EPS = 0;

template <typename T> int sgn(T x) { return x < -Point<T>::EPS ? -1 : x > Point<T>::EPS ? 1 : 0; }

}  // namespace geometry
#line 3 "src/geometry/ccw.hpp"

namespace geometry {

enum Position { COUNTER_CLOCKWISE = +1, CLOCKWISE = -1, ONLINE_BACK = +2, ONLINE_FRONT = -2, ON_SEGMENT = 0 };

template <typename T> Position ccw(const Point<T>& a, Point<T> b, Point<T> c) {
    b = b - a;
    c = c - a;
    if (sgn(b.det(c)) == 1) return COUNTER_CLOCKWISE;
    if (sgn(b.det(c)) == -1) return CLOCKWISE;
    if (sgn(b.dot(c)) == -1) return ONLINE_BACK;
    if (b.norm2() < c.norm2()) return ONLINE_FRONT;
    return ON_SEGMENT;
}

}  // namespace geometry
#line 4 "src/geometry/Polygon.hpp"

namespace geometry {

template <typename T> struct Polygon : std::vector<Point<T>> {
    using std::vector<Point<T>>::vector;

    Polygon(int n) : std::vector<Point<T>>(n) {}

    T area2() {
        T sum = 0;
        int n = this->size();
        for (int i = 0; i < n; i++) sum += (*this)[i].det((*this)[i + 1 == n ? 0 : i + 1]);
        return sum < 0 ? -sum : sum;
    }

    T area() { return area2() / 2; }

    bool is_convex() {
        int n = this->size();
        for (int j = 0; j < n; j++) {
            int i = (j == 0 ? n - 1 : j - 1), k = (j == n - 1 ? 0 : j + 1);
            if (ccw((*this)[i], (*this)[j], (*this)[k]) == CLOCKWISE) return false;
        }
        return true;
    }
};

}  // namespace geometry
#line 6 "src/geometry/convex_diameter.hpp"

namespace geometry {

template <typename T> std::pair<int, int> convex_diameter(const Polygon<T>& convex) {
    int n = convex.size();
    auto [si, sj] = [&] {
        auto [it_min, it_max] = std::minmax_element(begin(convex), end(convex));
        return std::pair<int, int>{it_min - begin(convex), it_max - begin(convex)};
    }();
    T max_dist = -1;
    std::pair<int, int> argmax{-1, -1};
    for (int i = si, j = sj; i != sj or j != si;) {
        T d = (convex[i] - convex[j]).norm2();
        if (max_dist < d) {
            max_dist = d;
            argmax = {i, j};
        }
        int ni = (i + 1 == n ? 0 : i + 1), nj = (j + 1 == n ? 0 : j + 1);
        if ((convex[ni] - convex[i]).det(convex[nj] - convex[j]) < 0)
            i = ni;
        else
            j = nj;
    }
    return argmax;
}

}  // namespace geometry
#line 3 "src/geometry/convex_hull.hpp"
#include <numeric>
#line 6 "src/geometry/convex_hull.hpp"

namespace geometry {

template <typename T> std::vector<int> convex_hull(const std::vector<Point<T>>& P, bool inclusive = false) {
    int n = P.size();
    if (n == 0) return {};
    if (n == 1) return {0};
    if (n == 2) return (P[0] != P[1] or inclusive ? std::vector<int>{0, 1} : std::vector<int>{0});
    std::vector<int> ord(n);
    std::iota(ord.begin(), ord.end(), 0);
    std::sort(ord.begin(), ord.end(), [&](int l, int r) { return P[l] < P[r]; });
    std::vector<int> ch(n + 1, -1);
    int s = 0, t = 0;
    for (int _ = 0; _ < 2; _++, s = --t, std::reverse(ord.begin(), ord.end())) {
        for (int& i : ord) {
            for (; t >= s + 2; t--) {
                auto det = (P[ch[t - 1]] - P[ch[t - 2]]).det(P[i] - P[ch[t - 2]]);
                if (inclusive ? det >= 0 : det > 0) break;
            }
            ch[t++] = i;
        }
    }
    while (not inclusive and t >= 2 and P[ch[0]] == P[ch[t - 1]]) t--;
    return {begin(ch), begin(ch) + t};
}

}  // namespace geometry
#line 4 "src/geometry/furthest_pair.hpp"

namespace geometry {

template <typename T> std::pair<int, int> furthest_pair(const std::vector<Point<T>>& P) {
    int n = P.size();
    assert(n >= 2);
    auto convex = convex_hull(P);
    n = convex.size();
    if (n == 1) return {0, 1};
    if (n == 2) return {convex[0], convex[1]};
    Polygon<T> Q;
    for (int i : convex) Q.emplace_back(P[i]);
    auto [i, j] = convex_diameter(Q);
    return {convex[i], convex[j]};
}

}  // namespace geometry
Back to top page