円
(geometry/circle.hpp)

• View this file on GitHub
• Last update: 2026-08-11 23:01:09+09:00
• Include: #include "geometry/circle.hpp"

二次元平面上の円に関する交点，接線，面積計算を扱う．

二次元平面上の円を中心 p と半径 r で表す．

• intersect(c1, c2)：二円の位置関係を返す．内包，内接，二点で交差，外接，離れている場合にそれぞれ $0,1,2,3,4$．
• cross_point_cl(c, l)：円 c と直線 l の交点を返す．接する場合は同じ点を二つ返す．
• cross_point_cc(c1, c2)：二円の交点を返す．接する場合は同じ点を二つ返す．
• tangent_points(c, p)：点 p から円 c に引いた接線の接点を返す．
• common_tangents(c1, c2)：二円の共通接線を列挙する．
• common_area(c1, c2)：二円の共通部分の面積を返す．

交点を返す関数には，交点が存在する入力を与える必要がある．

資料

• geometry/circle.hpp - Nyaan’s Library

Depends on

• 二次元幾何の基本要素 (geometry/geometry-base.hpp)
• 直線 (geometry/line.hpp)
• 多角形 (geometry/polygon.hpp)

Verified with

• verify/geometry/AOJ_CGL_7_A.test.cpp
• verify/geometry/AOJ_CGL_7_D.test.cpp
• verify/geometry/AOJ_CGL_7_E.test.cpp
• verify/geometry/AOJ_CGL_7_F.test.cpp
• verify/geometry/AOJ_CGL_7_G.test.cpp
• verify/geometry/AOJ_CGL_7_I.test.cpp

Code

#pragma once

#include "geometry/line.hpp"

struct Circle {
  Point p;
  Real r;

  Circle() = default;
  Circle(const Point& _p, Real _r) : p(_p), r(_r) {}
};

using Circles = vector<Circle>;

// 0: contained, 1: internally tangent, 2: intersecting,
// 3: externally tangent, 4: separate.
int intersect(Circle c1, Circle c2) {
  if (c1.r < c2.r) swap(c1, c2);
  Real d = abs(c1.p - c2.p);
  if (c1.r + c2.r < d - EPS) return 4;
  if (equals(c1.r + c2.r, d)) return 3;
  if (c1.r - c2.r < d - EPS) return 2;
  if (equals(c1.r - c2.r, d)) return 1;
  return 0;
}

pair<Point, Point> cross_point_cl(const Circle& c, const Line& l) {
  Point q = projection(l, c.p);
  Real h2 = max<Real>(0, c.r * c.r - norm(q - c.p));
  Point e = (l.b - l.a) / abs(l.b - l.a);
  Point d = e * sqrt(h2);
  return {q - d, q + d};
}

pair<Point, Point> cross_point_cc(const Circle& c1, const Circle& c2) {
  Point d = c2.p - c1.p;
  Real d2 = norm(d), len = sqrt(d2);
  Real x = (c1.r * c1.r - c2.r * c2.r + d2) / (2 * len);
  Real h = sqrt(max<Real>(0, c1.r * c1.r - x * x));
  Point q = c1.p + d * (x / len);
  Point v = d.rotate90() * (h / len);
  return {q - v, q + v};
}

// Compatibility with the reference implementation.
pair<Point, Point> crosspoint(const Circle& c1, const Circle& c2) {
  return cross_point_cc(c1, c2);
}

pair<Point, Point> tangent_points(const Circle& c, const Point& p) {
  Point d = p - c.p;
  Real d2 = norm(d);
  Real x = c.r * c.r / d2;
  Real y = c.r * sqrt(max<Real>(0, d2 - c.r * c.r)) / d2;
  Point q = c.p + d * x;
  Point v = d.rotate90() * y;
  return {q - v, q + v};
}

