1import math
2from typing import Optional, Union, Final
3from decimal import Decimal, getcontext
4from titan_pylib.algorithm.sort.merge_sort import merge_sort
5from titan_pylib.math.decimal_util import (
6 decimal_pi,
7 decimal_sin,
8 decimal_cos,
9 decimal_acos,
10)
11
12
[docs]
13class GeometryUtil:
14
15 getcontext().prec = 10
16 USE_DECIMAL: Final[bool] = False
17
18 EPS: Final[Union[Decimal, float]] = (
19 Decimal("1e-" + str(getcontext().prec // 2)) if USE_DECIMAL else 1e-10
20 )
21 GEO_INF: Final[Union[Decimal, float]] = (
22 Decimal("INF") if USE_DECIMAL else float("inf")
23 )
24 sin = decimal_sin if USE_DECIMAL else math.sin
25 cos = decimal_cos if USE_DECIMAL else math.cos
26 acos = decimal_acos if USE_DECIMAL else math.acos
27 pi = decimal_pi() if USE_DECIMAL else math.pi
28
[docs]
29 @classmethod
30 def sqrt(cls, v: Union[Decimal, float]) -> None:
31 return math.sqrt(v) if not cls.USE_DECIMAL else v.sqrt()
32
[docs]
33 @classmethod
34 def eq(cls, a: float, b: float) -> bool:
35 return abs(a - b) < cls.EPS
36
[docs]
37 @classmethod
38 def lt(cls, u: float, v: float) -> bool:
39 return u < v and u < max(v, v - cls.EPS)
40
[docs]
41 @classmethod
42 def gt(cls, u: float, v: float) -> bool:
43 return u > v and u > min(v, v - cls.EPS)
44
[docs]
45 @classmethod
46 def le(cls, u: float, v: float) -> bool:
47 return cls.eq(u, v) or cls.lt(u, v)
48
[docs]
49 @classmethod
50 def ge(cls, u: float, v: float) -> bool:
51 return cls.eq(u, v) or cls.gt(u, v)
52
53
[docs]
54class Point:
55
56 __slots__ = "x", "y"
57
58 def __init__(self, x: float, y: float) -> None:
59 self.x = x
60 self.y = y
61
62 def __str__(self) -> str:
63 return f"Point({self.x}, {self.y})"
64
65 __repr__ = __str__
66
[docs]
67 def show(self, precision: int = 16) -> None:
68 format_string = f"{{:.{precision}f}}"
69 print(f"{format_string.format(self.x)} {format_string.format(self.y)}")
70
71 def __add__(self, other: "Point") -> "Point":
72 if isinstance(other, Point):
73 return Point(self.x + other.x, self.y + other.y)
74 else:
75 return Point(self.x + other, self.y + other)
76
77 def __sub__(self, other: Union["Point", float]) -> "Point":
78 if isinstance(other, Point):
79 return Point(self.x - other.x, self.y - other.y)
80 else:
81 return Point(self.x - other, self.y - other)
82
83 def __mul__(self, other: Union["Point", float]) -> "Point":
84 if isinstance(other, Point):
85 return Point(self.x * other.x, self.y * other.y)
86 else:
87 return Point(self.x * other, self.y * other)
88
89 __rmul__ = __mul__
90
91 def __truediv__(self, other: float) -> "Point":
92 return Point(self.x / other, self.y / other)
93
94 def __le__(self, other: "Point") -> bool:
95 if self == other:
96 return True
97 if GeometryUtil.lt(self.x, other.x):
98 return True
99 if GeometryUtil.lt(other.x, self.x):
100 return False
101 return GeometryUtil.lt(self.y, other.y)
102
103 def __lt__(self, other: "Point") -> bool:
104 if self == other:
105 return False
106 if GeometryUtil.lt(self.x, other.x):
107 return True
108 if GeometryUtil.lt(other.x, self.x):
109 return False
110 return GeometryUtil.lt(self.y, other.y)
111
112 def __ge__(self, other: "Point") -> bool:
113 return not (self < other)
114
115 def __gt__(self, other: "Point") -> bool:
116 return not (self <= other)
117
118 def __eq__(self, other: "Point") -> bool:
119 return (
120 abs(self.x - other.x) < GeometryUtil.EPS
121 and abs(self.y - other.y) < GeometryUtil.EPS
122 )
123
124 def __neg__(self) -> "Point":
125 return Point(-self.x, -self.y)
126
127 def __abs__(self) -> float:
128 return GeometryUtil.sqrt(self.x * self.x + self.y * self.y)
129
[docs]
130 def norm(self) -> float:
131 """norm"""
132 return abs(self)
133
[docs]
134 def norm2(self) -> float:
135 """(x**2 + y**2)を返す(ルートの計算はしない)"""
136 return self.x * self.x + self.y * self.y
137
[docs]
138 def unit(self) -> "Point":
139 return self / abs(self)
140
[docs]
141 def rotate(self, theta: float) -> "Point":
142 return Point(
143 GeometryUtil.cos(theta) * self.x - GeometryUtil.sin(theta) * self.y,
144 GeometryUtil.sin(theta) * self.x + GeometryUtil.cos(theta) * self.y,
145 )
146
[docs]
147 def copy(self) -> "Point":
148 return Point(self.x, self.y)
149
150
[docs]
151class Line:
152 """ax + by + c = 0"""
153
154 __slots__ = "a", "b", "c", "p1", "p2"
155
156 def __init__(self, a: int, b: int, c: int, p1: Point, p2: Point) -> None:
157 self.a = a
158 self.b = b
159 self.c = c
160 self.p1 = p1
161 self.p2 = p2
162
[docs]
163 @classmethod
164 def from_points(cls, p1: Point, p2: Point) -> "Line":
165 a = p2.y - p1.y
166 b = p1.x - p2.x
167 c = -p1.y * b - p1.x * a
168 if isinstance(a, int) and isinstance(b, int) and isinstance(c, int):
169 g = math.gcd(a, b)
170 g = math.gcd(g, c)
171 a //= g
172 b //= g
173 c //= g
174 a, b, c = max((a, b, c), (-a, -b, -c))
175 return cls(a, b, c, p1, p2)
176
[docs]
177 @classmethod
178 def from_abc(cls, a: int, b: int, c: int) -> "Line":
179 raise NotImplementedError
180
[docs]
181 def get_y(self, x: float) -> Optional[float]:
182 return None if GeometryUtil.eq(self.b, 0) else -(self.c + self.a * x) / self.b
183
[docs]
184 def get_x(self, y: float) -> Optional[float]:
185 return None if GeometryUtil.eq(self.a, 0) else -(self.c + self.b * y) / self.a
186
187
[docs]
188class Polygon:
189 """多角形クラス"""
190
191 def __init__(self, ps: list[Point]) -> None:
192 self.n = len(ps)
193 self.ps = ps
194
[docs]
195 def area(self) -> float:
196 """面積を求める / :math:`O(n)`"""
197 res = 0
198 p = self.ps
199 for i in range(self.n - 1):
200 res += Geometry.cross(p[i], p[i + 1])
201 if p:
202 res += Geometry.cross(p[-1], p[0])
203 return res / 2
204
[docs]
205 def is_convex(self) -> bool:
206 """凸多角形かどうか / :math:`O(n)`"""
207 ps = self.ps
208 for i in range(self.n):
209 pre = (i - 1 + self.n) % self.n
210 now = i
211 nxt = (i + 1) % self.n
212 if Geometry.ccw(ps[pre], ps[now], ps[nxt]) == -1:
213 return False
214 return True
215
[docs]
216 def get_degree(self, radian):
217 return radian * (180 / GeometryUtil.pi)
218
[docs]
219 def contains(self, p: Point) -> int:
220 """点の包含関係を返す / O(n)
221
222 Returns:
223 int: `2`: `p` を含む
224 `1`: `p` が辺上にある
225 `0`: それ以外
226 """
227 ps = self.ps
228 deg = 0
229 for i in range(self.n):
230 if ps[i] == ps[(i + 1) % self.n]:
231 continue
232 if GeometryUtil.eq(
233 Geometry.dist_p_to_segmemt(p, Segment(ps[i], ps[(i + 1) % self.n])), 0
234 ):
235 return 1
236 ax = ps[i].x - p.x
237 ay = ps[i].y - p.y
238 bx = ps[(i + 1) % self.n].x - p.x
239 by = ps[(i + 1) % self.n].y - p.y
240 cos = (ax * bx + ay * by) / (
241 GeometryUtil.sqrt((ax * ax + ay * ay) * (bx * bx + by * by))
242 )
243 cos = min(1, max(-1, cos))
244 deg += self.get_degree(GeometryUtil.acos(cos))
245 return 2 if GeometryUtil.eq(deg, 360) else 0
246
247
[docs]
248class ConvexPolygon(Polygon):
249
250 def __init__(self, ps: list[Point], line_contains: bool = False) -> None:
251 self.ps = Geometry.convex_hull(ps, line_contains)
252 self.n = len(self.ps)
253
[docs]
254 def diameter(self) -> tuple[float, int, int]:
255 ps = self.ps
256 n = len(ps)
257 if n == 0:
258 return -1, -1, -1
259 if n == 1:
260 return 0, ps[0], ps[0]
261 if n == 2:
262 return abs(ps[0] - ps[1]), ps[0], ps[1]
263 u, v = 0, 0
264 up, vp = None, None
265 for i in range(n):
266 if ps[u] > ps[i]:
267 u = i
268 up = ps[i]
269 if ps[v] < ps[i]:
270 v = i
271 vp = ps[i]
272 d = 0
273 su, sv = u, v
274 loop = False
275 while (u != su or v != sv) or (not loop):
276 loop = True
277 if GeometryUtil.lt(d, abs(ps[u] - ps[v])):
278 d = abs(ps[u] - ps[v])
279 up = ps[u]
280 vp = ps[v]
281 if GeometryUtil.lt(
282 Geometry.cross(ps[(u + 1) % n] - ps[u], ps[(v + 1) % n] - ps[v]), 0
283 ):
284 u = (u + 1) % n
285 else:
286 v = (v + 1) % n
287 return d, up, vp
288
[docs]
289 def contains(self, p: Point) -> int:
290 """点の包含関係を返す / :math:`O(\\log{n})`
291
292 Returns:
293 int: `2`: `p` を含む
294 `1`: `p` が辺上にある
295 `0`: それ以外
296 """
297 ps = self.ps
298 n = len(ps)
299 b1 = Geometry.cross(ps[1] - ps[0], p - ps[0])
300 b2 = Geometry.cross(ps[-1] - ps[0], p - ps[0])
301 if GeometryUtil.lt(b1, 0) or GeometryUtil.gt(b2, 0):
302 return 0
303 l, r = 1, n - 1
304 while r - l > 1:
305 mid = (l + r) >> 1
306 c = Geometry.cross(p - ps[0], ps[mid] - ps[0])
307 if GeometryUtil.eq(c, 0) or GeometryUtil.gt(c, 0):
308 r = mid
309 else:
310 l = mid
311 b3 = Geometry.cross(ps[l] - p, ps[r] - p)
312 if GeometryUtil.eq(b3, 0):
313 return 1
314 if GeometryUtil.gt(b3, 0):
315 if GeometryUtil.eq(b1, 0) or GeometryUtil.eq(b2, 0):
316 return 1
317 else:
318 return 2
319 return 0
320
[docs]
321 def convex_cut(self, l: Line) -> "ConvexPolygon":
322 """直線 ``l`` で切断したときの左側の凸多角形を返す"""
323 ret = []
324 for i in range(self.n):
325 now = self.ps[i]
326 nxt = self.ps[(i + 1) % self.n]
327 if Geometry.ccw(l.p1, l.p2, now) != -1:
328 ret.append(now.copy())
329 if Geometry.ccw(l.p1, l.p2, now) * Geometry.ccw(l.p1, l.p2, nxt) < 0:
330 p = Geometry.cross_point_line(Line.from_points(now, nxt), l)
331 assert p is not None
332 ret.append(p.copy())
333 return ConvexPolygon(ret)
334
335
[docs]
336class Circle:
337
338 __slots__ = "p", "r"
339
340 def __init__(self, p: Point, r: float) -> None:
341 """円
342
343 Args:
344 p (Point): 直径を表す `Point`
345 r (float): 半径
346 """
347 self.p = p
348 self.r = r
349
[docs]
350 def area(self) -> float:
351 """面積を返す"""
352 return GeometryUtil.pi * self.r * self.r
353
354
[docs]
355class Segment:
356
357 def __init__(self, p1: Point, p2: Point) -> None:
358 assert p1 != p2, f"{p1=}, {p2=}"
359 if p1 > p2:
360 p1, p2 = p2, p1
361 self.p1 = p1
362 self.p2 = p2
363
364
[docs]
365class Geometry:
366 """ref:
367 - https://hcpc-hokudai.github.io/arpsive/geometry_004.pdf
368 - https://bakamono1357.hatenablog.com/entry/2020/04/29/025320
369 - https://tjkendev.github.io/procon-library/python/geometry/circles_associated_with_triangle.html
370 - https://atcoder.jp/contests/abc296/editorial/6109
371 """
372
[docs]
373 @classmethod
374 def dot(cls, p1: Point, p2: Point) -> float:
375 return p1.x * p2.x + p1.y * p2.y
376
[docs]
377 @classmethod
378 def cross(cls, p1: Point, p2: Point) -> float:
379 return p1.x * p2.y - p1.y * p2.x
380
[docs]
381 @classmethod
382 def ccw(cls, u: Point, v: Point, p: Point) -> int:
383 """u->vに対し、v->pの位置関係を求める
384
385 Returns:
386 int: `+1`: a->bに対し、b->cが半時計回りに進む
387 `-1`: a->bに対し、b->cが時計回りに進む
388 `+2`: a->bに対し、b->cが逆方向に進む
389 `-2`: a->bに対し、b->cが同一方向に進む
390 ` 0`: cが線分a->b上にある
391 """
392 if cls.cross(v - u, p - u) > GeometryUtil.EPS:
393 return 1
394 if cls.cross(v - u, p - u) < -GeometryUtil.EPS:
395 return -1
396 if cls.dot(v - u, p - u) < -GeometryUtil.EPS:
397 return 2
398 if GeometryUtil.lt(abs(v - u), abs(p - u)):
399 # if abs(v - u) + GeometryUtil.EPS < abs(p - u):
400 return -2
401 return 0
402
[docs]
403 @classmethod
404 def projection_point(cls, p: Point, l: Line) -> Point:
405 """直線 `l` に、点 `p` からおろした垂線の足の `Point` を返す"""
406 t = cls.dot(p - l.p1, l.p1 - l.p2) / (l.p1 - l.p2).norm2()
407 return l.p1 + (l.p1 - l.p2) * t
408
[docs]
409 @classmethod
410 def reflection_point(cls, p: Point, l: Line) -> Point:
411 """直線 `l` を対象軸として点 `p` と線対称の点を返す"""
412 return p + 2 * (cls.projection_point(p, l) - p)
413
[docs]
414 @classmethod
415 def is_orthogonal(cls, l1: Line, l2: Line) -> bool:
416 """直線 `l1, l2` が直行しているかどうか"""
417 return GeometryUtil.eq(cls.dot(l1.p2 - l1.p1, l2.p2 - l2.p1), 0)
418
[docs]
419 @classmethod
420 def is_parallel(cls, l1: Line, l2: Line) -> bool:
421 """直線 `l1, l2` が平行かどうか"""
422 return GeometryUtil.eq(cls.cross(l1.p2 - l1.p1, l2.p2 - l2.p1), 0)
423
[docs]
424 @classmethod
425 def is_intersect_linesegment(cls, s1: Segment, s2: Segment) -> bool:
426 """線分 `s1` と `s2` が交差しているかどうか判定する"""
427 return (
428 cls.ccw(s1.p1, s1.p2, s2.p1) * cls.ccw(s1.p1, s1.p2, s2.p2) <= 0
429 and cls.ccw(s2.p1, s2.p2, s1.p1) * cls.ccw(s2.p1, s2.p2, s1.p2) <= 0
430 )
431
[docs]
432 @classmethod
433 def is_intersect_circle(cls, c1: Circle, c2: Circle) -> int:
434 """2円の位置関係
435
436 Returns:
437 int: 共通接線の数
438 `0`: 内包している
439 `1`: 内接している
440 `2`: 2点で交わっている
441 `3`: 外接している
442 `4`: 離れている
443 """
444 d = abs(c1.p - c2.p)
445 if GeometryUtil.gt(d, c1.r + c2.r):
446 # if d > c1.r + c2.r + GeometryUtil.EPS:
447 return 4
448 if GeometryUtil.eq(d, c1.r + c2.r):
449 return 3
450 if GeometryUtil.eq(d, abs(c1.r - c2.r)):
451 return 1
452 if GeometryUtil.lt(d, abs(c1.r - c2.r)):
453 # if d < abs(c1.r - c2.r) - GeometryUtil.EPS:
454 return 0
455 return 2
456
[docs]
457 @classmethod
458 def cross_point_line(cls, l1: Line, l2: Line) -> Optional[Point]:
459 assert isinstance(l1, Line)
460 assert isinstance(l2, Line)
461 d1 = cls.cross(l1.p2 - l1.p1, l2.p2 - l2.p1)
462 d2 = cls.cross(l1.p2 - l1.p1, l1.p2 - l2.p1)
463 if GeometryUtil.eq(abs(d1), 0) and GeometryUtil.eq(abs(d2), 0):
464 return l2.p1
465 if GeometryUtil.eq(abs(d1), 0):
466 return None
467 return l2.p1 + (l2.p2 - l2.p1) * (d2 / d1)
468
[docs]
469 @classmethod
470 def cross_point_line_segment(cls, l: Line, s: Segment) -> Optional[Point]:
471 raise NotImplementedError
472 assert isinstance(l, Line)
473 assert isinstance(s, Segment)
474 assert s.p1.x <= s.p2.x
475 if (
476 (not GeometryUtil.eq(s.p1.x, s.p2.x))
477 and l.get_y(s.p1.x) is not None
478 and l.get_y(s.p2.x) is not None
479 ):
480 sl = Segment(
481 Point(s.p1.x, l.get_y(s.p1.x)),
482 Point(s.p2.x, l.get_y(s.p2.x)),
483 )
484 return cls.cross_point_segment(sl, s)
485 p1, p2 = s.p1, s.p2
486 if p1.y > p2.y:
487 p1, p2 = p2, p1
488 if (
489 (not GeometryUtil.eq(p1.y, p2.y))
490 and l.get_x(p1.y) is not None
491 and l.get_x(p2.y) is not None
492 ):
493 sl = Segment(
494 Point(l.get_x(p1.y), p1.y),
495 Point(l.get_x(p2.y), p2.y),
496 )
497 return cls.cross_point_segment(sl, s)
498 return False
499
[docs]
500 @classmethod
501 def cross_point_segment(cls, s1: Segment, s2: Segment) -> Optional[Point]:
502 assert isinstance(s1, Segment)
503 assert isinstance(s2, Segment)
504 d1 = cls.cross(s1.p2 - s1.p1, s2.p2 - s2.p1)
505 d2 = cls.cross(s1.p2 - s1.p1, s1.p2 - s2.p1)
506 if GeometryUtil.eq(abs(d1), 0) and GeometryUtil.eq(abs(d2), 0):
507 return s2.p1
508 if GeometryUtil.eq(abs(d1), 0):
509 return None
510 res = s2.p1 + (s2.p2 - s2.p1) * (d2 / d1)
511 mn_x = max(s1.p1.x, s2.p1.x)
512 mx_x = min(s1.p2.x, s2.p2.x)
513 if not (GeometryUtil.le(mn_x, res.x) and GeometryUtil.le(res.x, mx_x)):
514 return None
515 return None
516 s1_p1, s1_p2 = s1.p1, s1.p2
517 if s1_p1.y > s1_p2.y:
518 s1_p1, s1_p2 = s1_p2, s1_p1
519 s2_p1, s2_p2 = s2.p1, s2.p2
520 if s2_p1.y > s2_p2.y:
521 s2_p1, s2_p2 = s2_p2, s2_p1
522 mn_y = max(s1_p1.y, s2_p1.y)
523 mx_y = min(s1_p2.y, s2_p2.y)
524 if not (GeometryUtil.le(mn_y, res.y) and GeometryUtil.le(res.y, mx_y)):
525 return None
526 return res
527
[docs]
528 @classmethod
529 def cross_point_circle(cls, c1: Circle, c2: Circle) -> tuple[Point, ...]:
530 r = cls.is_intersect_circle(c1, c2)
531 if r == 4 or r == 0:
532 return ()
533 if r == 3: # 外接
534 t = c1.r / (c1.r + c2.r)
535 return (c1.p + (c2.p - c1.p) * t,)
536 d = abs(c1.p - c2.p)
537 if r == 1: # 内接
538 if GeometryUtil.lt(c2.r, c1.r):
539 # if c2.r < c1.r - GeometryUtil.EPS:
540 return (c1.p + (c2.p - c1.p) * (c1.r / d),)
541 else:
542 return (c2.p + (c1.p - c2.p) * (c2.r / d),)
543 # 交点が2つ
544 rcos = (c1.r * c1.r + d * d - c2.r * c2.r) / (2 * d)
545 rsin = GeometryUtil.sqrt(c1.r * c1.r - rcos * rcos)
546 if GeometryUtil.lt(c1.r - abs(rcos), 0):
547 # if c1.r - abs(rcos) < GeometryUtil.EPS:
548 rsin = 0
549 e12 = (c2.p - c1.p) / abs(c2.p - c1.p)
550 return tuple(
551 sorted(
552 [
553 c1.p + rcos * e12 + rsin * Point(-e12.y, e12.x),
554 c1.p + rcos * e12 + rsin * Point(e12.y, -e12.x),
555 ]
556 )
557 )
558
[docs]
559 @classmethod
560 def cross_point_circle_line(cls, c: Circle, l: Line) -> list[Point]:
561 res = []
562 d = cls.dist_p_to_line(c.p, l)
563 if d > c.r + GeometryUtil.EPS:
564 return res
565 h = cls.projection_point(c.p, l)
566 if GeometryUtil.eq(d, c.r):
567 res.append(h)
568 return res
569 e = (l.p2 - l.p1).unit()
570 ph = GeometryUtil.sqrt(c.r * c.r - d * d)
571 res.append(h - e * ph)
572 res.append(h + e * ph)
573 res.sort()
574 return res
575
[docs]
576 @classmethod
577 def dist_p_to_line(cls, p: Point, l: Line) -> float:
578 return abs(cls.cross(l.p2 - l.p1, p - l.p1)) / abs(l.p2 - l.p1)
579
[docs]
580 @classmethod
581 def dist_p_to_p(cls, p1: Point, p2: Point) -> float:
582 return ((p1.x - p2.x) ** 2 + (p1.y - p2.y) ** 2) ** 0.5
583
[docs]
584 @classmethod
585 def dist_p_to_segmemt(cls, p: Point, s: Segment) -> float:
586 if cls.dot(s.p2 - s.p1, p - s.p1) < GeometryUtil.EPS:
587 return abs(p - s.p1)
588 if cls.dot(s.p1 - s.p2, p - s.p2) < GeometryUtil.EPS:
589 return abs(p - s.p2)
590 q = cls.projection_point(p, Line.from_points(s.p1, s.p2))
591 if (min(s.p1.x, s.p2.x) <= q.x <= max(s.p1.x, s.p2.x)) or (
592 min(s.p1.y, s.p2.y) <= q.y <= max(s.p1.y, s.p2.y)
593 ):
594 return abs(cls.cross(s.p2 - s.p1, p - s.p1)) / abs(s.p2 - s.p1)
595 return min(cls.dist_p_to_p(p, s.p1), cls.dist_p_to_p(p, s.p2))
596
[docs]
597 @classmethod
598 def dist_segment_to_segment(cls, s1: Segment, s2: Segment) -> float:
599 if cls.is_intersect_linesegment(s1, s2):
600 return 0.0
601 ans = GeometryUtil.GEO_INF
602 ans = min(ans, cls.dist_p_to_segmemt(s1.p1, s2))
603 ans = min(ans, cls.dist_p_to_segmemt(s1.p2, s2))
604 ans = min(ans, cls.dist_p_to_segmemt(s2.p1, s1))
605 ans = min(ans, cls.dist_p_to_segmemt(s2.p2, s1))
606 return ans
607
[docs]
608 @classmethod
609 def triangle_incircle(cls, p1: Point, p2: Point, p3: Point) -> Circle:
610 a, b, c = abs(p2 - p3), abs(p1 - p3), abs(p1 - p2)
611 p = Point(a * p1.x + b * p2.x + c * p3.x, a * p1.y + b * p2.y + c * p3.y)
612 p /= a + b + c
613 d = cls.dist_p_to_segmemt(p, Line.from_points(p1, p2))
614 return Circle(p, d)
615
[docs]
616 @classmethod
617 def triangle_circumcircle(cls, p1: Point, p2: Point, p3: Point) -> Circle:
618 a = 2 * (p1.x - p2.x)
619 b = 2 * (p1.y - p2.y)
620 p = p1.x**2 - p2.x**2 + p1.y**2 - p2.y**2
621 c = 2 * (p1.x - p3.x)
622 d = 2 * (p1.y - p3.y)
623 q = p1.x**2 - p3.x**2 + p1.y**2 - p3.y**2
624 u = a * d - b * c
625 x = d * p - b * q
626 y = a * q - c * p
627 if GeometryUtil.lt(u, 0):
628 x = -x
629 y = -y
630 u = -u
631 x /= u
632 y /= u
633 r = GeometryUtil.sqrt((x - p1.x) ** 2 + (y - p1.y) ** 2)
634 return Circle(Point(x, y), r)
635
[docs]
636 @classmethod
637 def triangle_excircle(
638 cls, p1: Point, p2: Point, p3: Point
639 ) -> tuple[Circle, Circle, Circle]:
640 dx1 = p2.x - p1.x
641 dy1 = p2.y - p1.y
642 dx2 = p3.x - p1.x
643 dy2 = p3.y - p1.y
644 d1 = GeometryUtil.sqrt((p3.x - p2.x) ** 2 + (p3.y - p2.y) ** 2)
645 d3 = GeometryUtil.sqrt(dx1**2 + dy1**2)
646 d2 = GeometryUtil.sqrt(dx2**2 + dy2**2)
647 S2 = abs(dx1 * dy2 - dx2 * dy1)
648 dsum1 = -d1 + d2 + d3
649 r1 = S2 / dsum1
650 ex1 = (-p1.x * d1 + p2.x * d2 + p3.x * d3) / dsum1
651 ey1 = (-p1.y * d1 + p2.y * d2 + p3.y * d3) / dsum1
652
653 dsum2 = d1 - d2 + d3
654 r2 = S2 / dsum2
655 ex2 = (p1.x * d1 - p2.x * d2 + p3.x * d3) / dsum2
656 ey2 = (p1.y * d1 - p2.y * d2 + p3.y * d3) / dsum2
657
658 dsum3 = d1 + d2 - d3
659 r3 = S2 / dsum3
660 ex3 = (p1.x * d1 + p2.x * d2 - p3.x * d3) / dsum3
661 ey3 = (p1.y * d1 + p2.y * d2 - p3.y * d3) / dsum3
662
663 return (
664 Circle(Point(ex1, ey1), r1),
665 Circle(Point(ex2, ey2), r2),
666 Circle(Point(ex3, ey3), r3),
667 )
668
[docs]
669 @classmethod
670 def circle_tangent(cls, p: Point, c: Circle) -> tuple[Point, Point]:
671 return cls.cross_point_circle(
672 c, Circle(p, GeometryUtil.sqrt((c.p - p).norm2() - c.r * c.r))
673 )
674
[docs]
675 @classmethod
676 def circle_common_tangent(cls, c1: Circle, c2: Circle) -> list[Line]:
677 ret = []
678 d = abs(c1.p - c2.p)
679 if GeometryUtil.eq(d, 0):
680 return ret
681 u = (c2.p - c1.p).unit()
682 v = u.rotate(GeometryUtil.pi / 2)
683 for s in [-1, 1]:
684 h = (c1.r + c2.r * s) / d
685 if GeometryUtil.eq(h * h, 1):
686 ret.append(
687 Line.from_points(
688 c1.p + (u if h > 0 else -u) * c1.r,
689 c1.p + (u if h > 0 else -u) * c1.r + v,
690 )
691 )
692 elif 1 - h * h > 0:
693 U = u * h
694 V = v * GeometryUtil.sqrt(1 - h * h)
695 ret.append(
696 Line.from_points(c1.p + (U + V) * c1.r, c2.p - (U + V) * (c2.r * s))
697 )
698 ret.append(
699 Line.from_points(c1.p + (U - V) * c1.r, c2.p - (U - V) * (c2.r * s))
700 )
701 return ret
702
[docs]
703 @classmethod
704 def convex_hull(cls, ps: list[Point], line_contains: bool = False) -> list[Point]:
705 """凸包を求める / `O(n)`
706
707 Args:
708 line_contains (bool, optional): 一直線上の点を含める場合、 `True` .
709
710 Returns:
711 list[Point]: 凸包
712 """
713 ps = cls.sort_by_y(ps)
714 if not line_contains:
715 ps = cls.unique(ps)
716 if len(ps) <= 1:
717 return ps
718
719 def cross(a, p):
720 if line_contains:
721 return cls.cross(a[-1] - a[-2], p - a[-1]) < -GeometryUtil.EPS
722 else:
723 return cls.cross(a[-1] - a[-2], p - a[-1]) < GeometryUtil.EPS
724
725 lower = []
726 for p in ps:
727 while len(lower) > 1 and cross(lower, p):
728 lower.pop()
729 lower.append(p)
730
731 upper = []
732 for p in reversed(ps):
733 while len(upper) > 1 and cross(upper, p):
734 upper.pop()
735 upper.append(p)
736
737 return lower[:-1] + upper[:-1]
738
[docs]
739 @classmethod
740 def unique(cls, a: list[Point]) -> list[Point]:
741 """同一要素の削除
742
743 Note:
744 `a` は事前にソートされていなければいけません。
745 """
746 if not a:
747 return []
748 res = [a[0]]
749 for i in range(1, len(a)):
750 if a[i] == res[-1]:
751 continue
752 res.append(a[i])
753 return res
754
[docs]
755 @classmethod
756 def argumant_sort(cls, a: list[Point]) -> list[Point]:
757 upper = []
758 lower = []
759 zeros = []
760 for p in a:
761 if GeometryUtil.lt(p.y, 0) or (
762 GeometryUtil.eq(p.y, 0) and GeometryUtil.lt(p.x, 0)
763 ):
764 lower.append(p)
765 elif GeometryUtil.eq(p.x, 0) and GeometryUtil.eq(p.y, 0):
766 zeros.append(p)
767 else:
768 upper.append(p)
769
770 def cmp(s: Point, t: Point) -> bool:
771 return GeometryUtil.lt(s.y * t.x, s.x * t.y) or GeometryUtil.eq(
772 s.y * t.x, s.x * t.y
773 )
774
775 upper = merge_sort(upper, key=cmp)
776 lower = merge_sort(lower, key=cmp)
777 return lower + zeros + upper
778
[docs]
779 @classmethod
780 def sort_by_x(cls, a: list[Point]) -> list[Point]:
781 return merge_sort(a)
782
783 @classmethod
784 def _cmp_y(cls, u: Point, v: Point):
785 if GeometryUtil.eq(u, v):
786 return True
787 if GeometryUtil.lt(u.y, v.y):
788 return True
789 if GeometryUtil.lt(v.y, u.y):
790 return False
791 return GeometryUtil.lt(u.x, v.x)
792
[docs]
793 @classmethod
794 def sort_by_y(cls, a: list[Point]) -> list[Point]:
795 return merge_sort(a, cls._cmp_y)
796
[docs]
797 @classmethod
798 def furthest_pair(cls, a: list[Point]) -> tuple[float, int, int]:
799 p = ConvexPolygon(a)
800 d, up, vp = p.diameter()
801 ui, vi = -1, -1
802 for i, point in enumerate(a):
803 if point == up:
804 ui = i
805 for i, point in enumerate(a):
806 if i != ui and point == vp:
807 vi = i
808 assert ui != -1 and vi != -1
809 return d, ui, vi
810
[docs]
811 @classmethod
812 def closest_pair(cls, a: list[Point]) -> tuple[float, int, int]:
813 """最近点対を求める / `O(nlogn)`
814
815 Args:
816 a (list[Point]): 距離配列
817
818 Returns:
819 tuple[float, int, int]: 距離、最近点対
820 """
821
822 buff = [None] * len(a)
823 d, p, q = GeometryUtil.GEO_INF, -1, -1
824
825 def f(l: int, r: int) -> None:
826 nonlocal d, p, q
827 mid = (l + r) // 2
828 x = b[mid][0].x
829 if mid - l > 1:
830 f(l, mid)
831 if r - mid > 1:
832 f(mid, r)
833 i, j = l, mid
834 indx = 0
835 while i < mid and j < r:
836 if cls._cmp_y(b[i][0], b[j][0]):
837 buff[indx] = b[i]
838 i += 1
839 else:
840 buff[indx] = b[j]
841 j += 1
842 indx += 1
843 for k in range(i, mid):
844 b[l + indx + k - i] = b[k]
845 for k in range(l, l + indx):
846 b[k] = buff[k - l]
847
848 near_line = []
849 for i in range(l, r):
850 e_p, e_i = b[i]
851 v = abs(e_p.x - x)
852 if cls._gt(v, d) or GeometryUtil.eq(v, d):
853 continue
854 for j in range(len(near_line) - 1, -1, -1):
855 ne_p, ne_i = near_line[j]
856 dx = e_p.x - ne_p.x
857 dy = e_p.y - ne_p.y
858 if cls._gt(dy, d) or GeometryUtil.eq(dy, d):
859 break
860 cand = GeometryUtil.sqrt(dx * dx + dy * dy)
861 if GeometryUtil.lt(cand, d):
862 d, p, q = cand, ne_i, e_i
863 near_line.append(b[i])
864
865 b = merge_sort((e, i) for i, e in enumerate(a))
866 f(0, len(b))
867 return d, p, q
868
[docs]
869 @classmethod
870 def count_segment_intersection(cls, a: list[Segment]) -> int:
871 """線分の交点の個数"""
872 # from titan_pylib.data_structures.fenwick_tree.fenwick_tree import FenwickTree
873 segs = []
874 events = []
875 Y = []
876 for i, seg in enumerate(a):
877 p1 = seg.p1
878 p2 = seg.p2
879 if p1.x == p2.x:
880 if p1.y > p2.y:
881 p1, p2 = p2, p1
882 segs.append((p1.y, p2.y))
883 events.append((p1.x, 1, i))
884 Y.append(p1.y)
885 Y.append(p2.y)
886 else:
887 if p1.x > p2.x:
888 p1, p2 = p2, p1
889 segs.append((p1.y, p2.y))
890 events.append((p1.x, 0, i))
891 events.append((p2.x, 2, i))
892 Y.append(p1.y)
893
894 to_origin = sorted(set(Y))
895 to_zaatsu = {a: i for i, a in enumerate(to_origin)}
896
897 bt = FenwickTree(len(to_origin))
898
899 ans = 0
900 events.sort()
901 for _, q, i in events:
902 if q == 0:
903 y1, y2 = segs[i]
904 bt.add(to_zaatsu[y1], 1)
905 elif q == 1:
906 y1, y2 = segs[i]
907 ans += bt.sum(to_zaatsu[y1], to_zaatsu[y2] + 1)
908 else:
909 y1, y2 = segs[i]
910 bt.add(to_zaatsu[y1], -1)
911 return ans