minimum_enclosing_circle.hpp¶
求覆盖所有点的最小圆、边界支撑点或圆周点标记；整数判定的中间量须能用 128 位有符号整数表示，外接圆行列式为坐标差的四次式。
Complexity: Time: Expected \(O(n)\). Space: \(O(n)\) shuffled points.
AC 记录：minimum_enclosing_circle。
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: Expected O(n).
/// Space: O(n) shuffled points.
#include "noya/geometry/geometry_base.hpp"
#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <random>
#include <type_traits>
#include <vector>
namespace noya {
template <class T> struct minimum_enclosing_circle_result {
circle<long double> val;
std::array<point<T>, 3> sup{};
int cnt = 0;
};
namespace minimum_enclosing_circle_internal {
inline circle<long double> from_diameter(const point<long double> &lhs,
const point<long double> &rhs) {
point<long double> o = (lhs + rhs) / 2.0L;
return {o, distance(o, lhs)};
}
inline circle<long double> from_three(const point<long double> &lhs,
const point<long double> &rhs,
const point<long double> &z) {
point<long double> b = rhs - lhs, c = z - lhs;
long double den = 2.0L * cross(b, c);
if (std::abs(den) <= 1e-24L) {
circle<long double> res = from_diameter(lhs, rhs);
circle<long double> can = from_diameter(lhs, z);
if (can.r1 > res.r1) {
res = can;
}
can = from_diameter(rhs, z);
if (can.r1 > res.r1) {
res = can;
}
return res;
}
// Work relative to lhs to avoid subtracting large absolute squares.
long double bn = norm2(b), cn = norm2(c);
point<long double> v{(bn * c.y - cn * b.y) / den,
(b.x * cn - c.x * bn) / den};
return {lhs + v, std::sqrt(norm2(v))};
}
inline bool contains(const circle<long double> &val,
const point<long double> &can) {
long double r2 = val.r1 * val.r1;
long double tol = 1e-12L * std::max(1.0L, r2);
return norm2(can - val.o) <= r2 + tol;
}
} // namespace minimum_enclosing_circle_internal
/// @brief Return the minimum circle and one boundary support set. Randomly
/// ordering the points makes the expected number of constraint rebuilds
/// linear: an outside point must belong to the new boundary, reducing the
/// remaining problem successively to circles through one, two, then three
/// fixed points. Integral predicate intermediates must fit in signed 128 bits;
/// the circumcircle determinant has degree four in coordinate differences.
template <class T>
minimum_enclosing_circle_result<T>
minimum_enclosing_circle_with_support(const std::vector<point<T>> &in,
std::uint64_t sd = 712367821ULL) {
assert(!in.empty());
struct indexed_point {
point<long double> val;
point<T> org;
};
std::vector<indexed_point> pt;
pt.reserve(in.size());
for (const point<T> &val : in) {
pt.push_back(
{{static_cast<long double>(val.x), static_cast<long double>(val.y)},
val});
}
std::mt19937_64 rng(sd);
std::shuffle(pt.begin(), pt.end(), rng);
using namespace minimum_enclosing_circle_internal;
minimum_enclosing_circle_result<T> ans;
ans.val = {pt[0].val, 0};
ans.sup[0] = pt[0].org;
ans.cnt = 1;
auto has = [&](const indexed_point &can) {
if constexpr (!std::is_integral_v<T>) {
return contains(ans.val, can.val);
} else {
using wide = __int128;
const auto &a = ans.sup[0];
const auto &p = can.org;
if (ans.cnt == 1) {
return p == a;
}
const auto &b = ans.sup[1];
if (ans.cnt == 2) {
wide dx = wide(2) * p.x - a.x - b.x;
wide dy = wide(2) * p.y - a.y - b.y;
wide bx = wide(a.x) - b.x, by = wide(a.y) - b.y;
wide d2 = bx * bx + by * by;
return dx * dx + dy * dy <= d2;
}
const auto &c = ans.sup[2];
wide bx = wide(b.x) - a.x;
wide by = wide(b.y) - a.y;
wide cx = wide(c.x) - a.x;
wide cy = wide(c.y) - a.y;
wide px = wide(p.x) - a.x;
wide py = wide(p.y) - a.y;
wide det = (px * px + py * py) * (bx * cy - by * cx) -
(bx * bx + by * by) * (px * cy - py * cx) +
(cx * cx + cy * cy) * (px * by - py * bx);
// Compare signs without multiplying the fourth-degree determinant.
return (bx * cy - by * cx) > 0 ? det <= 0 : det >= 0;
}
};
for (int i = 1; i < int(pt.size()); i++) {
if (has(pt[i])) {
continue;
}
ans.val = {pt[i].val, 0};
ans.sup[0] = pt[i].org;
ans.cnt = 1;
for (int j = 0; j < i; j++) {
if (has(pt[j])) {
continue;
}
ans.val = from_diameter(pt[i].val, pt[j].val);
ans.sup[0] = pt[i].org;
ans.sup[1] = pt[j].org;
ans.cnt = 2;
for (int k = 0; k < j; k++) {
if (!has(pt[k])) {
long double ta = cross(pt[i].val, pt[j].val, pt[k].val);
if (std::abs(ta) <= 1e-24L) {
int lhs = i;
int rhs = j;
long double bst = norm2(pt[i].val - pt[j].val);
for (auto [l, r] : {std::pair{i, k}, std::pair{j, k}}) {
long double can = norm2(pt[l].val - pt[r].val);
if (can > bst) {
bst = can;
lhs = l;
rhs = r;
}
}
ans.val = from_diameter(pt[lhs].val, pt[rhs].val);
ans.sup[0] = pt[lhs].org;
ans.sup[1] = pt[rhs].org;
ans.cnt = 2;
} else {
ans.val = from_three(pt[i].val, pt[j].val, pt[k].val);
ans.sup[0] = pt[i].org;
ans.sup[1] = pt[j].org;
ans.sup[2] = pt[k].org;
ans.cnt = 3;
}
}
}
}
}
return ans;
}
/// @brief Return the minimum circle containing all points in expected O(n).
template <class T>
circle<long double> minimum_enclosing_circle(const std::vector<point<T>> &in,
std::uint64_t sd = 712367821ULL) {
return minimum_enclosing_circle_with_support(in, sd).val;
}
/// @brief Mark exactly which integral points lie on the minimum circle. Once
/// the randomized construction supplies two or three support points, equality
/// with their circle is tested after clearing all denominators: a doubled
/// midpoint equation for a diameter, or the integer circumcircle det
/// for three non-collinear points.
template <class T>
std::vector<bool>
minimum_enclosing_circle_boundary(const std::vector<point<T>> &in,
std::uint64_t sd = 712367821ULL) {
static_assert(std::is_integral_v<T>);
auto ans = minimum_enclosing_circle_with_support(in, sd);
std::vector<bool> res(in.size());
using wide = __int128;
const auto &a = ans.sup[0];
if (ans.cnt == 1) {
for (int i = 0; i < int(in.size()); i++) {
res[i] = in[i] == a;
}
} else if (ans.cnt == 2) {
const auto &b = ans.sup[1];
wide bx = wide(a.x) - b.x, by = wide(a.y) - b.y;
wide r4 = bx * bx + by * by;
for (int i = 0; i < int(in.size()); i++) {
wide dx = wide(2) * in[i].x - a.x - b.x;
wide dy = wide(2) * in[i].y - a.y - b.y;
res[i] = dx * dx + dy * dy == r4;
}
} else {
const auto &b = ans.sup[1];
const auto &c = ans.sup[2];
wide bx = wide(b.x) - a.x;
wide by = wide(b.y) - a.y;
wide cx = wide(c.x) - a.x;
wide cy = wide(c.y) - a.y;
wide bc = bx * cy - by * cx;
wide bn = bx * bx + by * by;
wide cn = cx * cx + cy * cy;
for (int i = 0; i < int(in.size()); i++) {
wide px = wide(in[i].x) - a.x;
wide py = wide(in[i].y) - a.y;
wide pn = px * px + py * py;
res[i] =
pn * bc - bn * (px * cy - py * cx) + cn * (px * by - py * bx) == 0;
}
}
return res;
}
} // namespace noya
// NOYA_MINIMUM_ENCLOSING_CIRCLE_HPP
#ifndef NOYA_MINIMUM_ENCLOSING_CIRCLE_HPP
#define NOYA_MINIMUM_ENCLOSING_CIRCLE_HPP 1
/// @complexity Time: Expected O(n).
/// Space: O(n) shuffled points.
#include "noya/geometry/geometry_base.hpp"
#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <random>
#include <type_traits>
#include <vector>
namespace noya {
template <class T> struct minimum_enclosing_circle_result {
circle<long double> val;
std::array<point<T>, 3> sup{};
int cnt = 0;
};
namespace minimum_enclosing_circle_internal {
inline circle<long double> from_diameter(const point<long double> &lhs,
const point<long double> &rhs) {
point<long double> o = (lhs + rhs) / 2.0L;
return {o, distance(o, lhs)};
}
inline circle<long double> from_three(const point<long double> &lhs,
const point<long double> &rhs,
const point<long double> &z) {
point<long double> b = rhs - lhs, c = z - lhs;
long double den = 2.0L * cross(b, c);
if (std::abs(den) <= 1e-24L) {
circle<long double> res = from_diameter(lhs, rhs);
circle<long double> can = from_diameter(lhs, z);
if (can.r1 > res.r1) {
res = can;
}
can = from_diameter(rhs, z);
if (can.r1 > res.r1) {
res = can;
}
return res;
}
// Work relative to lhs to avoid subtracting large absolute squares.
long double bn = norm2(b), cn = norm2(c);
point<long double> v{(bn * c.y - cn * b.y) / den,
(b.x * cn - c.x * bn) / den};
return {lhs + v, std::sqrt(norm2(v))};
}
inline bool contains(const circle<long double> &val,
const point<long double> &can) {
long double r2 = val.r1 * val.r1;
long double tol = 1e-12L * std::max(1.0L, r2);
return norm2(can - val.o) <= r2 + tol;
}
} // namespace minimum_enclosing_circle_internal
/// @brief Return the minimum circle and one boundary support set. Randomly
/// ordering the points makes the expected number of constraint rebuilds
/// linear: an outside point must belong to the new boundary, reducing the
/// remaining problem successively to circles through one, two, then three
/// fixed points. Integral predicate intermediates must fit in signed 128 bits;
/// the circumcircle determinant has degree four in coordinate differences.
template <class T>
minimum_enclosing_circle_result<T>
minimum_enclosing_circle_with_support(const std::vector<point<T>> &in,
std::uint64_t sd = 712367821ULL) {
assert(!in.empty());
struct indexed_point {
point<long double> val;
point<T> org;
};
std::vector<indexed_point> pt;
pt.reserve(in.size());
for (const point<T> &val : in) {
pt.push_back(
{{static_cast<long double>(val.x), static_cast<long double>(val.y)},
val});
}
std::mt19937_64 rng(sd);
std::shuffle(pt.begin(), pt.end(), rng);
using namespace minimum_enclosing_circle_internal;
minimum_enclosing_circle_result<T> ans;
ans.val = {pt[0].val, 0};
ans.sup[0] = pt[0].org;
ans.cnt = 1;
auto has = [&](const indexed_point &can) {
if constexpr (!std::is_integral_v<T>) {
return contains(ans.val, can.val);
} else {
using wide = __int128;
const auto &a = ans.sup[0];
const auto &p = can.org;
if (ans.cnt == 1) {
return p == a;
}
const auto &b = ans.sup[1];
if (ans.cnt == 2) {
wide dx = wide(2) * p.x - a.x - b.x;
wide dy = wide(2) * p.y - a.y - b.y;
wide bx = wide(a.x) - b.x, by = wide(a.y) - b.y;
wide d2 = bx * bx + by * by;
return dx * dx + dy * dy <= d2;
}
const auto &c = ans.sup[2];
wide bx = wide(b.x) - a.x;
wide by = wide(b.y) - a.y;
wide cx = wide(c.x) - a.x;
wide cy = wide(c.y) - a.y;
wide px = wide(p.x) - a.x;
wide py = wide(p.y) - a.y;
wide det = (px * px + py * py) * (bx * cy - by * cx) -
(bx * bx + by * by) * (px * cy - py * cx) +
(cx * cx + cy * cy) * (px * by - py * bx);
// Compare signs without multiplying the fourth-degree determinant.
return (bx * cy - by * cx) > 0 ? det <= 0 : det >= 0;
}
};
for (int i = 1; i < int(pt.size()); i++) {
if (has(pt[i])) {
continue;
}
ans.val = {pt[i].val, 0};
ans.sup[0] = pt[i].org;
ans.cnt = 1;
for (int j = 0; j < i; j++) {
if (has(pt[j])) {
continue;
}
ans.val = from_diameter(pt[i].val, pt[j].val);
ans.sup[0] = pt[i].org;
ans.sup[1] = pt[j].org;
ans.cnt = 2;
for (int k = 0; k < j; k++) {
if (!has(pt[k])) {
long double ta = cross(pt[i].val, pt[j].val, pt[k].val);
if (std::abs(ta) <= 1e-24L) {
int lhs = i;
int rhs = j;
long double bst = norm2(pt[i].val - pt[j].val);
for (auto [l, r] : {std::pair{i, k}, std::pair{j, k}}) {
long double can = norm2(pt[l].val - pt[r].val);
if (can > bst) {
bst = can;
lhs = l;
rhs = r;
}
}
ans.val = from_diameter(pt[lhs].val, pt[rhs].val);
ans.sup[0] = pt[lhs].org;
ans.sup[1] = pt[rhs].org;
ans.cnt = 2;
} else {
ans.val = from_three(pt[i].val, pt[j].val, pt[k].val);
ans.sup[0] = pt[i].org;
ans.sup[1] = pt[j].org;
ans.sup[2] = pt[k].org;
ans.cnt = 3;
}
}
}
}
}
return ans;
}
/// @brief Return the minimum circle containing all points in expected O(n).
template <class T>
circle<long double> minimum_enclosing_circle(const std::vector<point<T>> &in,
std::uint64_t sd = 712367821ULL) {
return minimum_enclosing_circle_with_support(in, sd).val;
}
/// @brief Mark exactly which integral points lie on the minimum circle. Once
/// the randomized construction supplies two or three support points, equality
/// with their circle is tested after clearing all denominators: a doubled
/// midpoint equation for a diameter, or the integer circumcircle det
/// for three non-collinear points.
template <class T>
std::vector<bool>
minimum_enclosing_circle_boundary(const std::vector<point<T>> &in,
std::uint64_t sd = 712367821ULL) {
static_assert(std::is_integral_v<T>);
auto ans = minimum_enclosing_circle_with_support(in, sd);
std::vector<bool> res(in.size());
using wide = __int128;
const auto &a = ans.sup[0];
if (ans.cnt == 1) {
for (int i = 0; i < int(in.size()); i++) {
res[i] = in[i] == a;
}
} else if (ans.cnt == 2) {
const auto &b = ans.sup[1];
wide bx = wide(a.x) - b.x, by = wide(a.y) - b.y;
wide r4 = bx * bx + by * by;
for (int i = 0; i < int(in.size()); i++) {
wide dx = wide(2) * in[i].x - a.x - b.x;
wide dy = wide(2) * in[i].y - a.y - b.y;
res[i] = dx * dx + dy * dy == r4;
}
} else {
const auto &b = ans.sup[1];
const auto &c = ans.sup[2];
wide bx = wide(b.x) - a.x;
wide by = wide(b.y) - a.y;
wide cx = wide(c.x) - a.x;
wide cy = wide(c.y) - a.y;
wide bc = bx * cy - by * cx;
wide bn = bx * bx + by * by;
wide cn = cx * cx + cy * cy;
for (int i = 0; i < int(in.size()); i++) {
wide px = wide(in[i].x) - a.x;
wide py = wide(in[i].y) - a.y;
wide pn = px * px + py * py;
res[i] =
pn * bc - bn * (px * cy - py * cx) + cn * (px * by - py * bx) == 0;
}
}
return res;
}
} // namespace noya
#endif // NOYA_MINIMUM_ENCLOSING_CIRCLE_HPP
/// @complexity Time: Expected O(n).
/// Space: O(n) shuffled points.
/// @complexity Time: O(1) per point/segment primitive, O(n) for polygon area or
/// point location, O(n log n) for convex hull. Space: O(1), except O(n) for hull.
#include <algorithm>
#include <cmath>
#include <optional>
#include <vector>
namespace noya {
/// @brief Two-dimensional point with vector arithmetic and lexicographic order.
template <class T> struct point {
T x{};
T y{};
point() = default;
point(T x_, T y_) : x(x_), y(y_) {}
point &operator+=(const point &rhs) {
x += rhs.x;
y += rhs.y;
return *this;
}
point &operator-=(const point &rhs) {
x -= rhs.x;
y -= rhs.y;
return *this;
}
point &operator*=(const T &scl) {
x *= scl;
y *= scl;
return *this;
}
point &operator/=(const T &scl) {
x /= scl;
y /= scl;
return *this;
}
friend point operator+(point l, const point &r) { return l += r; }
friend point operator-(point l, const point &r) { return l -= r; }
friend point operator*(point val, const T &scl) { return val *= scl; }
friend point operator*(const T &scl, point val) { return val *= scl; }
friend point operator/(point val, const T &scl) { return val /= scl; }
friend bool operator==(const point &, const point &) = default;
friend bool operator<(const point &l, const point &r) {
return l.x < r.x || (l.x == r.x && l.y < r.y);
}
};
/// @brief Circle represented by a center and a nonnegative radius.
template <class Real> struct circle {
point<Real> o;
Real r1{};
};
/// @brief Return the dot product of two vectors.
template <class T> T dot(const point<T> &a, const point<T> &b) {
return a.x * b.x + a.y * b.y;
}
/// @brief Return the signed cross product of two vectors.
template <class T> T cross(const point<T> &a, const point<T> &b) {
return a.x * b.y - a.y * b.x;
}
/// @brief Return cross(a - o1, b - o1).
template <class T>
T cross(const point<T> &o1, const point<T> &a, const point<T> &b) {
return cross(a - o1, b - o1);
}
/// @brief Return the squared Euclidean norm.
template <class T> T norm2(const point<T> &val) { return dot(val, val); }
/// @brief Compare a value with zero using an optional absolute tolerance.
template <class T> int sign(const T &val, const T &eps = T{}) {
return (val > eps) - (val < -eps);
}
/// @brief Return -1, 0, or 1 for a clockwise, collinear, or counter-clockwise
/// turn.
template <class T>
int orientation(const point<T> &a, const point<T> &b, const point<T> &c,
const T &eps = T{}) {
return sign(cross(a, b, c), eps);
}
/// @brief Test whether p lies on the closed segment [a, b].
template <class T>
bool on_segment(const point<T> &p, const point<T> &a, const point<T> &b,
const T &eps = T{}) {
if (orientation(a, b, p, eps) != 0) {
return false;
}
return std::min(a.x, b.x) - eps <= p.x && p.x <= std::max(a.x, b.x) + eps &&
std::min(a.y, b.y) - eps <= p.y && p.y <= std::max(a.y, b.y) + eps;
}
/// @brief Test whether the closed segments [a, b] and [c, d] intersect.
template <class T>
bool segments_intersect(const point<T> &a, const point<T> &b, const point<T> &c,
const point<T> &d, const T &eps = T{}) {
int s1 = orientation(a, b, c, eps);
int s2 = orientation(a, b, d, eps);
int s3 = orientation(c, d, a, eps);
int s4 = orientation(c, d, b, eps);
if (s1 == 0 && on_segment(c, a, b, eps)) {
return true;
}
if (s2 == 0 && on_segment(d, a, b, eps)) {
return true;
}
if (s3 == 0 && on_segment(a, c, d, eps)) {
return true;
}
if (s4 == 0 && on_segment(b, c, d, eps)) {
return true;
}
return s1 * s2 < 0 && s3 * s4 < 0;
}
/// @brief Intersect the infinite lines through (a, b) and (c, d), returning
/// nullopt when they are parallel or coincident.
template <class T>
std::optional<point<long double>>
line_intersection(const point<T> &a, const point<T> &b, const point<T> &c,
const point<T> &d, long double eps = 0) {
point<long double> lhs{static_cast<long double>(a.x),
static_cast<long double>(a.y)};
point<long double> b1{static_cast<long double>(b.x),
static_cast<long double>(b.y)};
point<long double> z{static_cast<long double>(c.x),
static_cast<long double>(c.y)};
point<long double> d1{static_cast<long double>(d.x),
static_cast<long double>(d.y)};
point<long double> da = b1 - lhs;
point<long double> db = d1 - z;
long double den = cross(da, db);
if (std::abs(den) <= eps) {
return std::nullopt;
}
long double rat = cross(z - lhs, db) / den;
return lhs + da * rat;
}
/// @brief Return the convex hull in counter-clockwise order without repetition.
/// Starts at the lexicographically smallest point; kc retains boundary points.
/// Coordinate differences and determinant arithmetic must fit in T.
template <class T>
std::vector<point<T>> convex_hull(std::vector<point<T>> pt, bool kc = false) {
std::sort(pt.begin(), pt.end());
pt.erase(std::unique(pt.begin(), pt.end()), pt.end());
if (pt.size() <= 1) {
return pt;
}
if (kc) {
bool col = true;
for (int i = 2; i < int(pt.size()) && col; i++) {
col = orientation(pt[0], pt[1], pt[i]) == 0;
}
if (col) return pt;
}
std::vector<point<T>> res;
res.reserve(pt.size() + 1);
for (const point<T> &p : pt) {
while (res.size() >= 2) {
int rot = orientation(res[res.size() - 2], res.back(), p);
if (rot > 0 || (kc && rot == 0)) {
break;
}
res.pop_back();
}
res.push_back(p);
}
auto k = res.size(); // Keep the lower chain while building the upper chain.
for (int i = int(pt.size()) - 2; i >= 0; i--) {
while (res.size() > k) {
int rot = orientation(res[res.size() - 2], res.back(), pt[i]);
if (rot > 0 || (kc && rot == 0)) {
break;
}
res.pop_back();
}
res.push_back(pt[i]);
}
res.pop_back();
return res;
}
/// @brief Return twice the signed area of a polygon.
/// Coordinate differences, products and the accumulated area must fit in T.
template <class T> T polygon_area2(const std::vector<point<T>> &pg) {
T res{};
if (pg.size() < 3) return res;
auto prev = pg[1] - pg[0];
for (int i = 2; i < int(pg.size()); i++) {
auto cur = pg[i] - pg[0];
res += cross(prev, cur);
prev = cur;
}
return res;
}
/// @brief Classify a point relative to a polygon: -1 outside, 0 boundary, 1
/// inside.
template <class T>
int point_in_polygon(const point<T> &p, const std::vector<point<T>> &pg) {
bool in = false;
for (int i = 0; i < int(pg.size()); i++) {
point<T> a = pg[i];
point<T> b = pg[(i + 1) % pg.size()];
if (on_segment(p, a, b)) {
return 0;
}
if (a.y <= p.y && p.y < b.y && orientation(a, b, p) > 0) {
in = !in;
}
if (b.y <= p.y && p.y < a.y && orientation(a, b, p) < 0) {
in = !in;
}
}
return in ? 1 : -1;
}
/// @brief Return the Euclidean distance between two points.
template <class T> long double distance(const point<T> &a, const point<T> &b) {
return std::hypot(
static_cast<long double>(a.x) - static_cast<long double>(b.x),
static_cast<long double>(a.y) - static_cast<long double>(b.y));
}
} // namespace noya
#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <random>
#include <type_traits>
#include <vector>
namespace noya {
template <class T> struct minimum_enclosing_circle_result {
circle<long double> val;
std::array<point<T>, 3> sup{};
int cnt = 0;
};
namespace minimum_enclosing_circle_internal {
inline circle<long double> from_diameter(const point<long double> &lhs,
const point<long double> &rhs) {
point<long double> o = (lhs + rhs) / 2.0L;
return {o, distance(o, lhs)};
}
inline circle<long double> from_three(const point<long double> &lhs,
const point<long double> &rhs,
const point<long double> &z) {
point<long double> b = rhs - lhs, c = z - lhs;
long double den = 2.0L * cross(b, c);
if (std::abs(den) <= 1e-24L) {
circle<long double> res = from_diameter(lhs, rhs);
circle<long double> can = from_diameter(lhs, z);
if (can.r1 > res.r1) {
res = can;
}
can = from_diameter(rhs, z);
if (can.r1 > res.r1) {
res = can;
}
return res;
}
// Work relative to lhs to avoid subtracting large absolute squares.
long double bn = norm2(b), cn = norm2(c);
point<long double> v{(bn * c.y - cn * b.y) / den,
(b.x * cn - c.x * bn) / den};
return {lhs + v, std::sqrt(norm2(v))};
}
inline bool contains(const circle<long double> &val,
const point<long double> &can) {
long double r2 = val.r1 * val.r1;
long double tol = 1e-12L * std::max(1.0L, r2);
return norm2(can - val.o) <= r2 + tol;
}
} // namespace minimum_enclosing_circle_internal
/// @brief Return the minimum circle and one boundary support set. Randomly
/// ordering the points makes the expected number of constraint rebuilds
/// linear: an outside point must belong to the new boundary, reducing the
/// remaining problem successively to circles through one, two, then three
/// fixed points. Integral predicate intermediates must fit in signed 128 bits;
/// the circumcircle determinant has degree four in coordinate differences.
template <class T>
minimum_enclosing_circle_result<T>
minimum_enclosing_circle_with_support(const std::vector<point<T>> &in,
std::uint64_t sd = 712367821ULL) {
assert(!in.empty());
struct indexed_point {
point<long double> val;
point<T> org;
};
std::vector<indexed_point> pt;
pt.reserve(in.size());
for (const point<T> &val : in) {
pt.push_back(
{{static_cast<long double>(val.x), static_cast<long double>(val.y)},
val});
}
std::mt19937_64 rng(sd);
std::shuffle(pt.begin(), pt.end(), rng);
using namespace minimum_enclosing_circle_internal;
minimum_enclosing_circle_result<T> ans;
ans.val = {pt[0].val, 0};
ans.sup[0] = pt[0].org;
ans.cnt = 1;
auto has = [&](const indexed_point &can) {
if constexpr (!std::is_integral_v<T>) {
return contains(ans.val, can.val);
} else {
using wide = __int128;
const auto &a = ans.sup[0];
const auto &p = can.org;
if (ans.cnt == 1) {
return p == a;
}
const auto &b = ans.sup[1];
if (ans.cnt == 2) {
wide dx = wide(2) * p.x - a.x - b.x;
wide dy = wide(2) * p.y - a.y - b.y;
wide bx = wide(a.x) - b.x, by = wide(a.y) - b.y;
wide d2 = bx * bx + by * by;
return dx * dx + dy * dy <= d2;
}
const auto &c = ans.sup[2];
wide bx = wide(b.x) - a.x;
wide by = wide(b.y) - a.y;
wide cx = wide(c.x) - a.x;
wide cy = wide(c.y) - a.y;
wide px = wide(p.x) - a.x;
wide py = wide(p.y) - a.y;
wide det = (px * px + py * py) * (bx * cy - by * cx) -
(bx * bx + by * by) * (px * cy - py * cx) +
(cx * cx + cy * cy) * (px * by - py * bx);
// Compare signs without multiplying the fourth-degree determinant.
return (bx * cy - by * cx) > 0 ? det <= 0 : det >= 0;
}
};
for (int i = 1; i < int(pt.size()); i++) {
if (has(pt[i])) {
continue;
}
ans.val = {pt[i].val, 0};
ans.sup[0] = pt[i].org;
ans.cnt = 1;
for (int j = 0; j < i; j++) {
if (has(pt[j])) {
continue;
}
ans.val = from_diameter(pt[i].val, pt[j].val);
ans.sup[0] = pt[i].org;
ans.sup[1] = pt[j].org;
ans.cnt = 2;
for (int k = 0; k < j; k++) {
if (!has(pt[k])) {
long double ta = cross(pt[i].val, pt[j].val, pt[k].val);
if (std::abs(ta) <= 1e-24L) {
int lhs = i;
int rhs = j;
long double bst = norm2(pt[i].val - pt[j].val);
for (auto [l, r] : {std::pair{i, k}, std::pair{j, k}}) {
long double can = norm2(pt[l].val - pt[r].val);
if (can > bst) {
bst = can;
lhs = l;
rhs = r;
}
}
ans.val = from_diameter(pt[lhs].val, pt[rhs].val);
ans.sup[0] = pt[lhs].org;
ans.sup[1] = pt[rhs].org;
ans.cnt = 2;
} else {
ans.val = from_three(pt[i].val, pt[j].val, pt[k].val);
ans.sup[0] = pt[i].org;
ans.sup[1] = pt[j].org;
ans.sup[2] = pt[k].org;
ans.cnt = 3;
}
}
}
}
}
return ans;
}
/// @brief Return the minimum circle containing all points in expected O(n).
template <class T>
circle<long double> minimum_enclosing_circle(const std::vector<point<T>> &in,
std::uint64_t sd = 712367821ULL) {
return minimum_enclosing_circle_with_support(in, sd).val;
}
/// @brief Mark exactly which integral points lie on the minimum circle. Once
/// the randomized construction supplies two or three support points, equality
/// with their circle is tested after clearing all denominators: a doubled
/// midpoint equation for a diameter, or the integer circumcircle det
/// for three non-collinear points.
template <class T>
std::vector<bool>
minimum_enclosing_circle_boundary(const std::vector<point<T>> &in,
std::uint64_t sd = 712367821ULL) {
static_assert(std::is_integral_v<T>);
auto ans = minimum_enclosing_circle_with_support(in, sd);
std::vector<bool> res(in.size());
using wide = __int128;
const auto &a = ans.sup[0];
if (ans.cnt == 1) {
for (int i = 0; i < int(in.size()); i++) {
res[i] = in[i] == a;
}
} else if (ans.cnt == 2) {
const auto &b = ans.sup[1];
wide bx = wide(a.x) - b.x, by = wide(a.y) - b.y;
wide r4 = bx * bx + by * by;
for (int i = 0; i < int(in.size()); i++) {
wide dx = wide(2) * in[i].x - a.x - b.x;
wide dy = wide(2) * in[i].y - a.y - b.y;
res[i] = dx * dx + dy * dy == r4;
}
} else {
const auto &b = ans.sup[1];
const auto &c = ans.sup[2];
wide bx = wide(b.x) - a.x;
wide by = wide(b.y) - a.y;
wide cx = wide(c.x) - a.x;
wide cy = wide(c.y) - a.y;
wide bc = bx * cy - by * cx;
wide bn = bx * bx + by * by;
wide cn = cx * cx + cy * cy;
for (int i = 0; i < int(in.size()); i++) {
wide px = wide(in[i].x) - a.x;
wide py = wide(in[i].y) - a.y;
wide pn = px * px + py * py;
res[i] =
pn * bc - bn * (px * cy - py * cx) + cn * (px * by - py * bx) == 0;
}
}
return res;
}
} // namespace noya
// NOYA_MINIMUM_ENCLOSING_CIRCLE_HPP