comp-geometry

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

View the Project on GitHub Luzhiled/comp-geometry

:heavy_check_mark: src/real-geometry/point-cloud/closest-pair.hpp

Depends on

Verified with

Code

#pragma once

#include "src/real-geometry/class/point.hpp"
#include "src/real-geometry/compare/compare-x.hpp"
#include "src/real-geometry/compare/compare-y.hpp"
#include "src/real-geometry/common/size-alias.hpp"
#include "src/real-geometry/utility/sign.hpp"

#include <algorithm>
#include <cmath>
#include <complex>
#include <iterator>
#include <limits>
#include <utility>

namespace geometry::internal {

  template< typename R >
  using closest_pair_result_t = std::pair< R, std::pair<point<R>, point<R> > >;

  // WARNING: pts are sorted by y after calling this function
  template< typename R >
  closest_pair_result_t<R> impl_closest_pair(points<R> &pts, usize l, usize r) {
    constexpr R inf = std::numeric_limits< R >::infinity();
    using result_t = closest_pair_result_t<R>;

    auto comp = [](const result_t &lhs, const result_t &rhs) {
      return lhs.first < rhs.first;
    };

    if (r - l <= 1) {
      return {inf, {}};
    }

    usize m = (l + r) / 2;
    R x = pts[m].x();
    result_t result = std::min(impl_closest_pair(pts, l, m), impl_closest_pair(pts, m, r), comp);

    auto f = pts.begin();
    std::inplace_merge(f + l, f + m, f + r, compare_y<R>);

    points<R> ps;
    for (usize i = l; i < r; i++) {
      if (sign(std::abs(pts[i].x() - x) - result.first) >= 0) continue;

      for (usize j = 0; j < ps.size(); j++) {
        R dy = pts[i].y() - (*std::next(ps.rbegin(), j)).y();
        if (sign(dy - result.first) >= 0) break;

        auto &u = pts[i];
        auto &v = *std::next(ps.rbegin(), j);
        result = std::min(result, {std::abs(u - v), std::make_pair(u, v)}, comp);
      }

      ps.emplace_back(pts[i]);
    }

    return result;
  }

}

namespace geometry {

  template< typename R >
  internal::closest_pair_result_t<R> closest_pair(points<R> pts) {
    constexpr R inf = std::numeric_limits< R >::infinity();
    if (pts.size() <= 1) {
      return {inf, {}};
    }

    std::sort(pts.begin(), pts.end(), compare_x<R>);

    return internal::impl_closest_pair(pts, 0, pts.size());
  }

}
#line 2 "src/real-geometry/point-cloud/closest-pair.hpp"

#line 2 "src/real-geometry/class/point.hpp"

#line 2 "src/real-geometry/class/vector.hpp"

#include <complex>
#include <iostream>

namespace geometry {

  template< typename R >
  class vec2d : public std::complex< R > {
    using complex = std::complex< R >;

   public:
    using complex::complex;

    vec2d(const complex &c): complex::complex(c) {}

    const R x() const { return this->real(); }
    const R y() const { return this->imag(); }

    friend vec2d operator*(const vec2d &v, const R &k) {
      return vec2d(v.x() * k, v.y() * k);
    }

    friend vec2d operator*(const R &k, const vec2d &v) {
      return vec2d(v.x() * k, v.y() * k);
    }

    friend std::istream &operator>>(std::istream &is, vec2d &v) {
      R x, y;
      is >> x >> y;
      v = vec2d(x, y);
      return is;
    }
 
  };

}
#line 4 "src/real-geometry/class/point.hpp"

#include <vector>

namespace geometry {

  template< typename R >
  using point = vec2d<R>;

  template< typename R >
  using points = std::vector< point< R > >;

}
#line 2 "src/real-geometry/compare/compare-x.hpp"

#line 2 "src/real-geometry/utility/equals/real-number.hpp"

#line 2 "src/real-geometry/utility/sign.hpp"

#line 2 "src/real-geometry/common/const/eps.hpp"

#line 2 "src/real-geometry/common/float-alias.hpp"

namespace geometry {

  using f80 = long double;
  using f64 = double;

}
#line 4 "src/real-geometry/common/const/eps.hpp"