Lines common_tangents(const Circle& c1, const Circle& c2) {
  Lines ret;
  Point d = c2.p - c1.p;
  Real d2 = norm(d);
  if (equals(d2, 0)) return ret;
  for (int s : {-1, 1}) {
    Real r = c1.r - s * c2.r;
    Real h2 = d2 - r * r;
    if (h2 < -EPS) continue;
    h2 = max<Real>(0, h2);
    for (int t : {-1, 1}) {
      Point v = (d * r + d.rotate90() * (sqrt(h2) * t)) / d2;
      ret.emplace_back(c1.p + v * c1.r, c2.p + v * (s * c2.r));
      if (equals(h2, 0)) break;
    }
  }
  return ret;
}

Real common_area(const Circle& c1, const Circle& c2) {
  Real d = abs(c1.p - c2.p);
  if (d >= c1.r + c2.r - EPS) return 0;
  if (d <= abs(c1.r - c2.r) + EPS) {
    Real r = min(c1.r, c2.r);
    return PI * r * r;
  }
  Real a1 = acos(clamp((d * d + c1.r * c1.r - c2.r * c2.r) / (2 * d * c1.r), (Real)-1, (Real)1));
  Real a2 = acos(clamp((d * d + c2.r * c2.r - c1.r * c1.r) / (2 * d * c2.r), (Real)-1, (Real)1));
  auto x_minus_sin = [](Real x) {
    if (abs(x) > 0.1L) return x - sin(x);
    Real x2 = x * x;
    return x * x2 * (1.0L / 6 - x2 / 120 + x2 * x2 / 5040 - x2 * x2 * x2 / 362880 + x2 * x2 * x2 * x2 / 39916800);
  };
  return (c1.r * c1.r * x_minus_sin(2 * a1) + c2.r * c2.r * x_minus_sin(2 * a2)) / 2;
}

/**
 * @brief 円
 * @docs docs/geometry/circle.md
 */

#line 2 "geometry/circle.hpp"

#line 2 "geometry/line.hpp"

#line 2 "geometry/geometry-base.hpp"

#include <bits/stdc++.h>

using Real = long double;
constexpr Real EPS = 1e-10;
constexpr Real PI = 3.141592653589793238462643383279L;
bool equals(Real x, Real y) { return fabs(x - y) < EPS; }
int sign(Real a) { return equals(a, 0) ? 0 : (a > 0 ? 1 : -1); }

template <class R>
struct PointBase {
  using P = PointBase;
  R x, y;
  PointBase() : x(0), y(0) {}
  PointBase(R _x, R _y) : x(_x), y(_y) {}
  template <typename T, typename U>
  PointBase(const pair<T, U>& p) : x(p.first), y(p.second) {}

  P operator+(const P& r) const { return P{x + r.x, y + r.y}; }
  P operator-(const P& r) const { return P{x - r.x, y - r.y}; }
  P operator-() const { return P{-x, -y}; }
  P operator*(R r) const { return P{x * r, y * r}; }
  P operator/(R r) const { return P{x / r, y / r}; }

  P& operator+=(const P& r) { return (*this) = (*this) + r; }
  P& operator-=(const P& r) { return (*this) = (*this) - r; }
  P& operator*=(R r) { return (*this) = (*this) * r; }
  P& operator/=(R r) { return (*this) = (*this) / r; }

  bool operator<(const P& r) const { return x != r.x ? x < r.x : y < r.y; }
  bool operator==(const P& r) const { return x == r.x and y == r.y; }
  bool operator!=(const P& r) const { return !((*this) == r); }

  P rotate(R rad) const {
    return {x * cos(rad) - y * sin(rad), x * sin(rad) + y * cos(rad)};
  }
  P rotate90() const { return {-y, x}; }

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

  friend P operator*(R r, const P& p) { return p * r; }
  friend R real(const P& p) { return p.x; }
  friend R imag(const P& p) { return p.y; }
  friend R dot(const P& l, const P& r) { return l.x * r.x + l.y * r.y; }
  friend R cross(const P& l, const P& r) { return l.x * r.y - l.y * r.x; }
  friend R abs(const P& p) { return sqrt(p.x * p.x + p.y * p.y); }
  friend R norm(const P& p) { return p.x * p.x + p.y * p.y; }
  friend R arg(const P& p) { return atan2(p.y, p.x); }

