sotanishy's code snippets for competitive programming
View the Project on GitHub sotanishy/cp-library-cpp
#include "geometry/triangle.hpp"
#pragma once #include "bisector.hpp" #include "geometry.hpp" #include "intersection.hpp" Vec centroid(const Vec& A, const Vec& B, const Vec& C) { assert(ccw(A, B, C) != 0); return (A + B + C) / T(3); } Vec incenter(const Vec& A, const Vec& B, const Vec& C) { assert(ccw(A, B, C) != 0); T a = std::abs(B - C); T b = std::abs(C - A); T c = std::abs(A - B); return (a * A + b * B + c * C) / (a + b + c); } Vec circumcenter(const Vec& A, const Vec& B, const Vec& C) { assert(ccw(A, B, C) != 0); return intersection(bisector(Segment(A, B)), bisector(Segment(A, C))); } // large error but beautiful // Vec circumcenter(const Vec& A, const Vec& B, const Vec& C) { // assert(ccw(A, B, C) != 0); // Vec p = C - B, q = A - C, r = B - A; // T a = std::norm(p) * dot(q, r); // T b = std::norm(q) * dot(r, p); // T c = std::norm(r) * dot(p, q); // return (a*A + b*B + c*C) / (a + b + c); // }
#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 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 4 "geometry/bisector.hpp" Line bisector(const Segment& s) { auto m = (s.p1 + s.p2) / T(2); return Line(m, m + Vec(-s.dir().imag(), s.dir().real())); } std::pair<Line, Line> bisector(const Line& l, const Line& m) { // parallel if (eq(cross(l.dir(), m.dir()), 0)) { auto n = Line(l.p1, l.p1 + perp(l.dir())); auto p = intersection(n, m); auto m = (l.p1 + p) / T(2); return {Line(m, m + l.dir()), Line()}; } auto p = intersection(l, m); T ang = (std::arg(l.dir()) + std::arg(m.dir())) / T(2); auto b1 = Line(p, p + std::polar(T(1), ang)); auto b2 = Line(p, p + std::polar(T(1), ang + PI / T(2))); return {b1, b2}; } #line 5 "geometry/triangle.hpp" Vec centroid(const Vec& A, const Vec& B, const Vec& C) { assert(ccw(A, B, C) != 0); return (A + B + C) / T(3); } Vec incenter(const Vec& A, const Vec& B, const Vec& C) { assert(ccw(A, B, C) != 0); T a = std::abs(B - C); T b = std::abs(C - A); T c = std::abs(A - B); return (a * A + b * B + c * C) / (a + b + c); } Vec circumcenter(const Vec& A, const Vec& B, const Vec& C) { assert(ccw(A, B, C) != 0); return intersection(bisector(Segment(A, B)), bisector(Segment(A, C))); } // large error but beautiful // Vec circumcenter(const Vec& A, const Vec& B, const Vec& C) { // assert(ccw(A, B, C) != 0); // Vec p = C - B, q = A - C, r = B - A; // T a = std::norm(p) * dot(q, r); // T b = std::norm(q) * dot(r, p); // T c = std::norm(r) * dot(p, q); // return (a*A + b*B + c*C) / (a + b + c); // }