This documentation is automatically generated by online-judge-tools/verification-helper
View the Project on GitHub ebi-fly13/Library
#define PROBLEM \ "https://onlinejudge.u-aizu.ac.jp/courses/library/4/CGL/7/CGL_7_F" #define ERROR 0.00000001 #include <algorithm> #include <cassert> #include <cstdint> #include <iomanip> #include <iostream> #include <vector> #include "geometry/circle.hpp" namespace ebi { using i64 = std::int64_t; void main_() { circle c; point p; std::cin >> p; std::cin >> c.c >> c.r; auto ps = tangent_to_circle(c, p); std::sort(ps.begin(), ps.end()); for (auto p : ps) { std::cout << p << '\n'; } } } // namespace ebi int main() { std::cout << std::fixed << std::setprecision(15); std::cin.tie(nullptr); std::ios::sync_with_stdio(false); ebi::main_(); }
#line 1 "test/geometry/tangent_to_circle.test.cpp" #define PROBLEM \ "https://onlinejudge.u-aizu.ac.jp/courses/library/4/CGL/7/CGL_7_F" #define ERROR 0.00000001 #include <algorithm> #include <cassert> #include <cstdint> #include <iomanip> #include <iostream> #include <vector> #line 2 "geometry/circle.hpp" #include <cmath> #line 5 "geometry/circle.hpp" #line 2 "geometry/line.hpp" #line 5 "geometry/line.hpp" #line 2 "geometry/point.hpp" #line 9 "geometry/point.hpp" namespace ebi { constexpr long double EPS = 1e-10; const long double PI = std::acos(-1); namespace internal { int sgn(long double a) { return (a < -EPS) ? -1 : (a > EPS) ? 1 : 0; } long double add(long double a, long double b) { if (std::abs(a + b) < EPS * (std::abs(a) + std::abs(b))) return 0; return a + b; } } // namespace internal long double arg_to_radian(long double arg) { return PI * arg / (long double)(180); } struct point { long double x, y; point() = default; point(long double x, long double y) : x(x), y(y) {} point &operator+=(const point rhs) noexcept { x = internal::add(x, rhs.x); y = internal::add(y, rhs.y); return *this; } point &operator-=(const point rhs) noexcept { x = internal::add(x, -rhs.x); y = internal::add(y, -rhs.y); return *this; } point &operator*=(const point rhs) noexcept { long double _x = internal::add(x * rhs.x, -y * rhs.y); long double _y = internal::add(x * rhs.y, y * rhs.x); x = _x; y = _y; return *this; } point &operator*=(const long double k) noexcept { x *= k; y *= k; return *this; } point &operator/=(const long double k) { assert(internal::sgn(k) != 0); x /= k; y /= k; return *this; } point operator+(const point &rhs) const noexcept { return point(*this) += rhs; } point operator-(const point &rhs) const noexcept { return point(*this) -= rhs; } point operator*(const point &rhs) const noexcept { return point(*this) *= rhs; } point operator*(const long double rhs) const noexcept { return point(*this) *= rhs; } point operator/(const long double rhs) const { return point(*this) /= rhs; } point operator-() const noexcept { return point(0, 0) - *this; } long double abs() const noexcept { return std::sqrt(internal::add(x * x, y * y)); } long double dot(const point rhs) const noexcept { return internal::add(x * rhs.x, y * rhs.y); } long double det(const point rhs) const noexcept { return internal::add(x * rhs.y, -y * rhs.x); } // arctan(y/x) (単位はラジアン) long double arg() const { return std::atan2(y, x); } // x昇順, その後y昇順 bool operator<(const point &rhs) const noexcept { if (internal::sgn(x - rhs.x)) return internal::sgn(x - rhs.x) < 0; return internal::sgn(y - rhs.y) < 0; } bool operator==(const point &rhs) const noexcept { if (internal::sgn(x - rhs.x) == 0 && internal::sgn(y - rhs.y) == 0) return true; else return false; } bool operator!=(const point &rhs) const noexcept { return !((*this) == rhs); } }; std::ostream &operator<<(std::ostream &os, const point &a) { return os << a.x << " " << a.y; } std::istream &operator>>(std::istream &os, point &a) { return os >> a.x >> a.y; } point conj(const point &a) { return point(a.x, -a.y); } // 点a をang(ラジアン)回転する point rot(const point &a, long double ang) { return point(std::cos(ang) * a.x - std::sin(ang) * a.y, std::sin(ang) * a.x + std::cos(ang) * a.y); } point rot90(const point &a) { return point(-a.y, a.x); } long double dot(const point &a, const point &b) { return a.dot(b); } long double det(const point &a, const point &b) { return a.det(b); } long double abs(const point &a) { return a.abs(); } long double norm(const point &a) { return internal::add(a.x * a.x, a.y * a.y); } int isp(const point &a, const point &b, const point &c) { int flag = internal::sgn(det(b - a, c - a)); if (flag == 0) { if (internal::sgn(dot(b - a, c - a)) < 0) return -2; if (internal::sgn(dot(a - b, c - b)) < 0) return +2; } return flag; } // 分割統治で最近点対を求める O(N log N) long double closest_pair(std::vector<point> p) { std::sort(p.begin(), p.end()); int n = p.size(); auto f = [&](auto &&self, int l, int r) -> long double { if (r - l == 1) { return 1e9; } int mid = (l + r) / 2; long double x = p[mid].x; long double d = std::min(self(self, l, mid), self(self, mid, r)); std::vector<point> b; b.reserve(r - l); int j = mid; for (int i = l; i < mid; i++) { while (j < r && p[j].y <= p[i].y) { b.emplace_back(p[j++]); } b.emplace_back(p[i]); } while (j < r) { b.emplace_back(p[j++]); } for (int i = 0; i < r - l; i++) { p[l + i] = b[i]; } b.clear(); for (int i = l; i < r; i++) { if (std::abs(p[i].x - x) >= d) continue; for (int j = int(b.size()) - 1; j >= 0; j--) { if (p[i].y - b[j].y >= d) break; d = std::min(d, abs(p[i] - b[j])); } b.emplace_back(p[i]); } return d; }; return f(f, 0, n); } // ∠ABCを求める(ラジアン) long double angle(const point &A, const point &B, const point &C) { long double a = (B - C).abs(), b = (C - A).abs(), c = (A - B).abs(); long double cos = internal::add(internal::add(a * a, c * c), -b * b) / (2.0 * c * a); return std::acos(cos); } // 符号付き long double calc_ang(point a, point b, point c) { long double cos = dot((a - b), (c - b)) / (abs(a - b) * abs(c - b)); long double sin = det((a - b), (c - b)) / (abs(a - b) * abs(c - b)); return std::atan2(sin, cos); } void arg_sort(std::vector<point> &a) { int n = a.size(); std::vector ps(4, std::vector<point>()); auto idx = [](point v) -> int { if (v.y >= 0) return (v.x >= 0) ? 0 : 1; else return (v.x >= 0) ? 3 : 2; }; for (auto p : a) { assert(!(p.x == 0 && p.y == 0)); ps[idx(p)].emplace_back(p); } a.clear(); a.reserve(n); for (int i = 0; i < 4; i++) { std::sort(ps[i].begin(), ps[i].end(), [](point &p1, point &p2) -> bool { int flag = internal::sgn(internal::add(p1.x * p2.y, -p2.x * p1.y)); return flag == 0 ? (norm(p1) < norm(p2)) : flag > 0; }); for (auto &p : ps[i]) a.emplace_back(p); } return; } template <class T> void arg_sort_ll(std::vector<std::pair<T, T>> &a) { using Point = std::pair<T, T>; int n = a.size(); std::vector ps(4, std::vector<Point>()); auto idx = [](Point v) -> int { if (v.second >= 0) return (v.first >= 0) ? 0 : 1; else return (v.first >= 0) ? 3 : 2; }; for (auto p : a) { assert(!(p.first == 0 && p.second == 0)); ps[idx(p)].emplace_back(p); } a.clear(); a.reserve(n); for (int i = 0; i < 4; i++) { std::sort(ps[i].begin(), ps[i].end(), [](Point &p1, Point &p2) -> bool { T flag = p1.first * p2.second - p2.first * p1.second; return flag == 0 ? (p1.first * p1.first + p1.second * p1.second < p2.first * p2.first + p2.second * p2.second) : flag > 0; }); for (auto &p : ps[i]) a.emplace_back(p); } return; } } // namespace ebi #line 7 "geometry/line.hpp" namespace ebi { struct line { point a, b; line(long double x1, long double y1, long double x2, long double y2) : a(x1, y1), b(x2, y2) {} line(const point &a, const point &b) : a(a), b(b) {} point proj(const point &p) const { return a + (b - a) * (dot(b - a, p - a) / norm(b - a)); } point relf(const point &p) const { return proj(p) * double(2) - p; } long double distance(const point &c) const { return std::abs(det(c - a, b - a) / abs(b - a)); } }; int intersection(const line &a, const line &b) { if (internal::sgn(det(a.b - a.a, b.a - b.b)) != 0) { if (internal::sgn(dot(a.b - a.a, b.b - b.a)) == 0) { // 垂直 return 1; } return 0; // 交差 } else if (internal::sgn(det(a.b - a.a, b.a - a.a)) != 0) { // 平行 return 2; } else { // 同一直線 return 3; } } point cross_point(const point &a, const point &b, const point &c, const point &d) { return a + (b - a) * det(c - a, d - c) / det(b - a, d - c); } // 交点があるか確認する! point cross_point(const line &s, const line &t) { assert(intersection(s, t) < 2); return s.a + (s.b - s.a) * det(t.a - s.a, t.b - t.a) / det(s.b - s.a, t.b - t.a); } // 直線aと点cの距離 long double distance(const line &a, const point &c) { return std::abs(det(c - a.a, a.b - a.a) / abs(a.b - a.a)); } long double distance(const line &a, const line &b) { if (intersection(a, b) < 2) { return 0; } else { return distance(a, b.a); } } } // namespace ebi #line 2 "geometry/line_segment.hpp" #line 5 "geometry/line_segment.hpp" #line 8 "geometry/line_segment.hpp" namespace ebi { struct line_segment { point a, b; line_segment() = default; line_segment(long double x1, long double y1, long double x2, long double y2) : a(x1, y1), b(x2, y2) {} line_segment(const point &a, const point &b) : a(a), b(b) {} }; // 線分ab, cd が交わるか判定 bool intersection_line_segment(const point &a, const point &b, const point &c, const point &d) { if (internal::sgn(isp(a, b, c) * isp(a, b, d)) <= 0 && internal::sgn(isp(c, d, a) * isp(c, d, b)) <= 0) { return true; } return false; } // 線分ab, cd が交わるか判定 bool intersection(const line_segment &a, const line_segment &b) { return intersection_line_segment(a.a, a.b, b.a, b.b); } bool intersection(const line &a, const line_segment &b) { if (internal::sgn(det(a.b - a.a, b.a - a.a)) * internal::sgn(det(a.b - a.a, b.b - a.a)) < 0) { return true; } else { return false; } } point cross_point(const line_segment &s, const line_segment &t) { assert(intersection(s, t)); return s.a + (s.b - s.a) * det(t.a - s.a, t.b - t.a) / det(s.b - s.a, t.b - t.a); } long double distance(const line_segment &a, const point &c) { if (internal::sgn(dot(a.b - a.a, c - a.a)) < 0) { return abs(c - a.a); } else if (internal::sgn(dot(a.a - a.b, c - a.b)) < 0) { return abs(c - a.b); } else { return std::abs(det(c - a.a, a.b - a.a) / abs(a.b - a.a)); } } long double distance(const line_segment &a, const line_segment &b) { if (intersection(a, b)) { return 0; } else { return std::min(std::min(distance(a, b.a), distance(a, b.b)), std::min(distance(b, a.a), distance(b, a.b))); } } long double distance(const line &a, const line_segment &b) { if (intersection(a, b)) { return 0; } else { return std::min(distance(a, b.a), distance(a, b.b)); } } } // namespace ebi #line 2 "geometry/polygon.hpp" #line 5 "geometry/polygon.hpp" #line 8 "geometry/polygon.hpp" namespace ebi { using Polygon = std::vector<point>; long double area(const Polygon &poly) { long double s = 0; int n = poly.size(); for (int i = 0; i < n; i++) { s = internal::add(s, det(poly[i], poly[(i + 1 != n) ? i + 1 : 0])); } s /= 2.0; return s; } // 凸多角形か判定. pに反時計回りで点が入っていると仮定 bool is_convex(const Polygon &poly) { int n = poly.size(); for (int i = 0; i < n; i++) { if (isp(poly[i], poly[(i + 1 != n) ? i + 1 : 0], poly[(i + 2 < n) ? i + 2 : (i + 2) % n]) == -1) { return false; } } return true; } enum { OUT, ON, IN }; int contains(const Polygon &poly, const point &p) { bool in = false; int n = poly.size(); for (int i = 0; i < n; i++) { point a = poly[i] - p; point b = poly[(i + 1 < n) ? i + 1 : 0] - p; if (a.y > b.y) std::swap(a, b); if (a.y <= 0 && 0 < b.y && det(a, b) < 0) in = !in; if (det(a, b) == 0 && dot(a, b) <= 0) return ON; } return in ? IN : OUT; } long double convex_diameter(const Polygon &convex_poly) { int n = convex_poly.size(); int is = 0, js = 0; for (int i = 1; i < n; i++) { if (convex_poly[i].y > convex_poly[is].y) is = i; if (convex_poly[i].y < convex_poly[js].y) js = i; } long double max = (convex_poly[is] - convex_poly[js]).abs(); int i, max_i, j, max_j; i = max_i = is; j = max_j = js; do { if (det(convex_poly[(i + 1 < n) ? i + 1 : 0] - convex_poly[i], convex_poly[(j + 1 < n) ? j + 1 : 0] - convex_poly[j]) >= 0) { j = (j + 1 < n) ? j + 1 : 0; } else { i = (i + 1 < n) ? i + 1 : 0; } if ((convex_poly[i] - convex_poly[j]).abs() > max) { max = (convex_poly[i] - convex_poly[j]).abs(); max_i = i; max_j = j; } } while (i != is || j != js); return max; } Polygon convex_polygon_cut(const Polygon &poly, const line &l) { Polygon ret; int n = poly.size(); for (int i = 0; i < n; i++) { const point &now = poly[i]; const point &nxt = poly[(i + 1 < n) ? i + 1 : 0]; long double cf = det(l.a - now, l.b - now); long double cs = det(l.a - nxt, l.b - nxt); if (internal::sgn(cf) >= 0) { ret.emplace_back(now); } if (internal::sgn(cf) * internal::sgn(cs) < 0) { ret.emplace_back(cross_point(line(now, nxt), l)); } } return ret; } } // namespace ebi #line 10 "geometry/circle.hpp" namespace ebi { struct circle { point c; long double r; circle() = default; circle(const point &c, long double r) : c(c), r(r) {} }; int intersection(const circle &c1, const circle &c2) { long double d = abs(c1.c - c2.c); long double r1 = c1.r, r2 = c2.r; if (r1 < r2) std::swap(r1, r2); if (d > internal::add(r1, r2)) { return 4; } else if (d == internal::add(r1, r2)) { return 3; } else if (d > internal::add(r1, -r2)) { return 2; } else if (d == internal::add(r1, -r2)) { return 1; } else { return 0; } } circle incircle_of_triangle(const point &A, const point &B, const point &C) { long double a = abs(B - C), b = abs(C - A), c = abs(A - B); point in = A * a + B * b + C * c; in /= (a + b + c); long double r = distance(line(A, B), in); return circle(in, r); } circle circumscribed_circle_of_triangle(const point &A, const point &B, const point &C) { line p((A + B) / 2, (A + B) / 2 + rot90(B - A)); line q((B + C) / 2, (B + C) / 2 + rot90(C - B)); point cross = cross_point(p, q); return circle(cross, abs(A - cross)); } std::vector<point> cross_point(const circle &c, const line &l) { std::vector<point> ps; long double d = distance(l, c.c); if (d == c.r) { ps.emplace_back(l.proj(c.c)); } else if (d < c.r) { point p = l.proj(c.c); point v = l.b - l.a; v = v * std::sqrt( std::max(internal::add(c.r * c.r, -d * d), (long double)0)) / v.abs(); ps.emplace_back(p + v); ps.emplace_back(p - v); } return ps; } std::vector<point> cross_point(const circle &c1, const circle &c2) { std::vector<point> ps; int flag = intersection(c1, c2); if (flag == 0 || flag == 4) { return ps; } long double d = (c2.c - c1.c).abs(); long double x = internal::add(internal::add(d * d, c1.r * c1.r), -c2.r * c2.r) / (2.0 * d); point p = c1.c + (c2.c - c1.c) * x / d; point v = rot90(c2.c - c1.c); if (flag == 1 || flag == 3) { ps.emplace_back(p); } else { v = v * std::sqrt( std::max(internal::add(c1.r * c1.r, -x * x), (long double)0)) / v.abs(); ps.emplace_back(p + v); ps.emplace_back(p - v); } return ps; } std::vector<point> cross_point(const circle &c, const line_segment &ls) { std::vector<point> ps = cross_point(c, line(ls.a, ls.b)); std::vector<point> ret; for (auto p : ps) { if (internal::sgn(distance(ls, p)) == 0) { ret.emplace_back(p); } } return ret; } std::vector<point> tangent_to_circle(const circle &c, const point &p) { std::vector<point> ps; point v = p - c.c; long double d = v.abs(); long double h = std::sqrt( std::max(internal::add(norm(v), -c.r * c.r), (long double)(0.0))); long double cos = c.r / d, sin = h / d; point f(internal::add(v.x * cos, -v.y * sin), internal::add(v.x * sin, v.y * cos)); point g(internal::add(v.x * cos, v.y * sin), internal::add(-v.x * sin, v.y * cos)); f = f * c.r / f.abs(); g = g * c.r / g.abs(); ps.emplace_back(c.c + f); ps.emplace_back(c.c + g); return ps; } std::vector<line> tangent(circle c1, circle c2) { std::vector<line> ret; int flag = intersection(c1, c2); if (flag == 2 || flag == 3 || flag == 4) { if (c1.r == c2.r) { point v = c2.c - c1.c; v = rot90(v * c1.r / abs(v)); ret.emplace_back(line(c1.c + v, c2.c + v)); ret.emplace_back(line(c1.c - v, c2.c - v)); } else { point v = c1.c * (-c2.r) + c2.c * c1.r; v = v / (c1.r - c2.r); auto bs = tangent_to_circle(c1, v); auto as = tangent_to_circle(c2, v); for (int i = 0; i < (int)std::min(bs.size(), as.size()); i++) { ret.emplace_back(line(bs[i], as[i])); } } } else if (flag == 1) { point v = c2.c - c1.c; if (c1.r > c2.r) v = v * c1.r / v.abs(); else v = v * (-c1.r) / v.abs(); point p = c1.c + v; ret.emplace_back(line(p, p + rot90(v))); } if (flag == 4) { point p = c1.c * c2.r + c2.c * c1.r; p = p / (c1.r + c2.r); auto bs = tangent_to_circle(c1, p); auto as = tangent_to_circle(c2, p); for (int i = 0; i < (int)std::min(bs.size(), as.size()); i++) { ret.emplace_back(line(bs[i], as[i])); } } else if (flag == 3) { point v = c2.c - c1.c; v = v * c1.r / v.abs(); point p = c1.c + v; ret.emplace_back(p, p + rot90(v)); } return ret; } long double calc_common_area(const circle &c, const Polygon &poly) { if (int(poly.size()) < 3) return 0.0; long double s = 0; int n = poly.size(); auto cross_area = [&](auto &&self, const point &a, const point b) -> long double { point va = c.c - a, vb = c.c - b; long double f = det(va, vb) / 2.0, ret = 0.0; if (internal::sgn(f) == 0) return 0.0; if (std::max(va.abs(), vb.abs()) < c.r + EPS) return f; if (internal::sgn( internal::add(distance(line_segment(a, b), c.c), -c.r)) >= 0) { return c.r * c.r * (vb * conj(va)).arg() / 2.0; } auto ps = cross_point(c, line_segment(a, b)); if (int(ps.size()) == 1) ps.emplace_back(ps.back()); std::vector<point> tot = {a, ps[0], ps[1], b}; std::sort(tot.begin(), tot.end()); if (b < a) { std::reverse(tot.begin(), tot.end()); } int sz = tot.size(); for (int i = 0; i + 1 < sz; i++) { ret += self(self, tot[i], tot[i + 1]); } return ret; }; for (int i = 0; i < n; i++) { s += cross_area(cross_area, poly[i], poly[(i + 1 < n) ? i + 1 : 0]); } return s; } long double common_area(const circle &c1, const circle &c2) { int flag = intersection(c1, c2); if (flag == 3 || flag == 4) return 0.0; else if (flag == 0 || flag == 1) { long double r = std::min(c1.r, c2.r); return PI * r * r; } else { long double d = (c1.c - c2.c).abs(); long double theta1 = std::acos( internal::add(internal::add(c1.r * c1.r, d * d), -c2.r * c2.r) / (2.0 * c1.r * d)); long double area1 = internal::add( c1.r * c1.r * theta1, -c1.r * c1.r * std::sin(theta1 * 2) / 2.0); long double theta2 = std::acos( internal::add(internal::add(c2.r * c2.r, d * d), -c1.r * c1.r) / (2.0 * c2.r * d)); long double area2 = internal::add( c2.r * c2.r * theta2, -c2.r * c2.r * std::sin(theta2 * 2) / 2.0); return internal::add(area1, area2); } } } // namespace ebi #line 13 "test/geometry/tangent_to_circle.test.cpp" namespace ebi { using i64 = std::int64_t; void main_() { circle c; point p; std::cin >> p; std::cin >> c.c >> c.r; auto ps = tangent_to_circle(c, p); std::sort(ps.begin(), ps.end()); for (auto p : ps) { std::cout << p << '\n'; } } } // namespace ebi int main() { std::cout << std::fixed << std::setprecision(15); std::cin.tie(nullptr); std::ios::sync_with_stdio(false); ebi::main_(); }