최소 외접원¶

다음 문제를 생각해 봅시다:

Library Checker - Minimum Enclosing Circle

$n \leq 10^5$개의 점 $p_i=(x_i, y_i)$가 주어집니다.

각 $p_i$에 대해, 그 점이 $\{p_1,\dots,p_n\}$의 최소 외접원의 원주 위에 있는지 판별하십시오.

여기서 최소 외접원(MEC)이란 $n$개의 점을 모두 원의 내부 또는 경계 위에 포함하는 원 중 반지름이 가능한 한 가장 작은 원을 의미합니다. 이 문제에는 간단한 무작위 풀이가 있는데, 언뜻 보면 $O(n^3)$에 동작할 것처럼 보이지만 실제로는 기대 시간 $O(n)$에 동작합니다.

아래의 논증을 더 잘 이해하기 위해, 이 문제의 해가 유일하다는 점을 먼저 짚고 넘어가야 합니다:

MEC는 왜 유일합니까?

다음과 같은 설정을 생각해 봅시다: $r$을 MEC의 반지름이라고 합시다. 각 점 $p_1,\dots,p_n$을 중심으로 반지름이 $r$인 원을 그립니다. 기하학적으로, 반지름이 $r$이면서 모든 점 $p_1,\dots,p_n$을 덮는 원들의 중심은 $n$개의 원 전체의 교집합을 이룹니다.

이제 교집합이 단 하나의 점이라면, 이것만으로 해가 유일하다는 것이 증명됩니다. 그렇지 않다면 교집합은 넓이가 0이 아닌 도형이므로, $r$을 아주 조금 줄여도 여전히 교집합이 비어 있지 않게 할 수 있으며, 이는 $r$이 외접원의 가능한 최소 반지름이라는 가정에 모순됩니다.

비슷한 논리로, MEC가 주어진 특정 점 $p_i$ 또는 두 점 $p_i$와 $p_j$를 지나야 한다는 조건을 추가로 요구하는 경우에도 MEC가 유일함을 보일 수 있습니다(이 경우에도 반지름이 원을 유일하게 결정하므로 유일합니다).

또는 MEC가 두 개 있다고 가정한 뒤, 두 원의 교집합(점 $p_1,\dots,p_n$을 이미 포함하는 영역)은 처음 원들보다 지름이 작아야 하므로 더 작은 원으로 덮을 수 있다는 점에 주목할 수도 있습니다.

Welzl 알고리즘¶

간결하게 쓰기 위해 $\{p_1,\dots,p_n\}$의 MEC를 $\operatorname{mec}(p_1,\dots,p_n)$으로 표기하고, $P_i = \{p_1,\dots,p_i\}$라고 합시다.

1991년 Welzl이 처음 제안한 알고리즘은 다음과 같습니다:

입력으로 주어진 점들의 순서에 무작위 순열을 적용합니다.
MEC의 현재 후보 $C$를 유지하며, $C = \operatorname{mec}(p_1, p_2)$에서 시작합니다.
$i=3..n$에 대해 반복하며 $p_i \in C$인지 확인합니다.
1. $p_i \in C$이면 $C$가 $P_i$의 MEC라는 뜻입니다.
2. 그렇지 않으면 $C = \operatorname{mec}(p_i, p_1)$로 두고, $j=2..i$에 대해 반복하며 $p_j \in C$인지 확인합니다.
  1. $p_j \in C$이면 $C$는 $p_i$를 지나는 원 중에서 $P_j$의 MEC입니다.
  2. 그렇지 않으면 $C=\operatorname{mec}(p_i, p_j)$로 두고, $k=1..j$에 대해 반복하며 $p_k \in C$인지 확인합니다.
    1. $p_k \in C$이면 $C$는 $p_i$와 $p_j$를 지나는 원 중에서 $P_k$의 MEC입니다.
    2. 그렇지 않으면 $C=\operatorname{mec}(p_i,p_j,p_k)$가 $p_i$와 $p_j$를 지나는 원 중에서 $P_k$의 MEC입니다.

여기서 각 중첩 단계마다 유지해야 할 불변식($C$가 추가로 주어진 $0$개, $1$개 또는 $2$개의 점도 지나는 원 중에서 MEC라는 것)이 있으며, 안쪽 루프가 끝날 때마다 그 불변식은 부모 루프의 현재 반복에서의 불변식과 동치가 됨을 알 수 있습니다. 이는 결과적으로 알고리즘 전체의 정확성을 보장합니다.

일부 기술적인 세부 사항은 일단 생략하면, 전체 알고리즘은 C++로 다음과 같이 구현할 수 있습니다:

struct point {...};

// Is represented by 2 or 3 points on its circumference
struct mec {...};

bool inside(mec const& C, point p) {
    return ...;
}