  friend istream& operator>>(istream& is, P& p) {
    R a, b;
    is >> a >> b;
    p = P{a, b};
    return is;
  }
  friend ostream& operator<<(ostream& os, const P& p) {
    return os << p.x << " " << p.y;
  }
};
using Point = PointBase<Real>;
using Points = vector<Point>;

// relative position of c from a->b
int ccw(const Point& a, const Point& b, const Point& c) {
  Point x = b - a, y = c - a;
  if (cross(x, y) > EPS) return +1;        // counter-clockwise
  if (cross(x, y) < -EPS) return -1;       // clockwise
  if (dot(x, y) < -EPS) return +2;         // collinear in the order c-a-b
  if (norm(x) + EPS < norm(y)) return -2;  // collinear in the order a-b-c
  return 0;                                // collinear in the order a-c-b
}

/**
 * @brief 二次元幾何の基本要素
 * @docs docs/geometry/geometry-base.md
 */
#line 2 "geometry/polygon.hpp"

#line 4 "geometry/polygon.hpp"

// 反時計周り
using Polygon = vector<Point>;

// 多角形の内部に点があるか
// OUT : 0, ON : 1, IN : 2
// 凸でなくてもよい, 自己交差はもたない
// crossing number algorithm : x 軸正方向に伸ばした半直線が奇数回交われば内部
int contains_polygon(const Polygon& P, const Point& p) {
  bool in = false;
  for (int i = 0; i < (int)P.size(); i++) {
    Point a = P[i] - p, b = P[(i + 1) % P.size()] - p;
    if (a.y > b.y) swap(a, b);
    if (equals(cross(a, b), 0) && sign(dot(a, b)) <= 0) return 1;
    if (sign(a.y) <= 0 && 0 < sign(b.y) && sign(cross(a, b)) < 0)
      in = !in;
  }
  return in ? 2 : 0;
}

// 多角形の面積
Real area(const Polygon& P) {
  Real A = 0;
  for (int i = 0; i < (int)P.size(); i++) {
    A += cross(P[i], P[(i + 1) % P.size()]);
  }
  return A * 0.5;
}

bool is_convex(const Polygon& P) {
  int n = P.size();
  if (n < 3) return true;
  for (int i = 0; i < n; ++i) {
    if (ccw(P[i], P[(i + 1) % n], P[(i + 2) % n]) == -1) return false;
  }
  return true;
}

// 頂点集合から凸包を生成
// boundary : 周上の点も列挙する場合 true
template <bool boundary = false>
Polygon convex_hull(vector<Point> ps) {
  sort(begin(ps), end(ps));
  ps.erase(unique(begin(ps), end(ps)), end(ps));
  int n = ps.size(), k = 0;
  if (n <= 2) return ps;
  if constexpr (boundary) {
    bool collinear = true;
    for (int i = 2; i < n; ++i) {
      if (!equals(cross(ps[1] - ps[0], ps[i] - ps[0]), 0)) {
        collinear = false;
        break;
      }
    }
    if (collinear) return ps;
  }
  vector<Point> ch(2 * n);
  // 反時計周り
  const Real th = boundary ? -EPS : +EPS;
  for (int i = 0; i < n; ch[k++] = ps[i++]) {
    while (k >= 2 && cross(ch[k - 1] - ch[k - 2], ps[i] - ch[k - 1]) < th) --k;
  }
  for (int i = n - 2, t = k + 1; i >= 0; ch[k++] = ps[i--]) {
    while (k >= t && cross(ch[k - 1] - ch[k - 2], ps[i] - ch[k - 1]) < th) --k;
  }
  ch.resize(k - 1);
  return ch;
}

