/
geometry.hpp
630 lines (581 loc) · 25.6 KB
/
geometry.hpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
#ifndef SUISEN_GEOMETRY
#define SUISEN_GEOMETRY
#include <algorithm>
#include <cassert>
#include <complex>
#include <iostream>
#include <optional>
#include <tuple>
#include <variant>
#include <vector>
namespace suisen {
namespace geometry {
using coordinate_t = long double;
using Point = std::complex<coordinate_t>;
coordinate_t getx(const Point& p) { return p.real(); }
coordinate_t gety(const Point& p) { return p.imag(); }
// operator
Point operator+(const Point &p, coordinate_t real) { return Point(p) + Point(real, 0); }
Point operator-(const Point &p, coordinate_t real) { return Point(p) - Point(real, 0); }
Point operator*(const Point &p, coordinate_t real) { return Point(p) * Point(real, 0); }
Point operator/(const Point &p, coordinate_t real) { return Point(p) / Point(real, 0); }
Point operator+(coordinate_t real, const Point &p) { return Point(real, 0) + Point(p); }
Point operator-(coordinate_t real, const Point &p) { return Point(real, 0) - Point(p); }
Point operator*(coordinate_t real, const Point &p) { return Point(real, 0) * Point(p); }
Point operator/(coordinate_t real, const Point &p) { return Point(real, 0) / Point(p); }
std::istream& operator>>(std::istream &in, Point &p) {
coordinate_t x, y;
in >> x >> y;
p = Point(x, y);
return in;
}
std::ostream& operator<<(std::ostream &out, const Point &p) {
return out << getx(p) << ' ' << gety(p);
}
}
namespace geometry {
// relations between three points X, Y, Z.
struct ISP {
static constexpr int L_CURVE = +1; // +---------------+ Z is in 'a' => ISP = +1
static constexpr int R_CURVE = -1; // |aaaaaaaaaaaaaaa| Z is in 'b' => ISP = -1
static constexpr int FRONT = +2; // |ddd X eee Y ccc| Z is in 'c' => ISP = +2
static constexpr int BACK = -2; // |bbbbbbbbbbbbbbb| Z is in 'd' => ISP = -2
static constexpr int MIDDLE = 0; // +---------------+ Z is in 'e' => ISP = 0
};
struct Sign {
static constexpr int NEGATIVE = -1;
static constexpr int ZERO = 0;
static constexpr int POSITIVE = +1;
};
enum class Containment {
OUT, ON, IN
};
constexpr Point ZERO = Point(0, 0);
constexpr Point ONE = Point(1, 0);
constexpr Point I = Point(0, 1);
constexpr coordinate_t EPS = 1e-9;
constexpr coordinate_t PI = 3.14159265358979323846264338327950288419716939937510L;
constexpr coordinate_t E = 2.71828182845904523536028747135266249775724709369995L;
constexpr auto XY_COMPARATOR = [](const Point &p, const Point &q) {
return getx(p) == getx(q) ? gety(p) < gety(q) : getx(p) < getx(q);
};
constexpr auto XY_COMPARATOR_GREATER = [](const Point &p, const Point &q) {
return getx(p) == getx(q) ? gety(p) > gety(q) : getx(p) > getx(q);
};
constexpr auto YX_COMPARATOR = [](const Point &p, const Point &q) {
return gety(p) == gety(q) ? getx(p) < getx(q) : gety(p) < gety(q);
};
constexpr auto YX_COMPARATOR_GREATER = [](const Point &p, const Point &q) {
return gety(p) == gety(q) ? getx(p) > getx(q) : gety(p) > gety(q);
};
int sgn(coordinate_t x) {
return x > EPS ? Sign::POSITIVE : x < -EPS ? Sign::NEGATIVE : Sign::ZERO;
}
int compare(coordinate_t x, coordinate_t y) {
return sgn(x - y);
}
auto cartesian(const coordinate_t real, const coordinate_t imag) {
return Point(real, imag);
}
auto polar(const coordinate_t rho, const coordinate_t theta) {
return Point(rho * std::cos(theta), rho * std::sin(theta));
}
auto cis(const coordinate_t theta) {
return Point(std::cos(theta), std::sin(theta));
}
auto conj(const Point &z) {
return Point(getx(z), -gety(z));
}
auto arg(const Point &z) {
return std::atan2(gety(z), getx(z));
}
auto square_abs(const Point &z) {
return getx(z) * getx(z) + gety(z) * gety(z);
}
auto abs(const Point &z) {
return std::sqrt(square_abs(z));
}
auto rot(const Point &z, const coordinate_t theta) {
return cis(theta) * z;
}
auto dot(const Point &a, const Point &b) {
return getx(a) * getx(b) + gety(a) * gety(b);
}
auto det(const Point &a, const Point &b) {
return getx(a) * gety(b) - gety(a) * getx(b);
}
bool equals(const Point &a, const Point &b) {
return sgn(getx(a) - getx(b)) == Sign::ZERO and sgn(gety(a) - gety(b)) == Sign::ZERO;
}
bool equals(coordinate_t a, coordinate_t b) {
return compare(a, b) == 0;
}
// https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_1_C
int isp(const Point &a, const Point &b, const Point &c) {
Point ab = b - a, ac = c - a;
int s = sgn(det(ab, ac));
if (s == Sign::POSITIVE) return ISP::L_CURVE;
if (s == Sign::NEGATIVE) return ISP::R_CURVE;
if (sgn(dot(ab, ac)) == Sign::NEGATIVE) return ISP::BACK;
Point ba = a - b, bc = c - b;
if (sgn(dot(ba, bc)) == Sign::NEGATIVE) return ISP::FRONT;
return ISP::MIDDLE;
}
struct Line {
Point a, b;
Line() : Line(ZERO, ZERO) {}
Line(const Point &from, const Point &to) : a(from), b(to) {}
// coef_x * x + coef_y * y + cnst = 0
Line(coordinate_t coef_x, coordinate_t coef_y, coordinate_t cnst) {
if (not equals(coef_x, 0.)) {
a = { (coef_y - cnst) / coef_x, -1. };
b = { (-coef_y - cnst) / coef_x, +1. };
} else {
a = { -1., (coef_x - cnst) / coef_y };
b = { +1., (-coef_x - cnst) / coef_y };
}
}
template <size_t I> const std::tuple_element_t<I, Line>& get() const {
static_assert(I < std::tuple_size_v<Line>);
if constexpr (I == 0) return a;
else if constexpr (I == 1) return b;
}
template <size_t I> std::tuple_element_t<I, Line>& get() {
static_assert(I < std::tuple_size_v<Line>);
if constexpr (I == 0) return a;
else if constexpr (I == 1) return b;
}
};
struct Ray {
Point a, b;
Ray() : Ray(ZERO, ZERO) {}
Ray(const Point &from, const Point &to) : a(from), b(to) {}
template <size_t I> const std::tuple_element_t<I, Ray>& get() const {
static_assert(I < std::tuple_size_v<Ray>);
if constexpr (I == 0) return a;
else if constexpr (I == 1) return b;
}
template <size_t I> std::tuple_element_t<I, Ray>& get() {
static_assert(I < std::tuple_size_v<Ray>);
if constexpr (I == 0) return a;
else if constexpr (I == 1) return b;
}
};
struct Segment {
Point a, b;
Segment() : Segment(ZERO, ZERO) {}
Segment(const Point &from, const Point &to) : a(from), b(to) {}
template <size_t I> const std::tuple_element_t<I, Segment>& get() const {
static_assert(I < std::tuple_size_v<Segment>);
if constexpr (I == 0) return a;
else if constexpr (I == 1) return b;
}
template <size_t I> std::tuple_element_t<I, Segment>& get() {
static_assert(I < std::tuple_size_v<Segment>);
if constexpr (I == 0) return a;
else if constexpr (I == 1) return b;
}
};
struct Circle {
Point center;
coordinate_t radius;
Circle() : Circle(ZERO, 0) {}
Circle(const Point &c, const coordinate_t &r) : center(c), radius(r) {}
};
// Triangle
coordinate_t signed_area(const Point &a, const Point &b, const Point &c) {
return det(b - a, c - a) / 2;
}
coordinate_t area(const Point &a, const Point &b, const Point &c) {
return std::abs(signed_area(a, b, c));
}
Point pG(const Point &a, const Point &b, const Point &c) {
return (a + b + c) / 3;
}
// https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_B
Circle pI(const Point &a, const Point &b, const Point &c) {
coordinate_t la = abs(b - c), lb = abs(c - a), lc = abs(a - b);
coordinate_t l = la + lb + lc;
la /= l, lb /= l, lc /= l;
Point center = la * a + lb * b + lc * c;
coordinate_t radius = 2 * area(a, b, c) / l;
return Circle(center, radius);
}
// https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_C
Circle pO(const Point &a, const Point &b, const Point &c) {
Point ab = b - a, bc = c - b, ca = a - c;
coordinate_t la = square_abs(bc), lb = square_abs(ca), lc = square_abs(ab);
coordinate_t s = la * (lb + lc - la), t = lb * (lc + la - lb), u = lc * (la + lb - lc);
coordinate_t l = s + t + u;
s /= l, t /= l, u /= l;
Point center = a * s + b * t + c * u;
return Circle(center, abs(center - a));
}
Point pH(const Point &a, const Point &b, const Point &c) {
return a + b + c - 2 * pO(a, b, c).center;
}
auto pIabc(const Point &a, const Point &b, const Point &c) {
return std::make_tuple(pI(-a, b, c), pI(a, -b, c), pI(a, b, -c));
}
// Line
// https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_2_A
template <typename line_t_1, typename line_t_2>
auto is_parallel(const line_t_1 &l1, const line_t_2 &l2) -> decltype(l1.a, l1.b, l2.a, l2.b, bool()) {
return sgn(det(l1.b - l1.a, l2.b - l2.a)) == Sign::ZERO;
}
// https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_2_A
template <typename line_t_1, typename line_t_2>
auto is_orthogonal(const line_t_1 &l1, const line_t_2 &l2) -> decltype(l1.a, l1.b, l2.a, l2.b, bool()) {
return sgn(dot(l1.b - l1.a, l2.b - l2.a)) == Sign::ZERO;
}
template <typename line_t_1, typename line_t_2>
auto on_the_same_line(const line_t_1 &l1, const line_t_2 &l2) -> decltype(l1.a, l1.b, l2.a, l2.b, bool()) {
return is_parallel(l1, l2) and sgn(det(l1.b - l1.a, l2.a - l1.a)) == Sign::ZERO;
}
// https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_1_A
template <typename line_t>
Point projection(const Point &p, const line_t &line) {
Point a = p - line.a;
Point b = line.b - line.a;
return line.a + dot(a, b) / square_abs(b) * b;
}
// https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_1_B
template <typename line_t>
Point reflection(const Point &p, const line_t &line) {
Point h = projection(p, line);
return p + (h - p) * 2;
}
// https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_2_D
coordinate_t square_distance(const Point &p, const Segment &l) {
Point h = projection(p, l);
if (isp(l.a, l.b, h) == ISP::MIDDLE) {
return square_abs(h - p);
} else {
return std::min(square_abs(p - l.a), square_abs(p - l.b));
}
}
coordinate_t square_distance(const Segment &l, const Point &p) {
return square_distance(p, l);
}
coordinate_t square_distance(const Point &p, const Ray &l) {
Point h = projection(p, l);
int dir = isp(l.a, l.b, h);
return dir == ISP::MIDDLE or dir == ISP::FRONT ? square_abs(h - p) : std::min(square_abs(p - l.a), square_abs(p - l.b));
}
coordinate_t square_distance(const Ray &l, const Point &p) {
return square_distance(p, l);
}
coordinate_t square_distance(const Point &p, const Line &l) {
return square_abs(projection(p, l) - p);
}
coordinate_t square_distance(const Line &l, const Point &p) {
return square_distance(p, l);
}
// https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_2_D
coordinate_t distance(const Point &p, const Segment &l) {
return std::sqrt(square_distance(p, l));
}
// https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_2_D
coordinate_t distance(const Segment &l, const Point &p) {
return distance(p, l);
}
coordinate_t distance(const Point &p, const Ray &l) {
return std::sqrt(square_distance(p, l));
}
coordinate_t distance(const Ray &l, const Point &p) {
return distance(p, l);
}
coordinate_t distance(const Point &p, const Line &l) {
return std::sqrt(square_distance(p, l));
}
coordinate_t distance(const Line &l, const Point &p) {
return distance(p, l);
}
Containment contains(const Segment &l, const Point &p) {
return sgn(distance(l, p)) == 0 ? Containment::ON : Containment::OUT;
}
Containment contains(const Ray &l, const Point &p) {
return sgn(distance(l, p)) == 0 ? Containment::ON : Containment::OUT;
}
Containment contains(const Line &l, const Point &p) {
return sgn(distance(l, p)) == 0 ? Containment::ON : Containment::OUT;
}
bool equals(const Line &l, const Line &m) {
return on_the_same_line(l, m);
}
bool equals(const Ray &l, const Ray &m) {
return on_the_same_line(l, m) and equals(l.a, m.a);
}
bool equals(const Segment &l, const Segment &m) {
return (equals(l.a, m.a) and equals(l.b, m.b)) or (equals(l.a, m.b) and equals(l.b, m.a));
}
// "https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_2_B"
bool has_common_point(const Segment &l1, const Segment &l2) {
int isp_1a = isp(l1.a, l1.b, l2.a), isp_1b = isp(l1.a, l1.b, l2.b);
if (isp_1a * isp_1b > 0) return false;
int isp_2a = isp(l2.a, l2.b, l1.a), isp_2b = isp(l2.a, l2.b, l1.b);
if (isp_2a * isp_2b > 0) return false;
return true;
}
namespace internal {
template <typename line_t_1, typename line_t_2>
Point cross_point(const line_t_1 &l1, const line_t_2 &l2) {
assert(not is_parallel(l1, l2));
Point u = l1.b - l1.a, v = l2.a - l2.b, c = l2.a - l1.a;
return l2.a - det(u, c) / det(u, v) * v;
}
}
// "https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_2_C"
std::variant<std::nullptr_t, Point, Segment> common_point(const Segment &l1, const Segment &l2) {
if (not has_common_point(l1, l2)) return nullptr;
if (not is_parallel(l1, l2)) return internal::cross_point(l1, l2);
std::vector<Point> ps { l1.a, l1.b, l2.a, l2.b };
for (int i = 0; i <= 2; ++i) for (int j = 2; j >= i; --j) {
if (XY_COMPARATOR(ps[j + 1], ps[j])) std::swap(ps[j], ps[j + 1]);
}
if (equals(ps[1], ps[2])) return ps[1];
return Segment(ps[1], ps[2]);
}
// https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_2_D
coordinate_t square_distance(const Segment &l1, const Segment &l2) {
if (has_common_point(l1, l2)) return 0;
return std::min({ square_distance(l1, l2.a), square_distance(l1, l2.b), square_distance(l1.a, l2), square_distance(l1.b, l2) });
}
// https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_2_D
coordinate_t distance(const Segment &l1, const Segment &l2) {
return std::sqrt(square_distance(l1, l2));
}
// Polygon
using Polygon = std::vector<Point>;
// https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_3_A
coordinate_t signed_area(const Polygon &poly) {
coordinate_t res = 0;
int sz = poly.size();
for (int i = 0; i < sz; ++i) {
int j = i + 1;
if (j == sz) j = 0;
res += signed_area(ZERO, poly[i], poly[j]);
}
return res;
}
// https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_3_A
auto area(const Polygon &poly) {
return std::abs(signed_area(poly));
}
// https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_3_B
template <bool accept_180_degree = true>
bool is_convex(const Polygon &poly) {
int sz = poly.size();
for (int i = 0; i < sz; ++i) {
int j = i + 1, k = i + 2;
if (j >= sz) j -= sz;
if (k >= sz) k -= sz;
int dir = isp(poly[i], poly[j], poly[k]);
if constexpr (accept_180_degree) {
if (dir == ISP::R_CURVE) return false;
} else {
if (dir != ISP::L_CURVE) return false;
}
}
return true;
}
// https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_3_C
Containment contains(const Polygon &poly, const Point &p) {
bool in = false;
int sz = poly.size();
for (int i = 0; i < sz; ++i) {
int j = i + 1;
if (j == sz) j -= sz;
Point a = poly[i] - p, b = poly[j] - p;
if (gety(a) > gety(b)) std::swap(a, b);
if (sgn(gety(a)) <= 0 and sgn(gety(b)) > 0 and sgn(det(a, b)) < 0) in = not in;
if (sgn(det(a, b)) == 0 and sgn(dot(a, b)) <= 0) return Containment::ON;
}
return in ? Containment::IN : Containment::OUT;
}
// https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_4_B
std::tuple<int, int, coordinate_t> convex_diameter(const Polygon& convex) {
const int sz = convex.size();
if (sz <= 2) return { 0, sz - 1, abs(convex.front() - convex.back()) };
auto [si, sj] = [&]{
auto [it_min, it_max] = std::minmax_element(convex.begin(), convex.end(), XY_COMPARATOR);
return std::pair<int, int> { it_min - convex.begin(), it_max - convex.begin() };
}();
coordinate_t max_dist = -1;
std::pair<int, int> argmax{ -1, -1 };
for (int i = si, j = sj; i != sj or j != si;) {
if (coordinate_t dij = square_abs(convex[j] - convex[i]); dij > max_dist) max_dist = dij, argmax = { i, j };
int ni = (i + 1) % sz, nj = (j + 1) % sz;
if (det(convex[ni] - convex[i], convex[nj] - convex[j]) < 0) i = ni;
else j = nj;
}
return { argmax.first, argmax.second, ::sqrtl(max_dist) };
}
// https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_4_C
auto convex_cut(const Polygon &convex, const Line &l) {
Polygon res;
int sz = convex.size();
for (int i = 0; i < sz; ++i) {
int j = i + 1;
if (j == sz) j -= sz;
const Point &a = convex[i], &b = convex[j];
int da = sgn(det(l.b - l.a, a - l.a));
if (da >= 0) res.push_back(a);
int db = sgn(det(l.b - l.a, b - l.a));
if (da * db < 0) res.push_back(internal::cross_point(l, Segment(a, b)));
}
return res;
}
// Circle
// https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_A
int tangent_num(const Circle &c1, const Circle &c2) {
coordinate_t r1 = c1.radius, r2 = c2.radius;
if (r1 > r2) return tangent_num(c2, c1);
coordinate_t d = abs(c1.center - c2.center);
int cp = compare(d, r1 + r2);
if (cp > 0) return 4;
if (cp == 0) return 3;
int cn = compare(d, r2 - r1);
if (cn > 0) return 2;
if (cn == 0) return 1;
return 0;
}
// https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_D
std::vector<Point> common_point(const Circle &c, const Line &l) {
Point h = projection(c.center, l);
coordinate_t d = abs(c.center - h);
int cp = compare(d, c.radius);
if (cp > 0) return {};
if (cp == 0) return { h };
auto v = (l.a - l.b) * (std::sqrt(c.radius * c.radius - d * d) / abs(l.a - l.b));
return { h - v, h + v };
}
// https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_H
std::vector<Point> common_point(const Circle &c, const Segment &l) {
auto ps = common_point(c, Line(l.a, l.b));
ps.erase(std::remove_if(ps.begin(), ps.end(), [&](const auto &p) { return contains(l, p) != Containment::ON; }), ps.end());
return ps;
}
// https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_E
std::vector<Point> common_point(const Circle &c1, const Circle &c2) {
coordinate_t r1 = c1.radius, r2 = c2.radius;
if (r1 > r2) return common_point(c2, c1);
coordinate_t d = abs(c1.center - c2.center);
int cp = compare(d, r1 + r2), cn = compare(d, r2 - r1);
if (cp > 0 or cn < 0) return {};
auto v = c1.center - c2.center;
coordinate_t lv = abs(v);
if (cp == 0 or cn == 0) {
return { c2.center + v * (r2 / lv) };
}
coordinate_t lp = d, ln = (r2 * r2 - r1 * r1) / d;
coordinate_t p = (lp + ln) / 2, x = sqrt(r2 * r2 - p * p);
auto h = c2.center + v * (p / lv);
auto t = v * I;
return { h + t * (x / lv), h - t * (x / lv) };
}
// https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_F
Containment contains(const Circle &c, const Point &p) {
coordinate_t d = abs(c.center - p);
int cp = compare(d, c.radius);
if (cp > 0) return Containment::OUT;
if (cp < 0) return Containment::IN;
return Containment::ON;
}
// https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_F
std::vector<Point> tangent_to_circle(const Circle &c, const Point &p) {
Containment cnt = contains(c, p);
if (cnt == Containment::IN) return {};
if (cnt == Containment::ON) return { p };
Point v = c.center - p;
coordinate_t r = c.radius, d = abs(v), l = sqrt(d * d - r * r);
coordinate_t t = std::asin(r / d);
return { p + rot(v, t) * (l / d), p + rot(v, -t) * (l / d) };
}
// https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_G
// returns { Line(p, q) | p is on c1, q is on c2, Line(p, q) is common tangent of c1 and c2 }
std::vector<Line> common_tangent(const Circle &c1, const Circle &c2) {
int num = tangent_num(c1, c2);
std::vector<Line> res;
if (num == 0) return res;
Point a = c1.center, b = c2.center, v = b - a;
coordinate_t r1 = c1.radius, r2 = c2.radius;
coordinate_t rp = r1 + r2, rm = r1 - r2, rd = r2 / r1;
coordinate_t sqxy = square_abs(v);
coordinate_t rtp = std::sqrt(std::max(sqxy - rp * rp, coordinate_t(0)));
coordinate_t rtm = std::sqrt(std::max(sqxy - rm * rm, coordinate_t(0)));
Point r = v * r1, u = r * Point(0, 1);
Point l12 = r * rp, r12 = u * rtp, l34 = r * rm, r34 = u * rtm;
Point p14 = (l34 + r34) / sqxy;
res.emplace_back(a + p14, b + p14 * rd);
if (num == 1) return res;
Point p13 = (l34 - r34) / sqxy;
res.emplace_back(a + p13, b + p13 * rd);
if (num == 2) return res;
Point p12 = (l12 + r12) / sqxy;
res.emplace_back(a + p12, b - p12 * rd);
if (num == 3) return res;
Point p11 = (l12 - r12) / sqxy;
res.emplace_back(a + p11, b - p11 * rd);
return res;
}
// https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_H
coordinate_t intersection_area(const Polygon &poly, const Circle &circle) {
int sz = poly.size();
coordinate_t r2 = circle.radius * circle.radius;
const Point &c = circle.center;
coordinate_t area = 0;
for (int i = 0; i < sz; i++) {
int j = i + 1;
if (j >= sz) j -= sz;
Point a = poly[i], b = poly[j];
bool in_a = contains(circle, a) == Containment::IN, in_b = contains(circle, b) == Containment::IN;
Point ca = a - c, cb = b - c;
if (in_a and in_b) {
area += det(ca, cb);
continue;
}
std::vector<Point> ps = common_point(circle, Segment(a, b));
if (ps.empty()) {
area += r2 * arg(cb / ca);
} else {
Point s = ps[0];
Point t = ps.size() == 1 ? s : ps[1];
if (compare(square_abs(t - a), square_abs(s - a)) < 0) std::swap(s, t);
Point cs = s - c, ct = t - c;
area += det(cs, ct);
area += in_a ? det(ca, cs) : r2 * arg(cs / ca);
area += in_b ? det(ct, cb) : r2 * arg(cb / ct);
}
}
return area / 2;
}
// https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_I
coordinate_t intersection_area(const Circle &c1, const Circle &c2) {
coordinate_t r = c1.radius, s = c2.radius;
if (r < s) return intersection_area(c2, c1);
Point a = c1.center, b = c2.center;
coordinate_t d = abs(a - b);
if (compare(d, r + s) >= 0) return 0;
if (compare(d, r - s) <= 0) return PI * s * s;
coordinate_t x = (d * d + r * r - s * s) / (2 * d);
coordinate_t h = std::sqrt(std::max(r * r - x * x, coordinate_t(0)));
coordinate_t a1 = r * r * std::acos(x / r);
coordinate_t a2 = s * s * std::acos((d - x) / s);
coordinate_t a12 = d * h;
return a1 + a2 - a12;
}
}
} // namespace suisen
namespace std {
template <> struct tuple_size<suisen::geometry::Segment> { static constexpr size_t value = 2; };
template <> struct tuple_element<0, suisen::geometry::Segment> { using type = suisen::geometry::Point; };
template <> struct tuple_element<1, suisen::geometry::Segment> { using type = suisen::geometry::Point; };
template <> struct tuple_size<suisen::geometry::Ray> { static constexpr size_t value = 2; };
template <> struct tuple_element<0, suisen::geometry::Ray> { using type = suisen::geometry::Point; };
template <> struct tuple_element<1, suisen::geometry::Ray> { using type = suisen::geometry::Point; };
template <> struct tuple_size<suisen::geometry::Line> { static constexpr size_t value = 2; };
template <> struct tuple_element<0, suisen::geometry::Line> { using type = suisen::geometry::Point; };
template <> struct tuple_element<1, suisen::geometry::Line> { using type = suisen::geometry::Point; };
}
#endif // SUISEN_GEOMETRY