This documentation is automatically generated by online-judge-tools/verification-helper
View the Project on GitHub ruthen71/rcpl
#define PROBLEM "http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_G" #define ERROR 0.00001 #include <iostream> #include <iomanip> #include <algorithm> #include "geometry/tangent.hpp" int main() { Circle<double> C1, C2; std::cin >> C1 >> C2; auto lines = tangent(C1, C2); std::vector<Point<double>> ans; for (auto&& l : lines) { if (is_intersect(C1, l.a)) { ans.push_back(l.a); } if (is_intersect(C1, l.b)) { ans.push_back(l.b); } } std::sort(ans.begin(), ans.end(), compare_x<double>); for (auto&& p : ans) { std::cout << std::fixed << std::setprecision(15) << p.x << ' ' << p.y << '\n'; } return 0; }
#line 1 "verify/geometry/tangent_cc.test.cpp" #define PROBLEM "http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_G" #define ERROR 0.00001 #include <iostream> #include <iomanip> #include <algorithm> #line 2 "geometry/tangent.hpp" #line 2 "geometry/point.hpp" #line 2 "geometry/geometry_template.hpp" #include <type_traits> // Constants (EPS, PI) // EPS の変更は Constants<T>::set_eps(new_eps) で template <class T> struct Constants { static T EPS; static void set_eps(const T e) { EPS = e; } static constexpr T PI = 3.14159'26535'89793L; }; template <> double Constants<double>::EPS = 1e-9; template <> long double Constants<long double>::EPS = 1e-12; template <> long long Constants<long long>::EPS = 0; // 汎用関数 template <class T> inline int sign(const T x) { return x < -Constants<T>::EPS ? -1 : (x > Constants<T>::EPS ? 1 : 0); } template <class T> inline bool equal(const T a, const T b) { return sign(a - b) == 0; } template <class T> inline T radian_to_degree(const T r) { return r * 180.0 / Constants<T>::PI; } template <class T> inline T degree_to_radian(const T d) { return d * Constants<T>::PI / 180.0; } // type traits template <class T> using is_geometry_floating_point = typename std::conditional<std::is_same<T, double>::value || std::is_same<T, long double>::value, std::true_type, std::false_type>::type; template <class T> using is_geometry_integer = typename std::conditional<std::is_same<T, long long>::value, std::true_type, std::false_type>::type; template <class T> using is_geometry = typename std::conditional<is_geometry_floating_point<T>::value || is_geometry_integer<T>::value, std::true_type, std::false_type>::type; #line 4 "geometry/point.hpp" #include <cmath> #include <cassert> // 点 template <class T> struct Point { T x, y; Point() = default; Point(const T x, const T y) : x(x), y(y) {} template <class U> Point(const Point<U> p) : x(p.x), y(p.y) {} 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 Point& p) { static_assert(is_geometry_floating_point<T>::value == true); return *this = Point(x * p.x - y * p.y, x * p.y + y * p.x); } Point& operator/=(const Point& p) { static_assert(is_geometry_floating_point<T>::value == true); return *this = Point(x * p.x + y * p.y, -x * p.y + y * p.x) / (p.x * p.x + p.y * p.y); } Point& operator*=(const T k) { x *= k, y *= k; return *this; } Point& operator/=(const T k) { static_assert(is_geometry_floating_point<T>::value == true); x /= k, y /= k; return *this; } Point operator+() const { return *this; } Point operator-() const { return Point(-x, -y); } friend Point operator+(const Point& a, const Point& b) { return Point(a) += b; } friend Point operator-(const Point& a, const Point& b) { return Point(a) -= b; } friend Point operator*(const Point& a, const Point& b) { return Point(a) *= b; } friend Point operator/(const Point& a, const Point& b) { return Point(a) /= b; } friend Point operator*(const Point& p, const T k) { return Point(p) *= k; } friend Point operator/(const Point& p, const T k) { return Point(p) /= k; } // 辞書式順序 friend bool operator<(const Point& a, const Point& b) { return a.x == b.x ? a.y < b.y : a.x < b.x; } friend bool operator>(const Point& a, const Point& b) { return a.x == b.x ? a.y > b.y : a.x > b.x; } friend bool operator==(const Point& a, const Point& b) { return a.x == b.x and a.y == b.y; } // I/O friend std::istream& operator>>(std::istream& is, Point& p) { return is >> p.x >> p.y; } friend std::ostream& operator<<(std::ostream& os, const Point& p) { return os << '(' << p.x << ' ' << p.y << ')'; } }; // 汎用関数 // 点の一致判定 template <class T> inline bool equal(const Point<T>& a, const Point<T>& b) { return equal(a.x, b.x) and equal(a.y, b.y); } // 内積 template <class T> inline T dot(const Point<T>& a, const Point<T>& b) { return a.x * b.x + a.y * b.y; } // 外積 template <class T> inline T cross(const Point<T>& a, const Point<T>& b) { return a.x * b.y - a.y * b.x; } // rad ラジアンだけ反時計回りに回転 template <class T> inline Point<T> rotate(const Point<T>& p, const T theta) { static_assert(is_geometry_floating_point<T>::value == true); return p * Point<T>(std::cos(theta), std::sin(theta)); } // (x, y) の辞書式順序 (誤差許容) template <class T> inline bool compare_x(const Point<T>& a, const Point<T>& b) { return equal(a.x, b.x) ? sign(a.y - b.y) < 0 : sign(a.x - b.x) < 0; } // (y, x) の辞書式順序 (誤差許容) template <class T> inline bool compare_y(const Point<T>& a, const Point<T>& b) { return equal(a.y, b.y) ? sign(a.x - b.x) < 0 : sign(a.y - b.y) < 0; } // 整数のまま行う偏角ソート // 無限の精度をもつ arg(p) = atan2(y, x) で比較し, 同じ場合は norm(p) で比較 (atan2(0, 0) = 0 とする) // 基本的に (-PI, PI] でソートされ, 点 (0, 0) は (-PI, 0) と [0, PI] の間に入る // https://ngtkana.hatenablog.com/entry/2021/11/13/202103 // https://judge.yosupo.jp/problem/sort_points_by_argument template <class T> inline bool compare_atan2(const Point<T>& a, const Point<T>& b) { static_assert(is_geometry_integer<T>::value == true); if ((Point<T>(a.y, -a.x) > Point<T>(0, 0)) == (Point<T>(b.y, -b.x) > Point<T>(0, 0))) { // a, b in (-PI, 0] or a, b in (0, PI] if (a.x * b.y != a.y * b.x) return a.x * b.y > a.y * b.x; // cross(a, b) != 0 return a == Point<T>(0, 0) ? b.x > 0 and b.y == 0 : (b == Point<T>(0, 0) ? a.y < 0 : norm(a) < norm(b)); } return Point<T>(a.y, -a.x) < Point<T>(b.y, -b.x); } // 絶対値の 2 乗 template <class T> inline T norm(const Point<T>& p) { return p.x * p.x + p.y * p.y; } // 絶対値 template <class T> inline T abs(const Point<T>& p) { static_assert(is_geometry_floating_point<T>::value == true); return std::sqrt(norm(p)); } // 偏角 template <class T> inline T arg(const Point<T>& p) { static_assert(is_geometry_floating_point<T>::value == true); return std::atan2(p.y, p.x); // (-PI, PI] } // 共役複素数 (x 軸について対象な点) template <class T> inline Point<T> conj(const Point<T>& p) { return Point(p.x, -p.y); } // 極座標系 -> 直交座標系 (絶対値が rho で偏角が theta ラジアン) template <class T> inline Point<T> polar(const T rho, const T theta = T(0)) { static_assert(is_geometry_floating_point<T>::value == true); assert(sign(rho) >= 0); return Point<T>(std::cos(theta), std::sin(theta)) * rho; } // ccw の戻り値 enum class Ccw { COUNTER_CLOCKWISE = 1, // a->b->c 反時計回り CLOCKWISE = -1, // a->b->c 時計回り ONLINE_BACK = 2, // c->a->b 直線 ONLINE_FRONT = -2, // a->b->c 直線 ON_SEGMENT = 0 // a->c->b 直線 }; // 点 a, b, c の位置関係 // http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_1_C template <class T> Ccw ccw(const Point<T>& a, const Point<T>& b, const Point<T>& c) { Point<T> ab = b - a; Point<T> ac = c - a; if (sign(cross(ab, ac)) == 1) return Ccw::COUNTER_CLOCKWISE; if (sign(cross(ab, ac)) == -1) return Ccw::CLOCKWISE; if (sign(dot(ab, ac)) == -1) return Ccw::ONLINE_BACK; if (sign(norm(ab) - norm(ac)) == -1) return Ccw::ONLINE_FRONT; return Ccw::ON_SEGMENT; } // 線分 a -> b から 線分 a -> c までの角度 (ラジアンで -PI より大きく PI 以下) template <class T> T get_angle(const Point<T>& a, const Point<T>& b, const Point<T>& c) { Point<T> ab = b - a; Point<T> ac = c - a; // a-bが x 軸になるように回転 ac *= conj(ab) / norm(ab); return arg(ac); // (-PI, PI] } #line 2 "geometry/circle.hpp" #line 4 "geometry/circle.hpp" // 円 template <class T> struct Circle { Point<T> o; T r; Circle() = default; Circle(const Point<T>& o, const T r) : o(o), r(r) {} friend std::istream& operator>>(std::istream& is, Circle& c) { return is >> c.o >> c.r; } friend std::ostream& operator<<(std::ostream& os, const Circle& c) { return os << c.o << ", " << c.r; } }; // 共通接線の本数 // http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_A template <class T> int tangent_number(Circle<T> c1, Circle<T> c2) { if (c1.r < c2.r) std::swap(c1, c2); const T d2 = norm(c1.o - c2.o); if (sign(d2 - (c1.r + c2.r) * (c1.r + c2.r)) == 1) return 4; // d > c1.r + c2.r and c1.r + c2.r >= 0 <=> d ^ 2 > (c1.r + c2.r) ^ 2 if (sign(d2 - (c1.r + c2.r) * (c1.r + c2.r)) == 0) return 3; // d = c1.r + c2.r and c1.r + c2.r >= 0 <=> d ^ 2 = (c1.r + c2.r) ^ 2 if (sign(d2 - (c1.r - c2.r) * (c1.r - c2.r)) == 1) return 2; // d > c1.r - c2.r and c1.r - c2.r >= 0 <=> d ^ 2 > (c1.r - c2.r) ^ 2 if (sign(d2 - (c1.r - c2.r) * (c1.r - c2.r)) == 0) return 1; // d = c1.r - c2.r and c1.r - c2.r >= 0 <=> d ^ 2 = (c1.r - c2.r) ^ 2 return 0; } #line 2 "geometry/cross_point.hpp" #line 2 "geometry/line.hpp" #line 4 "geometry/line.hpp" // 直線 template <class T> struct Line { Point<T> a, b; Line() = default; Line(const Point<T>& a, const Point<T>& b) : a(a), b(b) {} // Ax + By = C Line(const T A, const T B, const T C) { static_assert(is_geometry_floating_point<T>::value == true); assert(!(equal(A, T(0)) and equal(B, T(0)))); if (equal(A, T(0))) { a = Point<T>(T(0), C / B), b = Point<T>(T(1), C / B); } else if (equal(B, T(0))) { a = Point<T>(C / A, T(0)), b = Point<T>(C / A, T(1)); } else if (equal(C, T(0))) { a = Point<T>(T(0), T(0)), b = Point<T>(T(1), -A / B); } else { a = Point<T>(T(0), C / B), b = Point<T>(C / A, T(0)); } } friend std::istream& operator>>(std::istream& is, Line& p) { return is >> p.a >> p.b; } friend std::ostream& operator<<(std::ostream& os, const Line& p) { return os << p.a << "->" << p.b; } }; // 線分 template <class T> struct Segment : Line<T> { Segment() = default; Segment(const Point<T>& a, const Point<T>& b) : Line<T>(a, b) {} }; // 点 p から直線 l に下ろした垂線と直線 l の交点 // http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_1_A template <class T> Point<T> projection(const Line<T>& l, const Point<T>& p) { static_assert(is_geometry_floating_point<T>::value == true); T t = dot(p - l.a, l.b - l.a) / norm(l.b - l.a); return l.a + (l.b - l.a) * t; } // 直線 l に関して点 p と対象な位置にある点 // http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_1_B template <class T> Point<T> reflection(const Line<T>& l, const Point<T>& p) { static_assert(is_geometry_floating_point<T>::value == true); return p + (projection(l, p) - p) * T(2); } // 直線 l1, l2 の垂直判定 // http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_2_A template <class T> inline bool is_orthogonal(const Line<T>& l1, const Line<T>& l2) { return sign(dot(l1.b - l1.a, l2.b - l2.a)) == 0; } // 直線 l1, l2 の平行判定 // http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_2_A template <class T> inline bool is_parallel(const Line<T>& l1, const Line<T>& l2) { return sign(cross(l1.b - l1.a, l2.b - l2.a)) == 0; } #line 2 "geometry/is_intersect.hpp" #line 6 "geometry/is_intersect.hpp" // 交差判定 (直線, 線分, 円, 点) // 直線 l1, l2 の交差判定 template <class T> bool is_intersect(const Line<T>& l1, const Line<T>& l2) { Point<T> base = l1.b - l1.a; T d12 = cross(base, l2.b - l2.a); T d1 = cross(base, l1.b - l2.a); // sign(d12) != 0 -> 平行でないので交差する // sign(d12) == 0 and sign(d1) == 0 -> 一致するので交差する return sign(d12) != 0 or sign(d1) == 0; } // 直線 l, 点 p の交差判定 template <class T> inline bool is_intersect(const Line<T>& l, const Point<T>& p) { auto res = ccw(l.a, l.b, p); return res == Ccw::ONLINE_BACK or res == Ccw::ONLINE_FRONT or res == Ccw::ON_SEGMENT; } template <class T> bool is_intersect(const Point<T>& p, const Line<T>& l) { return is_intersect(l, p); } // 線分 s, 点 p の交差判定 template <class T> inline bool is_intersect(const Segment<T>& s, const Point<T>& p) { return ccw(s.a, s.b, p) == Ccw::ON_SEGMENT; } template <class T> inline bool is_intersect(const Point<T>& p, const Segment<T>& s) { return ccw(s.a, s.b, p) == Ccw::ON_SEGMENT; } // 直線 l, 線分 s の交差判定 template <class T> bool is_intersect(const Line<T>& l, const Segment<T>& s) { // s.a と s.b が直線 l に関して同じ側 (直線上を除き) にある場合に限り交差しない auto c1 = ccw(l.a, l.b, s.a); auto c2 = ccw(l.a, l.b, s.b); return !((c1 == c2) and (c1 == Ccw::CLOCKWISE or c1 == Ccw::COUNTER_CLOCKWISE)); } template <class T> bool is_intersect(const Segment<T>& s, const Line<T>& l) { return is_intersect(l, s); } // 線分 s1, s2 の交差判定 // http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_2_B template <class T> inline bool is_intersect(const Segment<T>& s1, const Segment<T>& s2) { auto c1 = ccw(s1.a, s1.b, s2.a); auto c2 = ccw(s1.a, s1.b, s2.b); auto c3 = ccw(s2.a, s2.b, s1.a); auto c4 = ccw(s2.a, s2.b, s1.b); // 平行な場合: 重なるなら 1 次元で考えると必ず端点のどれかがもう一方の線分上にある if (c1 == Ccw::ON_SEGMENT or c2 == Ccw::ON_SEGMENT or c3 == Ccw::ON_SEGMENT or c4 == Ccw::ON_SEGMENT) return true; // 平行でない場合: 一方の線分の両側にもう一方の線分の端点がある, を線分を入れ替えても成立 // 平行だが重ならない場合は以下の条件は成立しない bool ok1 = ((c1 == Ccw::CLOCKWISE and c2 == Ccw::COUNTER_CLOCKWISE) or (c1 == Ccw::COUNTER_CLOCKWISE and c2 == Ccw::CLOCKWISE)); bool ok2 = ((c3 == Ccw::CLOCKWISE and c4 == Ccw::COUNTER_CLOCKWISE) or (c3 == Ccw::COUNTER_CLOCKWISE and c4 == Ccw::CLOCKWISE)); return ok1 and ok2; } // 点 p1, p2 の交差判定 template <class T> inline bool is_intersect(const Point<T>& p1, const Point<T>& p2) { return equal(p1, p2); } // 円 c1, c2 の交差判定 template <class T> inline bool is_intersect(const Circle<T>& c1, const Circle<T>& c2) { const int num = tangent_number(c1, c2); return 1 <= num and num <= 3; } // 円 c, 点 p の交差判定 template <class T> inline bool is_intersect(const Circle<T>& c, const Point<T>& p) { return equal(norm(p - c.o), c.r * c.r); } template <class T> inline bool is_intersect(const Point<T>& p, const Circle<T>& c) { return equal(norm(p - c.o), c.r * c.r); } // 円 c, 直線 l の交差判定 template <class T> inline bool is_intersect(const Circle<T>& c, const Line<T>& l) { static_assert(is_geometry_floating_point<T>::value == true); // norm(c.o - projection(l, c.o)) が直線 l と 円 c の中心 c.o の距離の 2 乗 return sign(c.r * c.r - norm(c.o - projection(l, c.o))) >= 0; } template <class T> inline bool is_intersect(const Line<T>& l, const Circle<T>& c) { return is_intersect(c, l); } // 円 c, 線分 s の交差判定 template <class T> inline bool is_intersect(const Circle<T>& c, const Segment<T>& s) { static_assert(is_geometry_floating_point<T>::value == true); if (!is_intersect(c, Line(s.a, s.b))) return false; // 直線としても交差しない T d1 = abs(c.o - s.a), d2 = abs(c.o - s.b); if (sign(d1 - c.r) == -1 and sign(d2 - c.r) == -1) return false; // 端点がどちらも円の内側 if (sign(d1 - c.r) * sign(d2 - c.r) <= 0) return true; // 円周をまたいでいる const Point<T> h = projection(s, c.o); return ccw(s.a, s.b, h) == Ccw::ON_SEGMENT; // s.a -> h -> s.b の順で並んでいる } template <class T> inline bool is_intersect(const Segment<T>& s, const Circle<T>& c) { return is_intersect(c, s); } #line 7 "geometry/cross_point.hpp" #include <vector> #line 10 "geometry/cross_point.hpp" // 交点 (直線, 線分, 円) // 交点を持たない場合の挙動は未定義 // 直線 l1, l2 の交点 1 つ template <class T> Point<T> cross_point(const Line<T>& l1, const Line<T>& l2) { static_assert(is_geometry_floating_point<T>::value == true); Point<T> base = l1.b - l1.a; T d12 = cross(base, l2.b - l2.a); T d1 = cross(base, l1.b - l2.a); if (sign(d12) == 0) { assert(sign(d1) == 0); // 平行かつ一致しない return l2.a; } return l2.a + (l2.b - l2.a) * (d1 / d12); } // 線分 s1, s2 の交点 1 つ // http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_2_C template <class T> Point<T> cross_point(const Segment<T>& s1, const Segment<T>& s2) { static_assert(is_geometry_floating_point<T>::value == true); Point<T> base = s1.b - s1.a; T d12 = cross(base, s2.b - s2.a); T d1 = cross(base, s1.b - s2.a); if (sign(d12) == 0) { assert(sign(d1) == 0); // 平行かつ一致しない // 端点のどれかを返す if (is_intersect(s1, s2.a)) return s2.a; if (is_intersect(s1, s2.b)) return s2.b; if (is_intersect(s2, s1.a)) return s1.a; assert(is_intersect(s2, s1.b)); return s1.b; } return s2.a + (s2.b - s2.a) * (d1 / d12); } // 直線 l, 線分 s の交点 1 つ template <class T> Point<T> cross_point(const Line<T>& l, const Segment<T>& s) { // cross_point(l1, l2) は重なるとき l2.a を返すので OK return cross_point(l, Line<T>(s.a, s.b)); } template <class T> Point<T> cross_point(const Segment<T>& s, const Line<T>& l) { return cross_point(l, s); } // 円 c, 直線 l の交点 1 or 2 つ // http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_D template <class T> std::vector<Point<T>> cross_point(const Circle<T>& c, const Line<T>& l) { static_assert(is_geometry_floating_point<T>::value == true); auto pr = projection(l, c.o); if (equal(norm(pr - c.o), c.r * c.r)) return {pr}; Point<T> e = (l.b - l.a) / abs(l.b - l.a); auto k = sqrtl(c.r * c.r - norm(pr - c.o)); return {pr - e * k, pr + e * k}; } template <class T> std::vector<Point<T>> cross_point(const Line<T>& l, const Circle<T>& c) { return cross_point(c, l); } // 円 c, 線分 s の交点 1 or 2 つ template <class T> std::vector<Point<T>> cross_point(const Circle<T>& c, const Segment<T>& s) { assert(is_intersect(c, s)); auto ps = cross_point(c, Line<T>(s.a, s.b)); std::vector<Point<T>> res; for (auto&& p : ps) if (ccw(s.a, s.b, p) == Ccw::ON_SEGMENT) res.emplace_back(p); return res; } template <class T> std::vector<Point<T>> cross_point(const Segment<T>& s, const Circle<T>& c) { return cross_point(c, s); } // 円 c1, c2 の交点 1 or 2 つ // http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_E template <class T> std::vector<Point<T>> cross_point(const Circle<T>& c1, const Circle<T>& c2) { static_assert(is_geometry_floating_point<T>::value == true); assert(is_intersect(c1, c2)); T d = abs(c1.o - c2.o); T a = std::acos((c1.r * c1.r - c2.r * c2.r + d * d) / (T(2) * c1.r * d)); T t = arg(c2.o - c1.o); Point<T> p = c1.o + polar(c1.r, t + a); Point<T> q = c1.o + polar(c1.r, t - a); if (equal(p, q)) return {p}; return {p, q}; } #line 6 "geometry/tangent.hpp" #line 8 "geometry/tangent.hpp" // 円 c の点 p に対する接線 // http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_F template <class T> std::vector<Line<T>> tangent(const Circle<T>& c, const Point<T>& p) { static_assert(is_geometry_floating_point<T>::value == true); assert(sign(norm(c.o - p) - c.r * c.r) == 1); auto tangent_points = cross_point(c, Circle(p, std::sqrt(norm(c.o - p) - c.r * c.r))); std::vector<Line<T>> res; for (auto&& tp : tangent_points) { res.emplace_back(p, tp); } return res; } // 円 c1, c2 の共通接線 // http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_G template <class T> std::vector<Line<T>> tangent(Circle<T> c1, Circle<T> c2) { assert(!(equal(c1.r, c2.r) and equal(c1.o, c2.o))); // 2 つの円は同一なので接線が無数にある std::vector<Line<T>> res; if (c1.r < c2.r) std::swap(c1, c2); // c1.r >= c2.r T g2 = norm(c1.o - c2.o); T g = abs(c1.o - c2.o); if (equal(g2, T(0))) return res; // 中心は一致しているが半径が異なる Point<T> u = (c2.o - c1.o) / g; // c1.o から c2.o への単位ベクトル Point<T> v = rotate(u, Constants<T>::PI / 2); // u と直行する単位ベクトル // 中心間距離と半径の和/差の比を用いて, 各円の中心から接点へのベクトルを求める for (int s : {-1, 1}) { T h = (c1.r + s * c2.r) / g; if (equal(1 - h * h, T(0))) { // 円が内接/外接する場合 -> 接線は 1 本 res.emplace_back(c1.o + u * c1.r, c1.o + (u + v) * c1.r); } else if (sign(1 - h * h) == 1) { // uu + vv, uu - vv が中心から接点方向への単位ベクトル Point<T> uu = u * h, vv = v * std::sqrt(1 - h * h); res.emplace_back(c1.o + (uu + vv) * c1.r, c2.o - (uu + vv) * c2.r * s); res.emplace_back(c1.o + (uu - vv) * c1.r, c2.o - (uu - vv) * c2.r * s); } } return res; } #line 9 "verify/geometry/tangent_cc.test.cpp" int main() { Circle<double> C1, C2; std::cin >> C1 >> C2; auto lines = tangent(C1, C2); std::vector<Point<double>> ans; for (auto&& l : lines) { if (is_intersect(C1, l.a)) { ans.push_back(l.a); } if (is_intersect(C1, l.b)) { ans.push_back(l.b); } } std::sort(ans.begin(), ans.end(), compare_x<double>); for (auto&& p : ans) { std::cout << std::fixed << std::setprecision(15) << p.x << ' ' << p.y << '\n'; } return 0; }