Мінімальне охоплююче коло

Розглянемо таку задачу:

інформація

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)$.

Коли підходить цей алгоритм?

• Потрібне коло мінімального радіуса, що містить усі задані точки (всередині або на межі)?
• Тобі потрібна охоплююча фігура у вигляді кола, а не опуклого многокутника? (якщо ні → Побудова опуклої оболонки)
• Прийнятний рандомізований алгоритм з очікуваним лінійним часом $O(n)$?

Щоб краще зрозуміти міркування нижче, нам варто одразу зауважити, що розв'язок задачі єдиний:

Чому MEC єдине?

Розглянемо таку конструкцію. Нехай $r$ — радіус MEC. Навколо кожної з точок $p_1,\dots,p_n$ намалюємо коло радіуса $r$. Геометрично центри кіл, що мають радіус $r$ і покривають усі точки $p_1,\dots,p_n$, утворюють перетин усіх $n$ кіл.

Тепер, якщо цей перетин — лише одна точка, це вже доводить, що MEC єдине. Інакше перетин — це фігура ненульової площі, тож ми можемо трохи зменшити $r$ і все одно мати непорожній перетин, що суперечить припущенню про те, що $r$ був мінімально можливим радіусом охоплюючого кола.

За схожою логікою ми можемо також показати єдиність MEC, якщо додатково вимагати, щоб воно проходило через задану конкретну точку $p_i$ або через дві точки $p_i$ і $p_j$ (воно теж єдине, бо його радіус однозначно його визначає).

Інакше можна також припустити, що існують два MEC, і тоді зауважити, що їхній перетин (який уже містить точки $p_1,\dots,p_n$) повинен мати менший діаметр, ніж початкові кола, а отже, його можна покрити меншим колом.

Алгоритм Вельцля

Для стислості позначимо через $\operatorname{mec}(p_1,\dots,p_n)$ мінімальне охоплююче коло множини $\{p_1,\dots,p_n\}$, і нехай $P_i = \{p_1,\dots,p_i\}$.

Алгоритм, який спершу запропонував Вельцль 1991 року, виглядає так:

1. Застосуємо випадкову перестановку до вхідної послідовності точок.
2. Підтримуємо поточного кандидата на MEC $C$, починаючи з $C = \operatorname{mec}(p_1, p_2)$.
3. Проходимо по $i=3..n$ і перевіряємо, чи $p_i \in C$.
   1. Якщо $p_i \in C$, це означає, що $C$ є MEC множини $P_i$.
   2. Інакше присвоюємо $C = \operatorname{mec}(p_i, p_1)$ і проходимо по $j=2..i$, перевіряючи, чи $p_j \in C$.
      1. Якщо $p_j \in C$, то $C$ є MEC множини $P_j$ серед кіл, що проходять через $p_i$.
      2. Інакше присвоюємо $C=\operatorname{mec}(p_i, p_j)$ і проходимо по $k=1..j$, перевіряючи, чи $p_k \in C$.
         1. Якщо $p_k \in C$, то $C$ є MEC множини $P_k$ серед кіл, що проходять через $p_i$ і $p_j$.
         2. Інакше $C=\operatorname{mec}(p_i,p_j,p_k)$ є MEC множини $P_k$ серед кіл, що проходять через $p_i$ і $p_j$.

Ми бачимо, що кожен рівень вкладеності тут має інваріант, який треба підтримувати (а саме те, що $C$ є MEC серед кіл, які додатково проходять ще й через задані $0$, $1$ або $2$ точки), і щоразу, коли внутрішній цикл завершується, його інваріант стає еквівалентним інваріанту поточної ітерації батьківського циклу. Це, своєю чергою, забезпечує коректність алгоритму в цілому.

Опускаючи поки що деякі технічні деталі, увесь алгоритм можна реалізувати на C++ так:

C++

struct point {...};

// Представляється 2 або 3 точками на його колі
struct mec {...};

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

// Виберіть якийсь хороший генератор випадковості для перемішування
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;
}

Python

# point — це ...
# mec представляється 2 або 3 точками на його колі

def inside(C, p):
    return ...

# Виберіть якийсь хороший генератор випадковості для перемішування
rng = random.Random(...)

def enclosing_circle(p):
    n = len(p)
    rng.shuffle(p)
    C = (p[0], p[1])
    for i in range(n):
        if not inside(C, p[i]):
            C = (p[i], p[0])
            for j in range(i):
                if not inside(C, p[j]):
                    C = (p[i], p[j])
                    for k in range(j):
                        if not inside(C, p[k]):
                            C = (p[i], p[j], p[k])
    return C

TypeScript

// Point — це ...
// Mec представляється масивом із 2 або 3 точок на його колі
type Mec = Point[];

function inside(C: Mec, p: Point): boolean {
    return ...;
}

// Виберіть якийсь хороший генератор випадковості для перемішування
const rnd = makeRng(...);