// Choose some good generator of randomness for the shuffle
mt19937_64 gen(...);
mec enclosing_circle(vector<point> &p) {
    int n = p.size();
    ranges::shuffle(p, gen);
    auto C = mec{p[0], p[1]};
    for(int i = 0; i < n; i++) {
        if(!inside(C, p[i])) {
            C = mec{p[i], p[0]};
            for(int j = 0; j < i; j++) {
                if(!inside(C, p[j])) {
                    C = mec{p[i], p[j]};
                    for(int k = 0; k < j; k++) {
                        if(!inside(C, p[k])) {
                            C = mec{p[i], p[j], p[k]};
                        }
                    }
                }
            }
        }
    }
    return C;
}

점 $p_i$가 $2$개 또는 $3$개 점의 MEC 내부에 있는지 확인하는 작업은 $O(1)$에 할 수 있으리라 예상할 수 있습니다(이에 대해서는 뒤에서 다룹니다). 그렇더라도 위 알고리즘은 중첩된 루프들 때문에 최악의 경우 $O(n^3)$이 걸리는 것처럼 보입니다. 그렇다면 어떻게 기대 실행 시간이 선형이라고 주장할 수 있었을까요? 함께 알아봅시다!

복잡도 분석¶

가장 안쪽 루프($k$에 대한 루프)의 기대 실행 시간은 분명히 $O(j)$번의 연산입니다. 그렇다면 $j$에 대한 루프는 어떻습니까?

이 루프는 $p_j$가 점 $i$도 지나는 $P_j$의 MEC 경계 위에 있고, $p_j$를 제거하면 원이 더 작아지는 경우에만 다음 루프를 실행합니다. $P_j$의 모든 점 중에서 이러한 성질을 가진 점은 최대 $2$개뿐입니다. 경계 위에 $P_j$의 점이 $2$개보다 많이 있다면, 그중 어느 점을 제거하더라도 경계 위에 여전히 최소 $3$개의 점이 남아 원을 유일하게 결정하기에 충분하기 때문입니다.

다시 말해, 처음에 무작위로 섞은 후 최대 두 개인 운 나쁜 점 중 하나를 $p_j$로 뽑을 확률은 최대 $\frac{2}{j}$입니다. 이를 $1$부터 $i$까지의 모든 $j$에 대해 합하면, 기대 실행 시간은 다음과 같습니다.

$$ \sum\limits_{j=1}^i \frac{2}{j} \cdot O(j) = O(i). $$

이제 정확히 같은 방식으로 가장 바깥쪽 루프의 기대 실행 시간이 $O(n)$임도 증명할 수 있습니다.

점이 2개 또는 3개 점의 MEC 내부에 있는지 확인하기¶

이제 point와 mec의 구현 세부 사항을 알아봅시다. 이 문제에서는 점을 나타내는 클래스로 std::complex를 사용하는 것이 특히 유용합니다.

using ftype = int64_t;
using point = complex<ftype>;

다시 상기하자면, 복소수는 $x+yi$ 꼴의 수이며, 여기서 $i^2=-1$이고 $x, y \in \mathbb R$입니다. C++에서 이러한 복소수는 2차원 점 $(x, y)$로 표현됩니다. 복소수는 기본적인 성분별 선형 연산(덧셈, 실수배)을 이미 구현하고 있을 뿐만 아니라, 복소수의 곱셈과 나눗셈에도 특정한 기하학적 의미가 있습니다.

너무 자세히 들어가지 않고, 이 문제에 가장 중요한 성질만 짚겠습니다. 두 복소수를 곱하면 편각($Ox$에서 반시계 방향으로 잰 각)이 더해지고, 켤레를 취하면(즉, $z=x+yi$를 $\overline{z} = x-yi$로 바꾸면) 편각에 $-1$이 곱해집니다. 이를 이용하면 점 $z$가 특정한 $2$개 또는 $3$개 점의 MEC 내부에 있는지에 대한 매우 간단한 판정 기준을 세울 수 있습니다.

2개 점의 MEC¶

$2$개의 점 $a$와 $b$에 대해, 이들의 MEC는 단순히 중심이 $\frac{a+b}{2}$이고 반지름이 $\frac{|a-b|}{2}$인 원, 다시 말해 $ab$를 지름으로 하는 원입니다. $z$가 이 원 내부에 있는지 확인하려면 $za$와 $zb$ 사이의 각이 예각이 아닌지만 확인하면 됩니다.

내부의 각은 둔각이고, 외부의 각은 예각이며, 원주 위의 각은 직각입니다

이는 곧 다음 값이

$$ I_0=(b-z)\overline{(a-z)} $$

양의 실수 좌표를 갖지 않는지(편각이 $-90^\circ$와 $90^\circ$ 사이인 점에 해당합니다) 확인하는 것과 같습니다.

3개 점의 MEC¶

삼각형 $abc$에 $z$를 추가하면 사각형이 됩니다. 다음 식을 생각해 봅시다.

$$ \angle azb + \angle bca $$