// 凸包の内部に点があるか
// OUT : 0, ON : 1, IN : 2
int contains_convex(const Polygon& C, const Point& p) {
  int N = C.size();
  if (N == 0) return 0;
  if (N == 1) return abs(C[0] - p) < EPS ? 1 : 0;
  if (N == 2) return ccw(C[0], C[1], p) == 0 ? 1 : 0;
  auto b1 = cross(C[1] - C[0], p - C[0]);
  auto b2 = cross(C[N - 1] - C[0], p - C[0]);
  if (b1 < -EPS or b2 > EPS) return 0;
  int L = 1, R = N - 1;
  while (L + 1 < R) {
    int M = (L + R) / 2;
    (cross(p - C[0], C[M] - C[0]) >= 0 ? R : L) = M;
  }
  auto v = cross(C[L] - p, C[R] - p);
  if (equals(v, 0)) {
    return 1;
  } else if (v > 0) {
    return equals(b1, 0) or equals(b2, 0) ? 1 : 2;
  } else {
    return 0;
  }
}

// 凸多角形の最遠点対を返す
// 返り値：頂点番号のペア
pair<int, int> convex_polygon_diameter(const Polygon& C) {
  int N = (int)C.size();
  assert(N > 0);
  if (N == 1) return {0, 0};
  if (N == 2) return {0, 1};
  int is = 0, js = 0;
  for (int i = 1; i < N; i++) {
    if (C[i].y > C[is].y) is = i;
    if (C[i].y < C[js].y) js = i;
  }
  Real maxdis = norm(C[is] - C[js]);

  int maxi, maxj, i, j;
  i = maxi = is;
  j = maxj = js;
  do {
    if (cross(C[(i + 1) % N] - C[i], C[(j + 1) % N] - C[j]) >= 0) {
      j = (j + 1) % N;
    } else {
      i = (i + 1) % N;
    }
    if (norm(C[i] - C[j]) > maxdis) {
      maxdis = norm(C[i] - C[j]);
      maxi = i;
      maxj = j;
    }
  } while (i != is || j != js);
  return minmax(maxi, maxj);
}

// min argmax_i dot(C[i], p)
int argmax_dot(const Polygon& C, const Point& p) {
  int n = C.size();
  assert(n > 0);
  assert(p.x != 0 || p.y != 0);
  if (n <= 2) {
    if (n == 1 || dot(C[0], p) >= dot(C[1], p)) return 0;
    return 1;
  }
  auto arg_half = [](const Point& a) {
    return a.y < 0 || (a.y == 0 && a.x < 0);
  };
  auto arg_less = [&](const Point& a, const Point& b) {
    int ha = arg_half(a), hb = arg_half(b);
    return ha != hb ? ha < hb : cross(a, b) > 0;
  };
  auto edge = [&](int i) { return C[(i + 1) % n] - C[i]; };
  int l = 0, r = n - 1;
  while (l < r) {
    int m = (l + r) / 2;
    if (arg_less(edge(r), edge(m)))
      l = m + 1;
    else
      r = m;
  }
  int s = l;
  Point q = p.rotate90();
  l = 0, r = n;
  while (l < r) {
    int m = (l + r) / 2;
    if (arg_less(edge((s + m) % n), q))
      l = m + 1;
    else
      r = m;
  }
  int i = (s + (l == n ? 0 : l)) % n;
  Point e = edge(i);
  if (!arg_less(e, q) && !arg_less(q, e)) {
    int low = l;
    l = low, r = n;
    while (l < r) {
      int m = (l + r) / 2;
      if (!arg_less(q, edge((s + m) % n)))
        l = m + 1;
      else
        r = m;
    }
    int count = l - low;
    return i + count >= n ? 0 : i;
  }
  return i;
}

struct ConvexPolygonCutResult {
  int first = 0, last = 0;
  optional<Point> front, back;

  bool empty() const { return first == last && !front && !back; }
  Polygon to_polygon(const Polygon& C) const {
    Polygon ret;
    auto push = [&](const Point& p) { if (ret.empty() || abs(ret.back() - p) >= EPS) ret.push_back(p); };
    if (front) push(*front);
    for (int i = first; i < last; ++i) push(C[i % C.size()]);
    if (back) push(*back);
    if (ret.size() >= 2 && abs(ret.front() - ret.back()) < EPS) ret.pop_back();
    return ret;
  }
};

