sotanishy's code snippets for competitive programming
View the Project on GitHub sotanishy/cp-library-cpp
#define PROBLEM "http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=1283" #define ERROR 0.00001 #include "../../geometry/geometry.hpp" #include "../../geometry/polygon.hpp" #include <bits/stdc++.h> using namespace std; using ll = long long; int main() { ios_base::sync_with_stdio(false); cin.tie(nullptr); cout << fixed << setprecision(15); while (true) { int n; cin >> n; if (n == 0) break; vector<Vec> pts(n); for (auto& p : pts) cin >> p; double lb = 0, ub = 1e9; for (int k = 0; k < 100; ++k) { double m = (lb + ub) / 2; vector<pair<Vec, Vec>> hps; for (int i = 0; i < n; ++i) { auto d = pts[(i+1)%n] - pts[i]; auto e = Vec(-d.imag(), d.real()); auto q = pts[i] + e/abs(e)*m; hps.push_back({e, q}); } if (halfplane_intersection(hps).empty()) { ub = m; } else { lb = m; } } cout << lb << "\n"; } }
#line 1 "test/aoj/1283.test.cpp" #define PROBLEM "http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=1283" #define ERROR 0.00001 #line 2 "geometry/geometry.hpp" #include <algorithm> #include <cassert> #include <cmath> #include <complex> #include <iostream> #include <numbers> #include <numeric> #include <vector> // note that if T is of an integer type, std::abs does not work using T = double; using Vec = std::complex<T>; std::istream& operator>>(std::istream& is, Vec& p) { T x, y; is >> x >> y; p = {x, y}; return is; } T dot(const Vec& a, const Vec& b) { return (std::conj(a) * b).real(); } T cross(const Vec& a, const Vec& b) { return (std::conj(a) * b).imag(); } constexpr T PI = std::numbers::pi_v<T>; constexpr T eps = 1e-10; inline bool eq(T a, T b) { return std::abs(a - b) <= eps; } inline bool eq(Vec a, Vec b) { return std::abs(a - b) <= eps; } inline bool lt(T a, T b) { return a < b - eps; } inline bool leq(T a, T b) { return a <= b + eps; } struct Line { Vec p1, p2; Line() = default; Line(const Vec& p1, const Vec& p2) : p1(p1), p2(p2) {} Vec dir() const { return p2 - p1; } }; struct Segment : Line { using Line::Line; }; struct Circle { Vec c; T r; Circle() = default; Circle(const Vec& c, T r) : c(c), r(r) {} }; using Polygon = std::vector<Vec>; Vec rot(const Vec& a, T ang) { return a * Vec(std::cos(ang), std::sin(ang)); } Vec perp(const Vec& a) { return Vec(-a.imag(), a.real()); } Vec projection(const Line& l, const Vec& p) { return l.p1 + dot(p - l.p1, l.dir()) * l.dir() / std::norm(l.dir()); } Vec reflection(const Line& l, const Vec& p) { return T(2) * projection(l, p) - p; } // 0: collinear // 1: counter-clockwise // -1: clockwise int ccw(const Vec& a, const Vec& b, const Vec& c) { if (eq(cross(b - a, c - a), 0)) return 0; if (lt(cross(b - a, c - a), 0)) return -1; return 1; } void sort_by_arg(std::vector<Vec>& pts) { std::ranges::sort(pts, [&](auto& p, auto& q) { if ((p.imag() < 0) != (q.imag() < 0)) return (p.imag() < 0); if (cross(p, q) == 0) { if (p == Vec(0, 0)) return !(q.imag() < 0 || (q.imag() == 0 && q.real() > 0)); if (q == Vec(0, 0)) return (p.imag() < 0 || (p.imag() == 0 && p.real() > 0)); return (p.real() > q.real()); } return (cross(p, q) > 0); }); } #line 2 "geometry/polygon.hpp" #include <deque> #line 3 "geometry/intersect.hpp" bool intersect(const Segment& s, const Vec& p) { Vec u = s.p1 - p, v = s.p2 - p; return eq(cross(u, v), 0) && leq(dot(u, v), 0); } // 0: outside // 1: on the border // 2: inside int intersect(const Polygon& poly, const Vec& p) { const int n = poly.size(); bool in = 0; for (int i = 0; i < n; ++i) { auto a = poly[i] - p, b = poly[(i + 1) % n] - p; if (eq(cross(a, b), 0) && (lt(dot(a, b), 0) || eq(dot(a, b), 0))) return 1; if (a.imag() > b.imag()) std::swap(a, b); if (leq(a.imag(), 0) && lt(0, b.imag()) && lt(cross(a, b), 0)) in ^= 1; } return in ? 2 : 0; } int intersect(const Segment& s, const Segment& t) { auto a = s.p1, b = s.p2; auto c = t.p1, d = t.p2; if (ccw(a, b, c) != ccw(a, b, d) && ccw(c, d, a) != ccw(c, d, b)) return 2; if (intersect(s, c) || intersect(s, d) || intersect(t, a) || intersect(t, b)) return 1; return 0; } // true if they have positive area in common or touch on the border bool intersect(const Polygon& poly1, const Polygon& poly2) { const int n = poly1.size(); const int m = poly2.size(); for (int i = 0; i < n; ++i) { for (int j = 0; j < m; ++j) { if (intersect(Segment(poly1[i], poly1[(i + 1) % n]), Segment(poly2[j], poly2[(j + 1) % m]))) { return true; } } } return intersect(poly1, poly2[0]) || intersect(poly2, poly1[0]); } // 0: inside // 1: inscribe // 2: intersect // 3: circumscribe // 4: outside int intersect(const Circle& c1, const Circle& c2) { T d = std::abs(c1.c - c2.c); if (lt(d, std::abs(c2.r - c1.r))) return 0; if (eq(d, std::abs(c2.r - c1.r))) return 1; if (eq(c1.r + c2.r, d)) return 3; if (lt(c1.r + c2.r, d)) return 4; return 2; } #line 4 "geometry/dist.hpp" T dist(const Line& l, const Vec& p) { return std::abs(cross(p - l.p1, l.dir())) / std::abs(l.dir()); } T dist(const Segment& s, const Vec& p) { if (lt(dot(p - s.p1, s.dir()), 0)) return std::abs(p - s.p1); if (lt(dot(p - s.p2, -s.dir()), 0)) return std::abs(p - s.p2); return std::abs(cross(p - s.p1, s.dir())) / std::abs(s.dir()); } T dist(const Segment& s, const Segment& t) { if (intersect(s, t)) return T(0); return std::min( {dist(s, t.p1), dist(s, t.p2), dist(t, s.p1), dist(t, s.p2)}); } #line 4 "geometry/intersection.hpp" Vec intersection(const Line& l, const Line& m) { assert(!eq(cross(l.dir(), m.dir()), 0)); // not parallel Vec r = m.p1 - l.p1; return l.p1 + cross(m.dir(), r) / cross(m.dir(), l.dir()) * l.dir(); } std::vector<Vec> intersection(const Circle& c, const Line& l) { T d = dist(l, c.c); if (lt(c.r, d)) return {}; // no intersection Vec e1 = l.dir() / std::abs(l.dir()); Vec e2 = perp(e1); if (ccw(c.c, l.p1, l.p2) == 1) e2 *= -1; if (eq(c.r, d)) return {c.c + d * e2}; // tangent T t = std::sqrt(c.r * c.r - d * d); return {c.c + d * e2 + t * e1, c.c + d * e2 - t * e1}; } std::vector<Vec> intersection(const Circle& c1, const Circle& c2) { T d = std::abs(c1.c - c2.c); if (lt(c1.r + c2.r, d)) return {}; // outside Vec e1 = (c2.c - c1.c) / std::abs(c2.c - c1.c); Vec e2 = perp(e1); if (lt(d, std::abs(c2.r - c1.r))) return {}; // contain if (eq(d, std::abs(c2.r - c1.r))) return {c1.c + c1.r * e1}; // tangent T x = (c1.r * c1.r - c2.r * c2.r + d * d) / (2 * d); T y = std::sqrt(c1.r * c1.r - x * x); return {c1.c + x * e1 + y * e2, c1.c + x * e1 - y * e2}; } T area_intersection(const Circle& c1, const Circle& c2) { T d = std::abs(c2.c - c1.c); if (leq(c1.r + c2.r, d)) return 0; // outside if (leq(d, std::abs(c2.r - c1.r))) { // inside T r = std::min(c1.r, c2.r); return PI * r * r; } T ans = 0; T a; a = std::acos((c1.r * c1.r + d * d - c2.r * c2.r) / (2 * c1.r * d)); ans += c1.r * c1.r * (a - std::sin(a) * std::cos(a)); a = std::acos((c2.r * c2.r + d * d - c1.r * c1.r) / (2 * c2.r * d)); ans += c2.r * c2.r * (a - std::sin(a) * std::cos(a)); return ans; } #line 6 "geometry/polygon.hpp" T area(const Polygon& poly) { const int n = poly.size(); T res = 0; for (int i = 0; i < n; ++i) { res += cross(poly[i], poly[(i + 1) % n]); } return std::abs(res) / T(2); } bool is_convex(const Polygon& poly) { int n = poly.size(); for (int i = 0; i < n; ++i) { if (lt(cross(poly[(i + 1) % n] - poly[i], poly[(i + 2) % n] - poly[(i + 1) % n]), 0)) { return false; } } return true; } Polygon convex_cut(const Polygon& poly, const Line& l) { const int n = poly.size(); Polygon res; for (int i = 0; i < n; ++i) { auto p = poly[i], q = poly[(i + 1) % n]; if (ccw(l.p1, l.p2, p) != -1) { if (res.empty() || !eq(res.back(), p)) { res.push_back(p); } } if (ccw(l.p1, l.p2, p) * ccw(l.p1, l.p2, q) < 0) { auto c = intersection(Line(p, q), l); if (res.empty() || !eq(res.back(), c)) { res.push_back(c); } } } return res; } Polygon halfplane_intersection(std::vector<std::pair<Vec, Vec>> hps) { using Hp = std::pair<Vec, Vec>; // (normal vector, a point on the border) auto intersection = [&](const Hp& l1, const Hp& l2) -> Vec { auto d = l2.second - l1.second; return l1.second + (dot(d, l2.first) / cross(l1.first, l2.first)) * perp(l1.first); }; // check if the halfplane h contains the point p auto contains = [&](const Hp& h, const Vec& p) -> bool { return dot(p - h.second, h.first) > 0; }; constexpr T INF = 1e15; hps.emplace_back(Vec(1, 0), Vec(-INF, 0)); // -INF <= x hps.emplace_back(Vec(-1, 0), Vec(INF, 0)); // x <= INF hps.emplace_back(Vec(0, 1), Vec(0, -INF)); // -INF <= y hps.emplace_back(Vec(0, -1), Vec(0, INF)); // y <= INF std::ranges::sort(hps, {}, [&](auto& h) { return std::arg(h.first); }); std::deque<Hp> dq; int len = 0; for (auto& hp : hps) { while (len > 1 && !contains(hp, intersection(dq[len - 1], dq[len - 2]))) { dq.pop_back(); --len; } while (len > 1 && !contains(hp, intersection(dq[0], dq[1]))) { dq.pop_front(); --len; } // parallel if (len > 0 && eq(cross(dq[len - 1].first, hp.first), 0)) { // opposite if (lt(dot(dq[len - 1].first, hp.first), 0)) { return {}; } // same if (!contains(hp, dq[len - 1].second)) { dq.pop_back(); --len; } else continue; } dq.push_back(hp); ++len; } while (len > 2 && !contains(dq[0], intersection(dq[len - 1], dq[len - 2]))) { dq.pop_back(); --len; } while (len > 2 && !contains(dq[len - 1], intersection(dq[0], dq[1]))) { dq.pop_front(); --len; } if (len < 3) return {}; std::vector<Vec> poly(len); for (int i = 0; i < len - 1; ++i) { poly[i] = intersection(dq[i], dq[i + 1]); } poly[len - 1] = intersection(dq[len - 1], dq[0]); return poly; } class PolygonContainment { public: explicit PolygonContainment(Polygon poly) : poly(poly) {} // 0: outside // 1: on the border // 2: inside int query(Vec pt) { auto c1 = cross(poly[1] - poly[0], pt - poly[0]); auto c2 = cross(poly.back() - poly[0], pt - poly[0]); if (lt(c1, 0) || lt(0, c2)) return 0; int lb = 1, ub = (int)poly.size() - 1; while (ub - lb > 1) { int m = std::midpoint(lb, ub); if (leq(0, cross(poly[m] - poly[0], pt - poly[0]))) lb = m; else ub = m; } auto c = cross(poly[lb] - pt, poly[ub] - pt); if (lt(c, 0)) return 0; if (eq(c, 0) || eq(c1, 0) || eq(c2, 0)) return 1; return 2; } private: Polygon poly; }; #line 6 "test/aoj/1283.test.cpp" #include <bits/stdc++.h> using namespace std; using ll = long long; int main() { ios_base::sync_with_stdio(false); cin.tie(nullptr); cout << fixed << setprecision(15); while (true) { int n; cin >> n; if (n == 0) break; vector<Vec> pts(n); for (auto& p : pts) cin >> p; double lb = 0, ub = 1e9; for (int k = 0; k < 100; ++k) { double m = (lb + ub) / 2; vector<pair<Vec, Vec>> hps; for (int i = 0; i < n; ++i) { auto d = pts[(i+1)%n] - pts[i]; auto e = Vec(-d.imag(), d.real()); auto q = pts[i] + e/abs(e)*m; hps.push_back({e, q}); } if (halfplane_intersection(hps).empty()) { ub = m; } else { lb = m; } } cout << lb << "\n"; } }