function enclosingCircle(p: Point[]): Mec {
    const n = p.length;
    shuffle(p, rnd);
    let C: Mec = [p[0], p[1]];
    for (let i = 0; i < n; i++) {
        if (!inside(C, p[i])) {
            C = [p[i], p[0]];
            for (let j = 0; j < i; j++) {
                if (!inside(C, p[j])) {
                    C = [p[i], p[j]];
                    for (let k = 0; k < j; k++) {
                        if (!inside(C, p[k])) {
                            C = [p[i], p[j], p[k]];
                        }
                    }
                }
            }
        }
    }
    return C;
}

Go

// point — це ...
// mec представляється зрізом із 2 або 3 точок на його колі
type mec = []point

func inside(C mec, p point) bool {
    return ...
}

// Виберіть якийсь хороший генератор випадковості для перемішування
func enclosingCircle(p []point, rng *rand.Rand) mec {
    n := len(p)
    rng.Shuffle(n, func(i, j int) { p[i], p[j] = p[j], p[i] })
    C := mec{p[0], p[1]}
    for i := 0; i < n; i++ {
        if !inside(C, p[i]) {
            C = mec{p[i], p[0]}
            for j := 0; j < i; j++ {
                if !inside(C, p[j]) {
                    C = mec{p[i], p[j]}
                    for k := 0; k < j; k++ {
                        if !inside(C, p[k]) {
                            C = mec{p[i], p[j], p[k]}
                        }
                    }
                }
            }
        }
    }
    return C
}

Тепер можна сподіватися, що перевірку того, чи точка $p_i$ лежить усередині MEC з $2$ або $3$ точок, можна виконати за $O(1)$ (це ми обговоримо пізніше). Але навіть тоді наведений вище алгоритм виглядає так, ніби він займатиме $O(n^3)$ у найгіршому випадку — лише через усі ці вкладені цикли. То чому ж ми стверджували, що очікуваний час роботи лінійний? Розберімося!

Аналіз складності

Для найвнутрішнішого циклу (по $k$) очевидно, що його очікуваний час роботи становить $O(j)$ операцій. А що з циклом по $j$?

Він запускає наступний цикл лише тоді, коли $p_j$ лежить на межі MEC множини $P_j$, яке також проходить через точку $i$, і видалення $p_j$ ще більше зменшило б коло. Серед усіх точок $P_j$ може бути щонайбільше $2$ точки з такою властивістю, бо якщо на межі є більш ніж $2$ точки з $P_j$, це означає, що після видалення будь-якої з них на межі все одно залишиться щонайменше $3$ точки, чого достатньо, щоб однозначно визначити коло.

Іншими словами, після початкового випадкового перемішування ймовірність того, що однією з щонайбільше двох невдалих точок виявиться $p_j$, не перевищує $\frac{2}{j}$. Підсумувавши це по всіх $j$ від $1$ до $i$, отримуємо очікуваний час роботи

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

Точно так само ми тепер можемо довести, що найзовнішніший цикл має очікуваний час роботи $O(n)$.

Перевірка того, чи точка лежить у MEC з 2 або 3 точок

Розберімося тепер з деталями реалізації point і mec. У цій задачі виявляється особливо зручним використовувати std::complex як клас для точок:

C++

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

Python

# У Python вбудований тип complex уже представляє двовимірну точку x + yi
# (його координати — числа з рухомою комою, але для цілих координат
# проміжні обчислення з порядком O(A^4) лишаються точними).
Point = complex

TypeScript

// У TypeScript комплексних чисел немає, тому реалізуємо точку як пару (x, y)
// і власні операції множення/спряження нижче.
type Point = { x: number; y: number };

Go

// У Go вбудований тип complex128 уже представляє двовимірну точку x + yi
type point = complex128

Нагадаємо, що комплексне число — це число виду $x+yi$, де $i^2=-1$ і $x, y \in \mathbb R$. У C++ таке комплексне число представляється двовимірною точкою $(x, y)$. Комплексні числа вже реалізують базові покомпонентні лінійні операції (додавання, множення на дійсне число), але також їх множення й ділення несуть певний геометричний зміст.

Не вдаючись надто в деталі, відзначимо найважливішу для цієї конкретної задачі властивість: множення двох комплексних чисел додає їхні полярні кути (відраховані від $Ox$ проти годинникової стрілки), а взяття спряженого числа (тобто перетворення $z=x+yi$ на $\overline{z} = x-yi$) множить полярний кут на $-1$. Це дозволяє нам сформулювати дуже прості критерії того, чи лежить точка $z$ усередині MEC з $2$ або $3$ конкретних точок.

MEC з 2 точок

Для $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$).

MEC з 3 точок

Додавання $z$ до трикутника $abc$ перетворить його на чотирикутник. Розглянемо такий вираз:

$$\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$; інакше ми можемо визначити, чи лежить $z$ усередині або зовні охоплюючого кола трикутника $abc$, перевіривши знак уявної частини $I_1$. Додатна уявна частина відповідає додатним кутам, а від'ємна уявна частина — від'ємним кутам.

Але яка з них означає, що $z$ лежить усередині, а яка — що зовні кола? Як ми вже зауважили, коли $z$ усередині кола, це загалом збільшує модуль $\angle azb$, тоді як коли воно зовні кола — зменшує. Тож маємо такі 4 випадки:

1. $\angle bca > 0^\circ$, $c$ з того самого боку від $ab$, що й $z$. Тоді $\angle azb < 0^\circ$, і $\angle azb + \angle bca < 0^\circ$ для точок усередині кола.
2. $\angle bca < 0^\circ$, $c$ з того самого боку від $ab$, що й $z$. Тоді $\angle azb > 0^\circ$, і $\angle azb + \angle bca > 0^\circ$ для точок усередині кола.
3. $\angle bca > 0^\circ$, $c$ з протилежного боку від $ab$ відносно $z$. Тоді $\angle azb > 0^\circ$, і $\angle azb + \angle bca > 180^\circ$ для точок усередині кола.
4. $\angle bca < 0^\circ$, $c$ з протилежного боку від $ab$ відносно $z$. Тоді $\angle azb < 0^\circ$, і $\angle azb + \angle bca < 180^\circ$ для точок усередині кола.

Іншими словами, якщо $\angle bca$ додатний, точки всередині кола матимуть $\angle azb + \angle bca < 0^\circ$, інакше вони матимуть $\angle azb + \angle bca > 0^\circ$, за умови, що ми нормалізуємо кути в межах від $-180^\circ$ до $180^\circ$. Це, своєю чергою, можна перевірити за знаками уявних частин $I_2=(a-c)\overline{(b-c)}$ і $I_1 = I_0 I_2$.

Зауваження: оскільки для отримання $I_1$ ми перемножуємо чотири комплексні числа, проміжні коефіцієнти можуть бути аж порядку $O(A^4)$, де $A$ — найбільший за модулем координат у вхідних даних. З хорошого боку, якщо вхідні дані цілочислові, обидві перевірки вище можна виконати повністю в цілих числах.

Реалізація

Тепер, щоб реально реалізувати перевірку, нам спершу слід вирішити, як представляти MEC. Оскільки наші критерії працюють безпосередньо з точками, природний та ефективний спосіб зробити це — сказати, що MEC безпосередньо представляється парою або трійкою точок, що його визначають:

C++

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

Python

# mec — це просто кортеж із 2 або 3 точок; розрізняємо випадки за len()
# (наприклад, (a, b) або (a, b, c))

TypeScript

// mec — це масив із 2 або 3 точок; розрізняємо випадки за C.length
type Mec = Point[];

Go

// mec — це зріз із 2 або 3 точок; розрізняємо випадки за len(C)
type mec = []point

Тепер ми можемо використати std::visit, щоб ефективно опрацювати обидва випадки відповідно до наведених вище критеріїв:

C++

/* I < 0, якщо z всередині C,
   I > 0, якщо z зовні C,
   I = 0, якщо z на колі 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;
}

Python

# I < 0, якщо z всередині C; I > 0, якщо z зовні C; I = 0, якщо z на колі C
def indicator(C, z):
    a, b = C[0], C[1]
    I0 = (b - z) * (a - z).conjugate()
    if len(C) == 2:
        return I0.real
    c = C[2]
    I2 = (a - c) * (b - c).conjugate()
    I1 = I0 * I2
    return -I1.imag if I2.imag < 0 else I1.imag


def inside(C, p):
    return indicator(C, p) <= 0

TypeScript

// Операції над комплексними числами (множення додає полярні кути,
// спряження множить кут на -1)
const sub = (a: Point, b: Point): Point => ({ x: a.x - b.x, y: a.y - b.y });
const mul = (a: Point, b: Point): Point => ({
    x: a.x * b.x - a.y * b.y,
    y: a.x * b.y + a.y * b.x,
});
const conj = (a: Point): Point => ({ x: a.x, y: -a.y });

// I < 0, якщо z всередині C; I > 0, якщо z зовні C; I = 0, якщо z на колі C
function indicator(C: Mec, z: Point): number {
    const a = C[0], b = C[1];
    const I0 = mul(sub(b, z), conj(sub(a, z)));
    if (C.length === 2) {
        return I0.x;
    }
    const c = C[2];
    const I2 = mul(sub(a, c), conj(sub(b, c)));
    const I1 = mul(I0, I2);
    return I2.y < 0 ? -I1.y : I1.y;
}

const inside = (C: Mec, p: Point): boolean => indicator(C, p) <= 0;

Go

// I < 0, якщо z всередині C; I > 0, якщо z зовні C; I = 0, якщо z на колі C
func indicator(C mec, z point) float64 {
    a, b := C[0], C[1]
    I0 := (b - z) * cmplx.Conj(a-z)
    if len(C) == 2 {
        return real(I0)
    }
    c := C[2]
    I2 := (a - c) * cmplx.Conj(b-c)
    I1 := I0 * I2
    if imag(I2) < 0 {
        return -imag(I1)
    }
    return imag(I1)
}

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

Тепер ми можемо нарешті переконатися, що все працює, надіславши розв'язок задачі до Library Checker: #308668.

Задачі для практики

• Library Checker - Minimum Enclosing Circle
• BOI 2002 - Aliens