geometry

ソースコード

from titan_pylib.geometry.geometry import GeometryUtil
from titan_pylib.geometry.geometry import Point
from titan_pylib.geometry.geometry import Line
from titan_pylib.geometry.geometry import Polygon
from titan_pylib.geometry.geometry import ConvexPolygon
from titan_pylib.geometry.geometry import Circle
from titan_pylib.geometry.geometry import Segment
from titan_pylib.geometry.geometry import Geometry

view on github

展開済みコード

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

仕様

class Circle(p: Point, r: float)[source]

Bases: object

area() → float[source]

面積を返す

p
r
class ConvexPolygon(ps: list[Point], line_contains: bool = False)[source]

Bases: Polygon

contains(p: Point) → int[source]

点の包含関係を返す / \(O(\log{n})\)

Returns:

2: p を含む

1: p が辺上にある 0: それ以外

Return type:

int

convex_cut(l: Line) → ConvexPolygon[source]

直線 l で切断したときの左側の凸多角形を返す

diameter() → tuple[float, int, int][source]
class Geometry[source]

Bases: object

ref: - https://hcpc-hokudai.github.io/arpsive/geometry_004.pdf - https://bakamono1357.hatenablog.com/entry/2020/04/29/025320 - https://tjkendev.github.io/procon-library/python/geometry/circles_associated_with_triangle.html - https://atcoder.jp/contests/abc296/editorial/6109

classmethod argumant_sort(a: list[Point]) → list[Point][source]
classmethod ccw(u: Point, v: Point, p: Point) → int[source]

u->vに対し、v->pの位置関係を求める

Returns:

+1: a->bに対し、b->cが半時計回りに進む

-1: a->bに対し、b->cが時計回りに進む +2: a->bに対し、b->cが逆方向に進む -2: a->bに対し、b->cが同一方向に進む ` 0`: cが線分a->b上にある

Return type:

int

classmethod circle_common_tangent(c1: Circle, c2: Circle) → list[Line][source]
classmethod circle_tangent(p: Point, c: Circle) → tuple[Point, Point][source]
classmethod closest_pair(a: list[Point]) → tuple[float, int, int][source]

最近点対を求める / O(nlogn)

Parameters:

a (list[Point]) – 距離配列

Returns:

距離、最近点対

Return type:

tuple[float, int, int]

classmethod convex_hull(ps: list[Point], line_contains: bool = False) → list[Point][source]

凸包を求める / O(n)

Parameters:

line_contains (bool, optional) – 一直線上の点を含める場合、 True .

Returns:

凸包

Return type:

list[Point]

classmethod count_segment_intersection(a: list[Segment]) → int[source]

線分の交点の個数

classmethod cross(p1: Point, p2: Point) → float[source]
classmethod cross_point_circle(c1: Circle, c2: Circle) → tuple[Point, ...][source]
classmethod cross_point_circle_line(c: Circle, l: Line) → list[Point][source]
classmethod cross_point_line(l1: Line, l2: Line) → Point | None[source]
classmethod cross_point_line_segment(l: Line, s: Segment) → Point | None[source]
classmethod cross_point_segment(s1: Segment, s2: Segment) → Point | None[source]
classmethod dist_p_to_line(p: Point, l: Line) → float[source]
classmethod dist_p_to_p(p1: Point, p2: Point) → float[source]
classmethod dist_p_to_segmemt(p: Point, s: Segment) → float[source]
classmethod dist_segment_to_segment(s1: Segment, s2: Segment) → float[source]
classmethod dot(p1: Point, p2: Point) → float[source]
classmethod furthest_pair(a: list[Point]) → tuple[float, int, int][source]
classmethod is_intersect_circle(c1: Circle, c2: Circle) → int[source]

2円の位置関係

Returns:

共通接線の数

0: 内包している 1: 内接している 2: 2点で交わっている 3: 外接している 4: 離れている

Return type:

int

classmethod is_intersect_linesegment(s1: Segment, s2: Segment) → bool[source]

線分 s1 と s2 が交差しているかどうか判定する

classmethod is_orthogonal(l1: Line, l2: Line) → bool[source]

直線 l1, l2 が直行しているかどうか

classmethod is_parallel(l1: Line, l2: Line) → bool[source]

直線 l1, l2 が平行かどうか

classmethod projection_point(p: Point, l: Line) → Point[source]

直線 l に、点 p からおろした垂線の足の Point を返す

classmethod reflection_point(p: Point, l: Line) → Point[source]

直線 l を対象軸として点 p と線対称の点を返す

classmethod sort_by_x(a: list[Point]) → list[Point][source]
classmethod sort_by_y(a: list[Point]) → list[Point][source]
classmethod triangle_circumcircle(p1: Point, p2: Point, p3: Point) → Circle[source]
classmethod triangle_excircle(p1: Point, p2: Point, p3: Point) → tuple[Circle, Circle, Circle][source]
classmethod triangle_incircle(p1: Point, p2: Point, p3: Point) → Circle[source]
classmethod unique(a: list[Point]) → list[Point][source]

同一要素の削除

Note

a は事前にソートされていなければいけません。

class GeometryUtil[source]

Bases: object

EPS: Final[Decimal | float] = 1e-10
GEO_INF: Final[Decimal | float] = inf
USE_DECIMAL: Final[bool] = False
acos(x, /)

Return the arc cosine (measured in radians) of x.

The result is between 0 and pi.

cos(x, /)

Return the cosine of x (measured in radians).

classmethod eq(a: float, b: float) → bool[source]
classmethod ge(u: float, v: float) → bool[source]
classmethod gt(u: float, v: float) → bool[source]
classmethod le(u: float, v: float) → bool[source]
classmethod lt(u: float, v: float) → bool[source]
pi = 3.141592653589793
sin(x, /)

Return the sine of x (measured in radians).

classmethod sqrt(v: Decimal | float) → None[source]
class Line(a: int, b: int, c: int, p1: Point, p2: Point)[source]

Bases: object

ax + by + c = 0

a
b
c
classmethod from_abc(a: int, b: int, c: int) → Line[source]
classmethod from_points(p1: Point, p2: Point) → Line[source]
get_x(y: float) → float | None[source]
get_y(x: float) → float | None[source]
p1
p2
class Point(x: float, y: float)[source]

Bases: object

copy() → Point[source]
norm() → float[source]
norm2() → float[source]

(x**2 + y**2)を返す(ルートの計算はしない)

rotate(theta: float) → Point[source]
show(precision: int = 16) → None[source]
unit() → Point[source]
x
y
class Polygon(ps: list[Point])[source]

Bases: object

多角形クラス

area() → float[source]

面積を求める / \(O(n)\)

contains(p: Point) → int[source]

点の包含関係を返す / O(n)

Returns:

2: p を含む

1: p が辺上にある 0: それ以外

Return type:

int

get_degree(radian)[source]
is_convex() → bool[source]

凸多角形かどうか / \(O(n)\)

class Segment(p1: Point, p2: Point)[source]

Bases: object