namespace geometry {

  inline static f80 &eps() {
    static f80 EPS = 1e-10;
    return EPS;
  }

  void set_eps(f80 EPS) {
    eps() = EPS;
  }

}
#line 2 "src/real-geometry/numbers/sign.hpp"

#line 2 "src/real-geometry/common/int-alias.hpp"

namespace geometry {

  using i32 = int;
  using i64 = long long;

}
#line 4 "src/real-geometry/numbers/sign.hpp"

namespace geometry::number::sign {

  constexpr i32 PLUS  = +1;
  constexpr i32 ZERO  =  0;
  constexpr i32 MINUS = -1;

}
#line 5 "src/real-geometry/utility/sign.hpp"

namespace geometry {

  using namespace geometry::number::sign;

  template< typename R >
  inline int sign(R r) {
    if (r < -eps()) return MINUS;
    if (r > +eps()) return PLUS;
    return ZERO;
  }

}
#line 4 "src/real-geometry/utility/equals/real-number.hpp"

namespace geometry {

  template< typename R >
  bool equals(R a, R b) {
    return sign(a - b) == 0;
  }

}
#line 5 "src/real-geometry/compare/compare-x.hpp"

namespace geometry {

  template< typename R >
  bool compare_x(const point<R> &a, const point<R> &b) {
    return not equals(a.x(), b.x()) ? a.x() < b.x() : a.y() < b.y();
  }

}
#line 2 "src/real-geometry/compare/compare-y.hpp"

#line 5 "src/real-geometry/compare/compare-y.hpp"

namespace geometry {

  template< typename R >
  bool compare_y(const point<R> &a, const point<R> &b) {
    return not equals(a.y(), b.y()) ? a.y() < b.y() : a.x() < b.x();
  }

}
#line 2 "src/real-geometry/common/size-alias.hpp"

#include <cstddef>

namespace geometry {

  using isize = std::ptrdiff_t;
  using usize = std::size_t;

}
#line 8 "src/real-geometry/point-cloud/closest-pair.hpp"

#include <algorithm>
#include <cmath>
#line 12 "src/real-geometry/point-cloud/closest-pair.hpp"
#include <iterator>
#include <limits>
#include <utility>

namespace geometry::internal {

  template< typename R >
  using closest_pair_result_t = std::pair< R, std::pair<point<R>, point<R> > >;

  // WARNING: pts are sorted by y after calling this function
  template< typename R >
  closest_pair_result_t<R> impl_closest_pair(points<R> &pts, usize l, usize r) {
    constexpr R inf = std::numeric_limits< R >::infinity();
    using result_t = closest_pair_result_t<R>;

    auto comp = [](const result_t &lhs, const result_t &rhs) {
      return lhs.first < rhs.first;
    };

    if (r - l <= 1) {
      return {inf, {}};
    }

    usize m = (l + r) / 2;
    R x = pts[m].x();
    result_t result = std::min(impl_closest_pair(pts, l, m), impl_closest_pair(pts, m, r), comp);

    auto f = pts.begin();
    std::inplace_merge(f + l, f + m, f + r, compare_y<R>);

    points<R> ps;
    for (usize i = l; i < r; i++) {
      if (sign(std::abs(pts[i].x() - x) - result.first) >= 0) continue;

      for (usize j = 0; j < ps.size(); j++) {
        R dy = pts[i].y() - (*std::next(ps.rbegin(), j)).y();
        if (sign(dy - result.first) >= 0) break;

        auto &u = pts[i];
        auto &v = *std::next(ps.rbegin(), j);
        result = std::min(result, {std::abs(u - v), std::make_pair(u, v)}, comp);
      }

      ps.emplace_back(pts[i]);
    }

    return result;
  }

}

namespace geometry {

  template< typename R >
  internal::closest_pair_result_t<R> closest_pair(points<R> pts) {
    constexpr R inf = std::numeric_limits< R >::infinity();
    if (pts.size() <= 1) {
      return {inf, {}};
    }

    std::sort(pts.begin(), pts.end(), compare_x<R>);

    return internal::impl_closest_pair(pts, 0, pts.size());
  }

}
Back to top page