ConvexPolygonCutResult convex_polygon_cut_info(const Polygon& C, const Point& a, const Point& b) {
  int n = C.size();
  assert(n > 0 && a != b);
  Point d = b - a, normal = d.rotate90();
  auto side = [&](int i) { return sign(cross(d, C[i % n] - a)); };
  int mx = argmax_dot(C, normal), mn = argmax_dot(C, -normal);
  if (side(mx) < 0) return {};
  if (side(mn) >= 0) return {0, n, nullopt, nullopt};

  int l = mn, r = mx;
  if (r < l) r += n;
  while (l + 1 < r) {
    int m = (l + r) / 2;
    (side(m) >= 0 ? r : l) = m;
  }
  int first = r;
  l = mx, r = mn;
  if (r <= l) r += n;
  while (l + 1 < r) {
    int m = (l + r) / 2;
    (side(m) < 0 ? r : l) = m;
  }
  int last = r;
  if (last < first) last += n;
  int count = last - first;
  first %= n;
  ConvexPolygonCutResult ret{first, first + count, nullopt, nullopt};
  auto cross_point = [&](const Point& p, const Point& q) {
    return p + (q - p) * (cross(d, a - p) / cross(d, q - p));
  };
  if (side(first) > 0) ret.front = cross_point(C[(first + n - 1) % n], C[first]);
  int back = (ret.last - 1) % n;
  if (side(back) > 0) ret.back = cross_point(C[back], C[(back + 1) % n]);
  return ret;
}

Polygon convex_polygon_cut(const Polygon& C, const Point& a, const Point& b) {
  return convex_polygon_cut_info(C, a, b).to_polygon(C);
}

/**
 * @brief 多角形
 * @docs docs/geometry/polygon.md
 */
#line 5 "geometry/line.hpp"

struct Line {
  using P = Point;
  P a, b;

  Line() = default;
  Line(const P& _a, const P& _b) : a(_a), b(_b) {}

  // Ax + By = C
  Line(Real A, Real B, Real C) {
    assert(!equals(A, 0) || !equals(B, 0));
    if (!equals(B, 0)) {
      a = P{0, C / B};
      b = P{1, (C - A) / B};
    } else {
      a = P{C / A, 0};
      b = P{C / A, 1};
    }
  }

  friend istream& operator>>(istream& is, Line& l) { return is >> l.a >> l.b; }
  friend ostream& operator<<(ostream& os, const Line& l) {
    return os << l.a << " to " << l.b;
  }
};

using Lines = vector<Line>;

bool is_parallel(const Line& l, const Line& m) {
  return equals(cross(l.b - l.a, m.b - m.a), 0);
}

bool is_orthogonal(const Line& l, const Line& m) {
  return equals(dot(l.b - l.a, m.b - m.a), 0);
}

bool is_intersect_ll(const Line& l, const Line& m) {
  return !is_parallel(l, m) || equals(cross(l.b - l.a, m.a - l.a), 0);
}

// Parallel distinct lines must not be given.
Point cross_point_ll(const Line& l, const Line& m) {
  Point u = l.b - l.a, v = m.b - m.a;
  Real d = cross(u, v);
  if (equals(d, 0)) {
    assert(equals(cross(u, m.a - l.a), 0));
    return m.a;
  }
  return l.a + u * (cross(m.a - l.a, v) / d);
}

Point projection(const Line& l, const Point& p) {
  Point d = l.b - l.a;
  return l.a + d * (dot(p - l.a, d) / norm(d));
}

Point reflection(const Line& l, const Point& p) {
  return projection(l, p) * 2 - p;
}

Real distance_lp(const Line& l, const Point& p) {
  return abs(p - projection(l, p));
}

// Returns the left side of the directed line l.a -> l.b.
Polygon convex_polygon_cut(const Polygon& P, const Line& l) {
  return convex_polygon_cut(P, l.a, l.b);
}

ConvexPolygonCutResult convex_polygon_cut_info(const Polygon& P,
                                               const Line& l) {
  return convex_polygon_cut_info(P, l.a, l.b);
}

/**
 * @brief 直線
 * @docs docs/geometry/line.md
 */
#line 4 "geometry/circle.hpp"

