====== 기초 계산기하 ====== * 계산기하 기초 * 참고 * [[https://zigui.tistory.com/34]] * [[https://mathsciforstudent.tistory.com/206]] * 계산기하의 어려운 점 * 실수형을 반드시 사용해야 하는 경우가 많고, 그로 인해 실수 오차를 고려해야 한다 * 예외처리가 복잡하다. ===== 구현 유의사항 ===== ==== 실수 오차 줄이기 ==== * 우선 알아야 할것은 실수는 언제나 오차를 유발한다. * 파이썬은 기본적으로 실수형을 저장하는데에 64-bit 부동소숫점 표현을 사용한다. * [[https://docs.python.org/3/tutorial/floatingpoint.html]] * IEEE-754 “double precision" 을 따른다. * 이 방법은 소수부를 2진수로 근사해서 표현하기 때문에 작은 자릿수에서도 오차가 생길수 있다. * 예를 들어, 1.1 + 0.1 == 1.2 라는 식은 False를 리턴한다 * 1.1 + 0.1 == 1.2000000000000002 이다 * decimal 을 쓰는 것도 해결법이 되지는 않는다. 무한소수를 10진법으로 정확히 표현하는 것이 불가능하기 때문이다 * Decimal(1)/Decimal(3) 은 부동소숫점으로 표현되는 1/3 보다는 훨씬 더 많은 소숫점 이하 자리수를 저장하지만, 그렇다고 Decimal(1)/Decimal(3) * 3 == 1 이 True가 되지는 않는다. * Decimal(1)/Decimal(3) * 3 == '0.9999999999999999999999999999' 이다. * fraction을 쓰는 것은? * 실수를 요구하는 연산이 나눗셈 연산만이라면 주어진다면 가능한 해법이 될수도 있다 (속도 저하는 차치하고..) * 하지만 제곱근 연산 등이 추가되면 역시 해결 불가능한 오차가 생긴다. * 실수 오차는 계산이 누적될수록 점점 커진다. * 답이 실수로 나오더라도 문제에서 주어지는 값은 기본적으로 정수이다. * 문제를 출제시에 기본적으로 권장되는 가이드이다. * 결국 우리가 할수 있는 것은, 가능한 정수형만으로 계산하고, 실수로 결과가 나오는 연산은 최대한 마지막 부분으로 미루는 것이다. * 예를 들면 * ccw 체크를 할때, atan2로 각도를 진짜로 계산해서 비교하는게 아니라, 외적값의 부호를 이용 * 선분 교차 여부를 확인할때, 교점의 좌표를 구해보는 방식 보다는 점들간의 방향성을 이용 * 선분의 길이를 비교해야 할때에는, 실수로 표현되는 진짜 길이 보다는, 정수로 표현되는 길이의 제곱을 비교한다. ==== 가독성 vs 속도 ==== * 계산기하 알고리즘에서는 빈번히 쓰이는 연산들이 많은데 (벡터곱, CCW 등), 이런 연산들의 구현이 아주 짧은것은 또 아니다. 그래서 가독성을 위해서는 보조함수로 뽑아내는것이 당연한 선택이다. 그렇지만, 함수호출에 시간이 많이 걸리는 파이썬의 특성상, 이렇게 빈번하게 사용되는 연산들을 매번 함수 호출로 처리하는 것은 너무 느리다. * 그래서 실제 코드는 가독성을 좀 포기해도, 보조함수 사용을 좀 줄이도록 하겠다. * 이 문서에서 작성하는 코드는 설명 목적이므로, 가독성을 위해서 보조 함수를 사용해서 구현한다 * 실제 라이브러리에 구현하는 함수는 보조함수 대신에 인라이닝시켜 구현하겠다. * 예를 들면 * point_to_line_distance 함수를 이 문서에서는 보조함수들을 이용해서 이렇게 썼지만 * def vector_from_points(p0: Point, p1: Point): return (p1[0] - p0[0], p1[1] - p0[1]) def cross_product(u: Vector, v: Vector) -> int: return u[0] * v[1] - u[1] * v[0] def signed_parallelogram_area(p: Point, q: Point, r: Point) -> int: u = vector_from_points(p, q) v = vector_from_points(p, r) return cross_product(u, v) def point_to_line_distance(p0: Point, p1: Point, q: Point) -> float: return abs(signed_parallelogram_area(p0, p1, q)) / math.dist(p0, p1) * 실제 라이브러리에서는 이렇게 쓸 것이다. * def point_to_line_distance(p0: Point, p1: Point, q: Point) -> float: area = (p1[0] - p0[0]) * (q[1] - p0[1]) - (p1[1] - p0[1]) * (q[0] - p0[0]) return abs(area) / math.dist(p0, p1) * 그러니 필요하면 [[ps:teflib:geometry]]의 구현을 함께 참고하자. ===== 기본 자료형 ===== * 모든 자료형은 2차원 기준이다. * 3차원 이상의 좌표계에서의 기하학은 다른 토픽에서 다룬다 ==== 점 ==== * 대부분의 문제는 (x좌표, y좌표)를 입력으로 준다. 중간 계산 과정에서도 이 좌표계를 사용하는게 최선이다. * 대부분의 문제는 정수 좌표만을 사용해서 풀린다. tuple[int, int] 가 가장 자연스러운 선택이지만, 입력값을 point = [int(x) for x in input().split()] 처럼 파싱해서 list[int]로 받았을때 이것을 다시 tuple로 변환하는 것은 번거롭다. 두개를 다 허용하도록 하자. 타입 힌트는 Sequence[int]로 주면 된다. * Point: TypeAlias = Sequence[int] * 2차원에서의 점이라는 것을 분명히 하려고 Point2D 등으로 타입 이름을 붙이는 경우도 많다. 하지만 우리는 2차원 기하를 기본으로 다루기로 했으니까 2D라는 접미어는 쓰지 않겠다. (나중에 3차원 점이 필요하면 그때는 Point3D라고 쓰겠다). * (x, y)를 x+yi의 형태의 복소좌표로 간주해서, complex 타입을 사용해서 표현하는 트릭이 있다. * 점과 벡터를 모두 complex 타입으로 표현하면, 합과 차, 스칼라곱등을 따로 구현할 필요 없이 이미 built-in으로 정의된 연산을 이용해서 바로 사용할 수 있다. * 좀더 엄밀하게 말하면, 점과 벡터를 complex 타입으로 표현했다기 보다는, 벡터를 complex 타입으로 쓰고 점을 원점에서 출발하는 벡터로 표현했다는게 더 정확하긴 하다. * cpp에서는 실수부와 허수부의 값을 저장할 타입을 템플릿 형태로 지정할수 있고, 그래서 complex로 점과 벡터를 표현하는 것이 실제로 유용한 트릭이다 * 하지만 Python의 complex는 실수부와 허수부의 타입이 float로 고정된다. * 속도면에서는 문제가 없다. 원래 실수 연산은 정수 연산보다 느리지만, complex타입끼리의 built-in 연산이 list[int]에 대해서 직접 구현한 연산보다 더 빨랐다. * 문제는 정확도. complex에 정수 좌표만을 저장할 경우, 2^53 범위의 좌표까지만 정확하게 저장 가능하다. * 에러없이 표현가능한 유효숫자의 갯수는 15개이다 (sys.float_info.dig 으로 확인가능하다) * 그러므로 만약 32bit 정수형 (cpp기준 int)만 사용해서 모든 중간 결과를 다 계산할수 있는 문제라면, complex를 이용해도 상관 없기는 하다. 그렇지만 64bit 정수형 (cpp기준 long long)을 사용해서 풀어야만 하는 문제라면, complex로는 오차가 생길 여지가 있다. * complex를 사용하는 라이브러리와 Sequence[int]를 사용하는 라이브러리를 모두 만들어놓고, 문제에서 주어진 수의 범위에 따라서 골라서 쓰는것이 불가능한 아이디어는 아니지만 너무 에너지 낭비이다. * 그래서 python에서는 complex는 버리고 Sequence[int]로 통일시키는 것으로 결정했다. ==== 벡터 ==== * 2차원 벡터는 (x성분, y성분)으로 표현된다. 점과 동일한 차원이고, 같은 자료형으로 쓰면 되긴 한다. * 하지만, 그냥 똑같이 int의 페어로 쓰다보면 점과 개념적으로 헷갈릴때가 있는데 그래서 type alias를 구분해서 붙여두는것이 좋다 * Vector: TypeAlias = Sequence[int] * 점 p0에서 p1로 향하는 벡터는 두 점의 x좌표값의 차이와, y좌표값의 차이를 이용해서 구하면 된다. * def vector_from_points(p0: Point, p1: Point): return (p1[0] - p0[0], p1[1] - p0[1]) * 벡터를 complex 타입으로 쓸수 있었다면, p1 - p0로 쓰면 끝났을 텐데.. 아쉽다 ===== 벡터의 활용 ===== ==== 내적 (Inner Product) ==== * 두 벡터의 사이각 ==== 벡터곱 (Cross Product) ==== * 참고 * [[wp>ko:벡터곱]] * 원래 영문 용어는 cross product이고, 보통은 외적 또는 벡터곱이라고 부른다. 그러나 외적은 Outer product(텐서곱)을 의미하는 경우도 있기 때문에, 보다 명확히 하기 위해서 Cross product를 외적 대신에 벡터곱이라고 부르는 걸로 하겠다. * 원래 벡터곱은 3차원 벡터간에 정의되는 연산이다. 피연산자로 두개의 3차원벡터를 받아서, 연산결과로 3차원벡터가 나오는 연산 * 그런데도 계산기하에서는 2차원 벡터끼리 벡터곱을 취해서 스칼라값을 얻는 식으로 많이 사용되는데, 실제로는 z성분이 0인 3차원 벡터로 생각해서 연산하고, 연산결과의 z성분을 취해서 돌려주는 것이다. (x,y 성분은 모두 0이 된다) * def cross_product(u: Vector, v: Vector) -> int: return u[0] * v[1] - u[1] * v[0] * 이렇게 계산 결과로 얻어지는 값을 분석하면, * 부호는 두 벡터의 방향을 알려준다. * u×v>0 이면, u에서 v로의 회전이 반시계방향이다. (0,0)과 (ux, uy)를 잇는 직선을 기준으로 (vx, vy)가 왼쪽(위쪽)에 있다 * u×v<0 이면, u에서 v로의 회전이 시계방향이다. (0,0)과 (ux, uy)를 잇는 직선을 기준으로 (vx, vy)가 오른쪽(아래쪽)에 있다 * u×v=0 이면, u와 v는 평행해다 * 크기는, 두 벡터로 만들어지는 평행사변형의 넓이이다. |w| = |u||v|sinθ * 벡터곱의 부호와 크기, 두가지 모두 계산기하의 많은 알고리즘에 기본으로 활용된다. 게다가 입력 벡터값들이 모두 정수면, 연산결과도 정수형이어서 '실수를 쓰지 않고 계산하기'에 매우 적합하다. ==== 세 점의 방향 관계 ==== - 난이도: 골드 5 \\ - 기본 문제: [[ps:problems:boj:11758]] - 사전 지식: [[#벡터곱]] * 세 점 p,q,r이 시계방향(Clockwise; CW)으로 놓여있는지, 반시계방향(CounterClockwise; CCW)으로 놓여있는지를 판별하는 문제이다 * 흔히 그냥 CCW 이라고도 부른다. 인터넷에서 본 풀이에 'CCW을 써서 구한다' 이런 식의 표현이 있다면, 벡터곱을 이용해서 방향성을 확인한다는 의미라고 이해하면 된다. * p에서 q로 이동한 다음 r로 이동하려고 할때에 왼쪽 방향으로 회전해야 한다면 CCW, 오른쪽 방향으로 회전해야 한다면 CW 이다. * 계산하는 방법은 [[#벡터곱]]을 활용하면 간단하다. pq 벡터와 pr벡터의 벡터곱을 구한 뒤에, 값의 부호를 체크하면 방향성을 바로 판단할수 있다. * pq벡터와 pr벡터 대신에, pq벡터와 qr벡터의 벡터곱을 사용해도 결과는 동일하다. 이렇게 설명하는 자료도 많다. 내 경우는 pq와 pr로 생각하는 것이 더 이해하기 편해서 이쪽으로 사용했다. * def signed_parallelogram_area(p: Point, q: Point, r: Point) -> int: u = vector_from_points(p, q) v = vector_from_points(p, r) return cross_product(u, v) def is_ccw(p: Point, q: Point, r: Point) -> bool: return signed_parallelogram_area(p, q, r) > 0 def is_cw(p: Point, q: Point, r: Point) -> bool: return signed_parallelogram_area(p, q, r) < 0 def is_collinear(p: Point, q: Point, r: Point) -> bool: return signed_parallelogram_area(p, q, r) == 0 * 보통 계산기하의 입문으로 배우는 알고리즘이다. 이 방법을 모르는 상태에서는, 선분들의 실제 각도를 삼각함수 등을 써서 직접 구해서 비교하려다가, 실수 오차와 수직선등의 예외처리에 고생하기 쉽상이다. ==== 삼각형의 넓이 ==== * p, q, r을 꼭짓점으로 하는 삼각형의 넓이. * 위에서 말했듯 pq 벡터와 pr벡터의 벡터곱에 절댓값을 씌우면, pq와 pr로 만들어지는 평행사변형의 넓이가 나온다. 이 값을 2로 나누면 삼각형의 넓이이다. * def triangle_area(p: Point, q: Point, r: Point) -> float: return abs(signed_parallelogram_area(p, q, r)) / 2 ==== 점과 직선 사이의 거리 ==== * p0, p1을 연결하는 직선과 점 q와의 수직 거리 * 벡터를 배우기 전인 중학수학에서는 직선의 방정식을 ax+by+c=0 꼴로 표현했을 때의 점과 직선과의 거리공식을 유도해서 사용하지만. 벡터와 벡터곱을 이용하면 훨씬 간단하게 계산 가능하다. * def point_to_line_distance(p0: Point, p1: Point, q: Point) -> float return abs(signed_parallelogram_area(p0, p1, q)) / math.dist(p0, p1) ==== 점과 선분 사이의 최단 거리 ==== * p1, p2을 연결하는 선분과 점 q와의 최단 거리. * 점 q에서 선분 (p1, p2) 에 내린 수선의 발이 선분 안에 있는 지를 먼저 확인해야 한다. * 수선의 발이 선분 안에 있을 경우에는 그 수선의 길이가 최단 거리이고, 이 값은 [[#점과 직선 사이의 거리]]를 구한 것과 동일하게 구할수 있다 * 수선의 발이 선분 밖에 있을 경우에는 q에서 p1까지의 길이와 q에서 p2까지의 길이 중에서 더 짧은 쪽이 최단 거리가 된다. * 점 q에서 선분 (p1, p2) 에 내린 수선의 발이 선분 안에 있는지를 확인하기 위해서는, 각 q-p1-p2 와 각 q-p2-p1 이 모두 90도 보다 작은지를 확인하면 된다. * 실제로는 p1=>q 벡터와 p1=>p2 벡터의 내적이 0보다 크고, p2=>q 벡터와 p2=>p1 벡터의 내적이 0보다 큰지 확인하면 된다.