원에 내접하는 사각형에서 $c$와 $z$가 $ab$에 대해 같은 쪽에 있으면 두 각은 같으며, 부호를 고려하여(즉, 반시계 방향이면 양수, 시계 방향이면 음수로) 더하면 합이 $0^\circ$가 됩니다. 마찬가지로, $c$와 $z$가 서로 반대쪽에 있으면 두 각의 합은 $180^\circ$가 됩니다.

이웃한 원주각은 서로 같고, 마주 보는 각은 합이 180도가 됩니다

복소수의 관점에서 보면, $\angle azb$는 $(b-z)\overline{(a-z)}$의 편각이고 $\angle bca$는 $(a-c)\overline{(b-c)}$의 편각임을 알 수 있습니다. 따라서 $\angle azb + \angle bca$는 다음 값의 편각이라고 결론지을 수 있습니다

$$ I_1 = (b-z) \overline{(a-z)} (a-c) \overline{(b-c)} $$

각이 $0^\circ$ 또는 $180^\circ$라면 $I_1$의 허수부가 $0$이라는 뜻이고, 그렇지 않다면 $I_1$의 허수부의 부호를 확인하여 $z$가 $abc$의 외접원 안쪽에 있는지 바깥쪽에 있는지 알아낼 수 있습니다. 양의 허수부는 양의 각에 대응하고, 음의 허수부는 음의 각에 대응합니다.

그렇다면 이 중 어느 쪽이 $z$가 원의 안쪽에 있음을, 어느 쪽이 바깥쪽에 있음을 의미합니까? 앞서 살펴보았듯이, $z$가 원 안쪽에 있으면 일반적으로 $\angle azb$의 크기가 커지고, 원 바깥쪽에 있으면 크기가 작아집니다. 따라서 다음과 같은 4가지 경우가 있습니다:

$\angle bca > 0^\circ$이고, $c$가 $ab$를 기준으로 $z$와 같은 쪽에 있는 경우입니다. 그러면 원 안쪽의 점에 대해 $\angle azb < 0^\circ$이고, $\angle azb + \angle bca < 0^\circ$입니다.
$\angle bca < 0^\circ$이고, $c$가 $ab$를 기준으로 $z$와 같은 쪽에 있는 경우입니다. 그러면 원 안쪽의 점에 대해 $\angle azb > 0^\circ$이고, $\angle azb + \angle bca > 0^\circ$입니다.
$\angle bca > 0^\circ$이고, $c$가 $ab$를 기준으로 $z$의 반대쪽에 있는 경우입니다. 그러면 원 안쪽의 점에 대해 $\angle azb > 0^\circ$이고 $\angle azb + \angle bca > 180^\circ$입니다.
$\angle bca < 0^\circ$이고, $c$가 $ab$를 기준으로 $z$의 반대쪽에 있는 경우입니다. 그러면 원 안쪽의 점에 대해 $\angle azb < 0^\circ$이고 $\angle azb + \angle bca < 180^\circ$입니다.

다시 말해, 각을 $-180^\circ$에서 $180^\circ$ 사이로 정규화한다고 가정하면, $\angle bca$가 양수일 때 원 안쪽의 점은 $\angle azb + \angle bca < 0^\circ$를 만족하고, 그렇지 않을 때는 $\angle azb + \angle bca > 0^\circ$를 만족합니다. 이는 다시 $I_2=(a-c)\overline{(b-c)}$와 $I_1 = I_0 I_2$의 허수부 부호로 확인할 수 있습니다.

참고: $I_1$을 구하기 위해 복소수 네 개를 곱하므로 중간 계수는 $O(A^4)$만큼 커질 수 있습니다. 여기서 $A$는 입력에서 가장 큰 좌표의 절댓값입니다. 다행히 입력이 정수라면 위의 두 검사를 모두 정수 연산만으로 수행할 수 있습니다.

구현¶

이제 이 검사를 실제로 구현하려면, 먼저 MEC를 어떻게 표현할지 정해야 합니다. 우리의 판정 기준은 점을 직접 다루므로, 자연스럽고 효율적인 방법은 MEC를 그것을 정의하는 두 점 또는 세 점으로 직접 표현하는 것입니다:

using mec = variant<
    array<point, 2>,
    array<point, 3>
>;

이제 std::visit를 사용하여 위의 판정 기준에 따라 두 경우를 모두 효율적으로 처리할 수 있습니다:

/* I < 0 if z inside C,
   I > 0 if z outside C,
   I = 0 if z on the circumference of C */
ftype indicator(mec const& C, point z) {
    return visit([&](auto &&C) {
        point a = C[0], b = C[1];
        point I0 = (b - z) * conj(a - z);
        if constexpr (size(C) == 2) {
            return real(I0);
        } else {
            point c = C[2];
            point I2 = (a - c) * conj(b - c);
            point I1 = I0 * I2;
            return imag(I2) < 0 ? -imag(I1) : imag(I1);
        }
    }, C);
}

bool inside(mec const& C, point p) {
    return indicator(C, p) <= 0;
}

이제 마지막으로 Library Checker에 문제를 제출하여 모든 것이 올바르게 동작하는지 확인할 수 있습니다: #308668.