splitting geometry content

code/geometry/2d.cpp

@@ -87,31 +87,6 @@ T rect_overlapping_area(pair<T, T> l1, pair<T, T> r1, pair<T, T> l2, pair<T, T>
   return a1 + a2 - ai;
 }
 
-template <typename T = double> struct Line {
-  T a, b, c; // ax + by = c
-
-  Line() : a(0), b(0), c(0) {}
-  Line(T a_, T b_, T c_) : a(a_), b(b_), c(c_) {}
-  Line(const Point<T> &p, const Point<T> &q) {
-    a = p.y - q.y;
-    b = q.x - p.x;
-    c = p.cross(q);
-  }
-
-  bool operator==(const Line &l) const {
-    return eq(a, l.a) and eq(b, l.b) and eq(c, l.c);
-  }
-
-  bool parallel(const Line &l) const { return eq(a * l.b, b * l.a); }
-
-  optional<Point<T>> intersection(const Line &l) const {
-    if (parallel(l))
-      return {};
-    auto det = a * l.b - b * l.a;
-    return Point<T>((b * l.c - c * l.b) / det, (c * l.a - a * l.c) / det);
-  }
-};
-
 template <typename T = double> struct Segment {
   Point<T> p, q;
 

code/geometry/circle.cpp

@@ -0,0 +1,42 @@
+template <typename T> struct Circle {
+  Point<T> c; T r;
+  Circle() : c(), r(0) {}
+  Circle(Point<T> c_, T r_) : c(c_), r(r_) {}
+  template <typename U> bool contains(const Point<U> &p) { return c.dist(p) <= r; }
+  pair<Point<T>, Point<T>> tangent_points(const Point<T> &p) {
+    double d1 = (p - c).norm(), theta = asin(r / d1);
+    assert(d1 >= r);
+    auto p1 = (c-p).rotated(theta) * sqrt(d1*d1 - r*r) / d1 + p;
+    auto p2 = (c-p).rotated(-theta) * sqrt(d1*d1 - r*r) / d1 + p;
+    return {p1, p2};
+  }
+  friend ostream& operator<<(ostream &os, const Circle &cr) { return os << cr.c << ' ' << cr.r; }
+  friend istream& operator>>(istream &is, Circle &cr) { return is >> cr.c >> cr.r; }
+};
+
+template <typename T> Circle<T> circumcircle(const Point<T> &a, const Point<T> &b, const Point<T> &c) {
+  Point<T> u = {(b-a).y, (a-b).x};
+  Point<T> v = {(c-a).y, (a-c).x};
+  Point<T> n = (c-b)/2.;
+  auto t = u.cross(n) / v.cross(u);
+  Point<T> center = ((a+c)/2.)+v*t;
+  return {center, center.dist(a)};
+}
+
+template <typename T> Circle<double> minimum_enclosing_circle(vector<Point<T>> &p) {
+  random_shuffle(begin(p), end(p));
+  Circle<double> res(p[0], 0);
+  for (int i = 0; i < (int)size(p); i++) {
+    if (res.contains(p[i])) continue;
+    res = {p[i], 0};
+    for (int j = 0; j < i; j++) {
+      if (res.contains(p[j])) continue;
+      res = {(p[j]+p[i])/2., p[j].dist(p[i])/2.};
+      for (int k = 0; k < j; k++) {
+        if (res.contains(p[k])) continue;
+        res = circumcircle(p[j], p[i], p[k]);
+      }
+    }
+  }
+  return res;
+}

code/geometry/line.cpp

@@ -0,0 +1,21 @@
+template <typename T = double> struct Line {
+  T a, b, c; // ax + by = c
+
+  Line() : a(0), b(0), c(0) {}
+  Line(T a_, T b_, T c_) : a(a_), b(b_), c(c_) {}
+  Line(const Point<T> &p, const Point<T> &q) {
+    a = p.y - q.y;
+    b = q.x - p.x;
+    c = p.cross(q);
+  }
+  double eval(T x) const { return (c-a*x)/b; }
+  double dist(const Point<T> &p) const { return abs(a*p.x+b*p.y+c)/(sqrt(a*a+b*b)); }
+  bool parallel(const Line &l) const { return eq(a * l.b, b * l.a); }
+  optional<Point<T>> intersection(const Line &l) const {
+    if (parallel(l))
+      return {};
+    auto det = a * l.b - b * l.a;
+    return Point<T>((b * l.c - c * l.b) / det, (c * l.a - a * l.c) / det);
+  }
+  bool operator==(const Line &l) const { return eq(a, l.a) and eq(b, l.b) and eq(c, l.c); }
+};

code/geometry/point.cpp

@@ -1,46 +1,40 @@
 template<typename T=double> struct Point {
   T x, y;
-
   Point() : x(0), y(0) {}
   Point(T x_, T y_) : x(x_), y(y_) {}
-
-  inline Point operator+(const Point &p) const { return Point(x + p.x, y + p.y); }
-  inline Point operator-(const Point &p) const { return Point(x - p.x, y - p.y); }
-  inline Point operator+(const T &k) const { return Point(x + k, y + k); }
-  inline Point operator-(const T &k) const { return Point(x - k, y - k); }
-  inline Point operator*(const T &k) const { return Point(x * k, y * k); }
-  inline Point operator/(const T &k) const { return Point(x / k, y / k); }
-
-  inline Point& operator+=(const Point &p) { x += p.x, y += p.y; return *this; }
-  inline Point& operator-=(const Point &p) { x -= p.x, y -= p.y; return *this; }
-  inline Point& operator+=(const T &k) { x += k, y += k; return *this; }
-  inline Point& operator-=(const T &k) { x -= k, y -= k; return *this; }
-  inline Point& operator*=(const T &k) { x *= k, y *= k; return *this; }
-  inline Point& operator/=(const T &k) { x /= k, y /= k; return *this; }
-
-  inline bool operator==(const Point &p) const { return eq(x, p.x) and eq(y, p.y); }
-  inline bool operator<(const Point &p) const { return eq(x, p.x) ? y < p.y : x < p.x; }
-  inline bool operator>(const Point &p) const { return eq(x, p.x) ? y > p.y : x > p.x; }
-  inline bool operator<=(const Point &p) const { return *this == p or *this < p; }
-  inline bool operator>=(const Point &p) const { return *this == p or *this > p; }
-
+  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+(const T &k) const { return Point(x + k, y + k); }
+  Point operator-(const T &k) const { return Point(x - k, y - k); }
+  Point operator*(const T &k) const { return Point(x * k, y * k); }
+  Point operator/(const T &k) const { return Point(x / k, y / k); }
+  Point& operator+=(const Point &p) { x += p.x, y += p.y; return *this; }
+  Point& operator-=(const Point &p) { x -= p.x, y -= p.y; return *this; }
+  Point& operator+=(const T &k) { x += k, y += k; return *this; }
+  Point& operator-=(const T &k) { x -= k, y -= k; return *this; }
+  Point& operator*=(const T &k) { x *= k, y *= k; return *this; }
+  Point& operator/=(const T &k) { x /= k, y /= k; return *this; }
+  bool operator==(const Point &p) const { return eq(x, p.x) and eq(y, p.y); }
+  bool operator<(const Point &p) const { return eq(x, p.x) ? y < p.y : x < p.x; }
+  bool operator>(const Point &p) const { return eq(x, p.x) ? y > p.y : x > p.x; }
+  bool operator<=(const Point &p) const { return *this == p or *this < p; }
+  bool operator>=(const Point &p) const { return *this == p or *this > p; }
   friend ostream& operator<<(ostream &os, const Point &p) { return os << p.x << ' ' << p.y; }
   friend istream& operator>>(istream &is, Point &p) { return is >> p.x >> p.y; }
-
   template<typename U> void rotate(U rad) {
     tie(x, y) = make_pair(x * cos(rad) - y * sin(rad), x * sin(rad) + y * cos(rad));
   }
   template<typename U> Point<U> rotated(U rad) const {
     return Point<U>(x * cos(rad) - y * sin(rad), x * sin(rad) + y * cos(rad));
   }
-  inline T dot(const Point &p) const { return x * p.x + y * p.y; }
-  inline T cross(const Point &p) const { return x * p.y - y * p.x; }
-  inline T cross(const Point &a, const Point &b) const { return (a - *this).cross(b - *this); }
-  inline T dist2() const { return x * x + y * y; }
-  inline double dist() const { return hypot(x, y); }
-  inline double angle() const { return atan2(y, x); }
-  inline double norm() const { return sqrt(dot(*this)); }
-  inline Point rot90() const { return Point(-y, x); }
-  inline Point to(const Point &p) const { return p - *this; }
+  T dot(const Point &p) const { return x * p.x + y * p.y; }
+  T cross(const Point &p) const { return x * p.y - y * p.x; }
+  T cross(const Point &a, const Point &b) const { return (a - *this).cross(b - *this); }
+  T dist2() const { return x * x + y * y; }
+  double dist() const { return hypot(x, y); }
+  double dist(const Point &p) const { return to(p).dist(); }
+  double angle() const { return atan2(y, x); }
+  double norm() const { return sqrt(dot(*this)); }
+  Point rot90() const { return Point(-y, x); }
+  Point to(const Point &p) const { return p - *this; }
 };
-