struct Circle {
  Point p;
  Real r;

  Circle() = default;
  Circle(const Point& _p, Real _r) : p(_p), r(_r) {}
};

using Circles = vector<Circle>;

// 0: contained, 1: internally tangent, 2: intersecting,
// 3: externally tangent, 4: separate.
int intersect(Circle c1, Circle c2) {
  if (c1.r < c2.r) swap(c1, c2);
  Real d = abs(c1.p - c2.p);
  if (c1.r + c2.r < d - EPS) return 4;
  if (equals(c1.r + c2.r, d)) return 3;
  if (c1.r - c2.r < d - EPS) return 2;
  if (equals(c1.r - c2.r, d)) return 1;
  return 0;
}

pair<Point, Point> cross_point_cl(const Circle& c, const Line& l) {
  Point q = projection(l, c.p);
  Real h2 = max<Real>(0, c.r * c.r - norm(q - c.p));
  Point e = (l.b - l.a) / abs(l.b - l.a);
  Point d = e * sqrt(h2);
  return {q - d, q + d};
}

pair<Point, Point> cross_point_cc(const Circle& c1, const Circle& c2) {
  Point d = c2.p - c1.p;
  Real d2 = norm(d), len = sqrt(d2);
  Real x = (c1.r * c1.r - c2.r * c2.r + d2) / (2 * len);
  Real h = sqrt(max<Real>(0, c1.r * c1.r - x * x));
  Point q = c1.p + d * (x / len);
  Point v = d.rotate90() * (h / len);
  return {q - v, q + v};
}

// Compatibility with the reference implementation.
pair<Point, Point> crosspoint(const Circle& c1, const Circle& c2) {
  return cross_point_cc(c1, c2);
}

pair<Point, Point> tangent_points(const Circle& c, const Point& p) {
  Point d = p - c.p;
  Real d2 = norm(d);
  Real x = c.r * c.r / d2;
  Real y = c.r * sqrt(max<Real>(0, d2 - c.r * c.r)) / d2;
  Point q = c.p + d * x;
  Point v = d.rotate90() * y;
  return {q - v, q + v};
}

Lines common_tangents(const Circle& c1, const Circle& c2) {
  Lines ret;
  Point d = c2.p - c1.p;
  Real d2 = norm(d);
  if (equals(d2, 0)) return ret;
  for (int s : {-1, 1}) {
    Real r = c1.r - s * c2.r;
    Real h2 = d2 - r * r;
    if (h2 < -EPS) continue;
    h2 = max<Real>(0, h2);
    for (int t : {-1, 1}) {
      Point v = (d * r + d.rotate90() * (sqrt(h2) * t)) / d2;
      ret.emplace_back(c1.p + v * c1.r, c2.p + v * (s * c2.r));
      if (equals(h2, 0)) break;
    }
  }
  return ret;
}

Real common_area(const Circle& c1, const Circle& c2) {
  Real d = abs(c1.p - c2.p);
  if (d >= c1.r + c2.r - EPS) return 0;
  if (d <= abs(c1.r - c2.r) + EPS) {
    Real r = min(c1.r, c2.r);
    return PI * r * r;
  }
  Real a1 = acos(clamp((d * d + c1.r * c1.r - c2.r * c2.r) / (2 * d * c1.r), (Real)-1, (Real)1));
  Real a2 = acos(clamp((d * d + c2.r * c2.r - c1.r * c1.r) / (2 * d * c2.r), (Real)-1, (Real)1));
  auto x_minus_sin = [](Real x) {
    if (abs(x) > 0.1L) return x - sin(x);
    Real x2 = x * x;
    return x * x2 * (1.0L / 6 - x2 / 120 + x2 * x2 / 5040 - x2 * x2 * x2 / 362880 + x2 * x2 * x2 * x2 / 39916800);
  };
  return (c1.r * c1.r * x_minus_sin(2 * a1) + c2.r * c2.r * x_minus_sin(2 * a2)) / 2;
}

/**
 * @brief 円
 * @docs docs/geometry/circle.md
 */

Back to top page