Алгоритм Форчуна на C++ для побудови діаграми Вороного на площині

Вітаю, шановні читачі даної статті! У статті я дам опис імплементації алгоритму Форчуна (англ. Fortune's algorithm), для побудови діаграми Вороного (англ. Voronoi diagram) з використанням нативних збалансованих двійкових дерев пошуку (для унікальних елементів) (англ. BST, binary search tree), передбачених стандартом C++, — асоціативних впорядкованих контейнерів
std::map
і
std::set
.

Діаграма Вороного для 15 точок
Вибачте мене за все стилістичні, граматичні та фактичні помилки. Прошу вас вказати на них в особистих повідомленнях. Все виправлю. Якщо є якісь упущення, неясності, двозначності — також прошу вас повідомити мені.
Репозиторій
https://github.com/tomilov/sweepline/
Поки немає релізу, але гілка
develop
містить у даний момент стабільну версію. Потім з'явиться і
master
, і релізи.
Визначення діаграми Вороного
Для деякого кінцевого набору попарно різних точок на площині (тут і далі
N
— кількість точок) діаграма Вороного представляє з себе розбиття площини на області — осередку (англ. cell) діаграми Вороного. Кожна комірка містить в собі лише одну точку вихідного набору точок, званих сайтами діаграми Вороного. Всі точки осередку знаходяться ближче до відповідного сайту, ніж до будь-якого іншого. Для звичайного визначення відстані на площині — метрики Евкліда (корінь квадратний із суми квадратів різниць координат точок — з теореми Піфагора) — форма комірок у загальному випадку є опуклим гратки. Звичайно ж існують крайні випадки, як якщо у вихідній множині 1 (комірка — вся площина), 2 точки (комірка — полуплоскость), або точки (
N >
) знаходяться на одній прямій (в цьому випадку внутрішніми осередками будуть смуги, а зовнішніми — півплощини). Крайні з безлічі комірок у загальному випадку є частиною опуклих багатокутників з двома сторонами, що йдуть в нескінченність, то є паралельними, або розбіжними променями. Сторони многокутників обмежують осередку (але при цьому не є їх частиною) — це ребра (англ. edge) діаграми Вороного. Вони можуть бути відрізками (англ. line segment), променями (англ. ray), або прямими (англ. line). Кожна ребро — це безліч точок, рівновіддалених рівному від двох сайтів, тобто лежить на середньому перпендикулярі для двох сайтів. Вершини многокутників, що обмежують осередку, так і називаються — вершини діаграми Вороного (англ. мн. ч. vertices, од. ч. vertex). Вершини є точками, рівновіддаленими від трьох і більше сайтів, тобто є центрами (англ. circumcenter) описаних кіл (англ. circumcircle) багатокутників, які можна було б побудувати на сайти примикають комірок, на вершинах.
Неповне опис алгоритму Форчуна
Раніше, в опублікованих тут іншими авторами статтях ((1, 2, 3)), вже було наведено опис алгоритму Форчуна і деяких структур даних, які використовуються для його реалізації. Коротенько, суть алгоритму в наступному: по площині рухається замітаюча пряма (англ. sweepline). Рухається стрибками від точки до точки. Перша частина цих точок — це точки з введення, які стають сайти. Друга частина — це "віртуальні" точки, крайні по ходу руху замітаючої прямої точки згаданих описаних кіл. При русі (паралельному перенесенні) замітаючої прямої вона стосується будь-такий описаної окружності двічі — другий раз еквівалентний події, при якому діаграма Вороного добудовується: до неї додається вершина, одне або більше ребер закінчуються цією новою вершиною і одне або два нових ребра виходять з неї.
Не думаю, що наведене вище опис алгоритму досить добре, тим більше, в голові у мене своя картина того, як це все відбувається. Картина ця сформувалася під впливом мови C++ і частини його стандартної бібліотеки — стандартної бібліотеки шаблонів — STL (а саме, алгоритмів і контейнерів). І тут ми плавно переходимо до наступної частини нашої розповіді.
Алгоритм і структури даних
"Натягнути" алгоритм Форчуна STL — непросте завдання. Точніше, зовсім не очевидно, як саме це зробити. Пов'язано це з хитрою структурою використовується в алгоритмі, яка називається "берегова лінія" (англ. beachline). Вона описує геометричне безліч точок рівновіддалених від замітаючої прямої і сайтів, осередку яких ще не були замкнуті. В нашому випадку (площина з звичної Евклідовою метрикою) берегова лінія складається з дуг (англ. arc) — обмежених з однієї або з двох сторін частин парабол (англ. parabola), фокусами (англ. од. ч. focus, мн. ч. focii) яких є сайти. В ході переміщення замітаючої прямої ці шматки парабол ростуть, діляться (у разі якщо замітаюча пряма перетинає черговий сайт і перпендикуляр до замітаючої прямої з цього місця потрапляє в шматок параболи десь посередині, то цей шматок параболи ділиться на два і між ними "вклинюється" ще один — з фокусом у новому сайт), і стискаються в точку (в цьому випадку з'являється нова вершина і одне або два нових ребра). Шматки парабол формують берегову лінію межують між собою endpoint-ах (укр. кінець відрізка). Кожен даний endpoint рухається по прямій (тобто ребру) і є точкою, яка рівновіддалена в кожен момент: від замітаючої прямої і від кожного з двох сайтів, які є в свою чергу фокусами межують через цей endpoint двійки шматків парабол.
Так от, ця структура, по-перше, повинна бути впорядкованою, по-друге, вона зазвичай є у авторів існуючих імплементацій гетерогенною, тобто містить як шматки парабол, так і endpoint-и в якості ключів/елементів. Для того, щоб імплементація мала середню і максимальну асимптотическую складність
O(N * log(N))
, з необхідністю доведеться використовувати саме збалансовані сортують дерева. Для зберігання гетерогенних даних можна використовувати
std::variant
. Останній, втім, мені не знадобиться.
std::map і неточний користувальницький "прозорий" компаратор
Здавалося б (!) контейнери STL
std::set
та
std::map
не дуже підходять для зберігання берегової лінії, але, як виявляється, їх можна використовувати на повну котушку. Гетерогенності даних у вузлах дерев можна повністю уникнути. Вся справа в тому, що в якості ключів можна зберігати саме endpoint-и — впорядковані пари покажчиків (або ітераторів, що вказують) на сайти — "лівий" і "правий" (у деякому сенсі, який стане ясним далі). У цьому полягає неминуча надмірність: кожна пара сусідніх endpoint-ов містить однакову пару покажчиків (ітераторів, що вказують) на сайт.
Шматки парабол ростуть і стискаються, множаться при русі замітаючої прямої, але моя імплементація стежить за тим, щоб завжди в будь-який момент інваріанти контейнера (унікальність/неэквивалентность і абсолютна упорядкованість елементів в ньому) зберігалися. endpoint-и видаляються до того, як абсолютна впорядкованість контейнера може бути порушена.
Сучасні реалізації стандартних бібліотек C++ (принаймні
libc++
та
libstdc++
) мають імплементацію цікавить мене контейнера
std::map
, яка дозволяє виробляти гетерогенний пошук (англ. heterogeneous lookup), а саме, функції:
count
,
find
,
equal_range
,
lower_bound
upper_bound
, — мають серед перевантажень (англ. overloadings) шаблон функції, який дозволяє проводити пошук елементів, еквівалентних значенням, що має тип відмінний від типу значення (або від типу ключа значення, у разі відображення) зберігається в контейнері. Дана перевантаження кожної з функцій може бути використана тільки якщо компаратором (функтором, використовуваним для порівняння елементів) є "прозорий компаратор", тобто такий компаратор, який у своєму просторі імен визначає тип з іменем
is_transparent
, будь то таг (ім'я класу (
struct
,
class
,
union
) або перерахування (
enum
або
enum class
)),
typedef
або псевдонім типу (англ. type alias), як це наведено в наступний код:
struct less
{
using is_transparent = void; // int, std::true_type, std::false_type...
// або
// typedef void is_transparent;
// або
//struct is_transparent;
//class is_transparent;
//union is_transparent;
// або
//enum is_transparent;
//enum class is_transparent;
bool operator () (const T & lhs, const T & rhs) const
{
return /* повернути істину, якщо lhs строго менше rhs */;
}
};

Як видно, тип імені (англ. qualified-id)
less::is_transparent
може бути незавершеним (англ. incomplete). Наприклад,
void
або оголошення імені може бути forward-declaration класу/перерахування.
Для того щоб компаратор (
less
у прикладі) міг порівнювати значення типу
T
, що знаходяться в контейнері, зі значеннями деякого іншого типу
U
, необхідно довизначити компаратор
less
ще парою операцій, еквівалентних
less::operator () (T, U) const
та
less::operator () (U, T) const
для кожного такого типу
U
.
Еквівалентність з допуском
eps

Інтерфейс (англ. custom) компаратор тут ще необхідний у зв'язку з наступними міркуваннями: у випадку, якщо для проміжних обчисленні використовуються числа з плаваючою комою (англ. floating-point numbers) мають обмежену точність, то алгоритм параметризуется позитивною константою
eps
, який задає точність порівняння. Числа, які відрізняються на величину меншу ніж
eps
вважаються еквівалентними:
double eps = 0.2;
double a = 2.0;
double b = 2.1;
assert(!(a + eps < b) && !(b + eps < a));

assert
-іоні стоїть якраз затвердження (визначення) еквівалентності
a
та
b
. В даному додатку — еквівалентності з точністю
eps
.
Ось приклад виводу списку значень ключів (дійсних (з нек. застереженнями) чисел) всіх елементів, що містяться в контейнері
std::map
або
std::set
та еквівалентних цілому числу
1
(взято звідси):
#include <set>
#include <map>
#include <iterator>
#include < algorithm>
#include < iostream>

const double eps = 0.2;

struct less
{
bool operator () (double l, double r) const { return, l < r; }
using is_transparent = void;
bool operator () (int l, double r) const { return l + eps < r; }
bool operator () (double l, int r) const { return l + eps < r; }
};

int main()
{
{
std::set< double, less > s{0.0, 0.9, 1.0, 1.1, 2.0};
auto p = s.equal_range(1);
std::copy(p.first, p.second, std::ostream_iterator< double >(std::cout, " "));
std::cout << std::endl;
}
{
struct A {}; // irrelevant
std::map< double A, less > m{{0.0, {}}, {0.9, {}}, {1.0, {}}, {1.1, {}}, {2.0, {}}};
auto p = m.equal_range(1);
for (auto it = p.first; it != p.second; ++it) {
std::cout << it->first << ' ';
}
}
}

Має вивести
0.9 1 1.1
для кожного контейнера. Ясно, що навіть в контейнері, що містить унікальні значення якогось типу
T
, може виявитися більш одного елемента, еквівалентного деякого значення деякого типу
U
, відмінного від
T
. При цьому для роботи алгоритмів
std::map::{count,upper_bound,lower_bound,equal_range}
має виконуватися умова — елементи контейнера повинні бути хоча б частково впорядковані (англ. partitioned) при порівнянні зі значенням типу
U
, як це описано у цьому відповіді на SO.
Вставка з підказкою
Корисна особливість реалізації асоціативних впорядкованих контейнерів STL — це можливість вставити елемент з підказкою (англ. hint). Функції вставки
std::map::insert
передається параметром крім вставляється значення ще й ітератор на елемент
hint
, перед яким потрібно вставити елемент. У разі успішної вставки стандартом дається гарантія, що буде вироблено лише обмежене порівняння елементів контейнера (м/у собою (хоча чому?) і зі вставленим). Тут я, зізнаюся, використовую невизначений поведінка (англ. UB, undefined behaviour): справа в тому, що є деякі вимоги до компаратора — Comparator concept — які повинні виконуватися завжди, що б ми не робили з контейнером. Проте ж, логічно можна укласти (і так на практиці і виходить з
libc++
та
libstdc++
), що при вставці унікального елемента в правильне місце за підказкою необхідно максимум два порівняння додаванні елемента з попереднім
*std::prev(hint)
елементом і з елементом, на який вказує ітератор
hint
, якщо такі є. Якщо контейнер порожній, то порівнянь не потрібно, а якщо вставити в початок або в кінець, то тільки лише одне порівняння. Жодна з двох реалізацій стандартної бібліотеки, протестованих мною, ніколи не намагається порівняти ні один з елементів сам з собою або, у випадку зі вставкою з підказкою, — виконувати відмінні від абсолютно логічно необхідних для підтвердження правильності вставки порівняння елементів. Загалом скрізь і всюди торжествує здоровий глузд. Але це все ж UB, якщо ви в своєму коді покладаєтеся на це. Я попередив.
Лексикографічна впорядкованість
Для того, щоб імплементація алгоритму Форчуна володіла внутрішньою красою, необхідно, щоб на вхід алгоритму подавалися точки лексикографічно впорядковані (англ. lexicographically ordered). Лексикографічна впорядкованість означає, що спочатку дві точки порівнюються за абсциссе (x, "ікс"), а потім, у разі якщо їх перші координати еквівалентні, то за ординате (y, "ігрек"). Причому в моєму випадку — з допуском
eps
.
Дві точки площини з координатами
(lx, ly)
та
(rx, ry)
можна порівняти лексикографічно з допуском
eps
:
bool less(double l,
double eps,
double r)
{
return l + eps < r;
}

bool less(double lx, double ly,
double eps,
double rx, double ry)
{
if (less(lx, eps, rx)) {
return true;
} else if (less(rx, eps, lx)) {
return false;
} else {
if (less(ly, eps, ry)) {
return true;
} else {
return false;
}
}
}

Лексикографічно без допуску можна порівняти дві точки, використовуючи готову реалізацію з STL так:
#include <tuple>

bool less(double lx, double ly,
double rx, double ry)
{
return std::tie(lx, ly) < std::tie(rx, ry);
}

Детальніше про це тут. Як
std::tuple
та
std::pair
мають саме такі операції
operator <
.
Введення і виведення, формати даних
Весь передбачається введення лексикографічно впорядковані і не містить еквівалентних з точністю
eps
точок. Так як сортування введення — це теж
O(N * log(N))
, то цю (логічно цілком окрему) частину я виношу за межі імплементації. Це дозволяє використовувати на вході проксі-ітератори (є в прикладі), тобто сортувати ітератори вказують на великовагові точки, а не самі точки.
Алгоритм виглядає як алгоритми STL — на вхід приймає пару ітераторів
[first,last)
, що задають діапазон точок (сайтів), для обробки. Формат точки:
struct point { value_type x, y; };

Де
value_type
— тип значення координати і, одночасно, тип для зберігання будь-яких проміжних значень в ході дій алгоритму.
Результатом є граф: вершини містяться в (автоматично упорядкованому в ході роботи алгоритму) списку
vertices_
, ребра — у двосторонній черзі
edges_
. Контейнери для вершин (
std::list
) і ребер (
std::deque
) обрані такими, щоб не було проблем з недійсними итераторами при вставки та видалення.
Формат ребра:
struct edge { site l, r; pvertex b, e; };

Де
site
— покажчик (або ітератор, який вказує на) сайт,
pvertex
— ітератор, який вказує на вершину. Самі вектора в парі
(l, r)
та
(b, e)
упорядковані за годинниковою стрілкою, тобто так, що
[(l, r) x (b, e)] < 0
(тут — векторний добуток). Останнє дає можливість все правильно намалювати в кінці і, взагалі, зберігати всю повноту інформації про графа.
Вершина
nov = std::end(vertices_)
вказує на нескінченність. Копіювати і переміщати цю пару контейнерів так просто не можна. Так як ні при копіюванні (англ. copy), ні при переміщенні (англ. move) дійсні (англ. valid) ітератори, які не можна розв'язати (англ. not dereferencable), в їх числі і end-ітератор
nov
, не зберігаються. Можливо тільки глибоке копіювання (англ. deep-copying, cloning (є в прикладі)), яке тут має квадратичну складність. Хоча є один обхідний маневр. На самому початку зберегти фейковий елемент
vertices_
і використовувати його ітератор як покажчик на нескінченність. Я не пробував, але це принципово можливо, але тільки для копіювання/переміщення, не для клонування.
Інші особливості імплементації
Зазвичай алгоритм Форчуна веде чергу для подій двох типів. Це так звані події кола (англ. circle event) і події точки (англ. point event).
Події точки виникають, коли замітаюча пряма перетинає черговий сайт. Вони неминучі.
Події кола виникають, коли схлопывается якийсь шматок параболи. В чергу вони додаються, коли три сайту
a
,
b
та
c
в порядку їх перерахування у структурі берегової лінії утворюють двійку векторів, упорядковану за годинниковою стрілкою:
[(a, b) x (b, c)] < 0
(тут — векторний добуток). Але події кола можуть бути видалені, якщо схлопывается якийсь із сусідніх шматків парабол, або один з шматків парабол розчленовується шматком параболи, знову створеним для пересекаемого замітаючої прямої сайту. Тому, замість найчастіше використовуваної обгортки контейнерів
std::priority_queue
, я використовую все той же контейнер
std::map
або
std::set
. Складності потрібних операцій для цих контейнерів ті ж самі, що і для
std::priority_queue
, але вони дозволяють видалити не тільки найбільший/найменший елемент, але і будь-який інший (за ітератора — за
O(1)
часу). Зазвичай події з черги з пріоритетами не видаляють, а позначають як недійсне (англ. disabled). Коли така подія трапляється — його пропускають.
Я використовую тільки для подій кола, так як точки на вхід вже надходять лексикографічному порядку. Події точки мають пріоритет перед подіями кола у разі їх еквівалентності. Чергу подій кола має своїм ключем потенційну вершину, яка буде скопійована в список вершин діаграми Вороного у разі виникнення події кола. Структура описує вершину виглядає так:
struct vertex { struct point { value_type x, y; } c; value_type R; };

Де
(x, y)
— власне сама вершина, а
R
— відстань до сайтів сусідніх комірок.
Для того, щоб використовувати лексикографічну впорядкованість всього і вся в моєму алгоритмі, замітаюча пряма паралельна осі ординат і рухається строго в напрямку осі абсцис. Таким чином всі події кола виникають при координатах
c.x + R
для відповідних вершин.
Особливі випадки
Імплементація спрямована на правильну обробку складних особливих випадків. Їм було приділено підвищену увагу.
Особливі випадки, які потребують особливої уваги — це випадки низької розмірності введення, або локальні симетрії і одночасні події:
  • 0, 1, 2 точки
  • 2 і більш перші точки з введення лежать на одній вертикальній прямій (
    x = const
    ) (при додаванні кожного наступного такого сайту створюється тільки одне ребро і тільки один відповідний endpoint додається берегову лінію).
  • більше, ніж три точки розташовані концентрично (при події сходиться більш ніж два endpoint-а, створюється одне нове ребро і відповідний endpoint додається берегову лінію)
  • при русі замітаючої прямої черговий сайт потрапляє в один endpoint (автоматично виникає і завершується подія кола породжується нова вершина і два нових ребра виходять з неї, два відповідних endpoint-а додаються берегову лінію)
  • при русі замітаючої прямої черговий сайт потрапляє в безліч endpoint-ів, які б на наступному кроці алгоритму зійшлися в новій вершині — в цьому випадку новий сайт еквівалентний вже лежить у черзі подій кола самого найближчого події (подія автоматично завершується, породжуючи нову вершину і два нових ребра, що виходять з неї, два відповідних endpoint-а додаються берегову лінію)
Два останніх особливих випадку якраз виявляються гетерогенним порівнянням сайту endpoint-ами в структурі берегової лінії. За допомогою алгоритму
std::map::equal_range
за
log(size(endpoints_))
часу можна знайти весь діапазон endpoint-ов, що сходяться в вершині, еквівалентної поточного сайту. Якщо діапазон (
[first, last)
) не порожній, то ми потрапили в ребро (
assert(std::next(first) == last);
), або в кілька ребер (
assert(1 < std::distance(first, last));
), а якщо порожній (
assert(first == last)
), то ми просто потрапили в шматок параболи десь посередині. Сам цей шматок розчленовується на дві частини, між якими вставляється новий цілісний шматок параболи (яка в цей момент выроджена в промінь, перпендикулярний замітаючої прямої) фокусом в даному сайт. При цьому створюється тільки одне ребро і два відповідних йому endpoint-а додається до берегову лінію — це найбільш часто зустрічається в загальному випадку варіант. Лише в цьому (і в другому особливому випадку) створюється ребро без хоча б одного вже відомого кінця (тобто асоційованої вершини).
При будь-яких подіях, які змінюють структуру берегової лінії всі залучені сусідні пари endpoint-ів перевіряються на предмет того, чи можуть вони надалі зійтися у вершині. Якщо можуть, то додається відповідна подія в чергу подій кола.
Вставити endpoint-ов
Як я вже згадав вище, мій користувальницький компаратор, який використовується в дереві представляє берегову лінію, є/призводить до UB, при порівнянні endpoint-ов. Для гетерогенної частині ж — все в суворій відповідності зі Стандартом — language lawyer носа не підточить.
Алгоритм влаштований так, що при вставці за підказкою
hint
, нові endpoint-и вставляються справа наліво тільки правильні місця (виявлені на попередніх кроках), і, так як у випадку вставки тільки одного endpoint-а покажчики (ітератори) на хоча б один разграничиваемый сайт присутня в будь-якому з сусідніх endpoint-ах (або
hint
або
std::prev(hint)
), то завжди можна зрозуміти тільки порівнюючи ці покажчики (ітератори), знаходиться порівнюваний endpoint лівіше (нижче) або праворуч (вище) у структурі берегової лінії (прим. подумки я завжди повертаю всю картину на +90°).
Проте ж, є і варіант, коли ми вставляємо два endpoint-а, які виникають в самому останньому особливому випадку. Серед чотирьох endpoint-ов, залучених в операцію вставки за підказкою (вставляються два плюс два граничних), безумовно є такі, які не мають спільних разграничиваемых сайтів. Але і тут є простий вихід. Коли ми вставляємо endpoint-и (як зазвичай: зверху вниз), то кожен з вставляються endpoint-ов має одним з сайтів той, який і викликав це подія точки. Лексикографічно порівнюючи по одному сайту, лексикографическому більшому з двох для кожного endpoint-а, ми можемо з'ясувати, який з endpoint-ов "правіше", а "лівіше" у структурі берегової лінії.
Примітка
також UB в компараторі, імплементація алгоритму покладається на так звані SCARY-ітератори ((4, 5)), які допустимі починаючи з C++11, але не обов'язкові. Вони тут дійсно необхідні для крос-референсинга черзі подій кола і берегової лінії. Інакше (без SCARY-ітераторів) не можна красиво отримати еталонну асимптотическую складність
O(N * log(N))
. Використовувати стирання типу (англ. type erasure) тут для таких цілей виглядало б абсолютно протиприродним.
Пучки
У разі концентрично розташованих сайтів, в списку, який зберігає покажчики endpoint-и, зрослі при відповідному подію, самі ці покажчики знаходяться в деякому безладді. Щоб "дешево" зрозуміти який з сходяться endpoint-ів самий "лівий", а який самий "правий", можна, користуючись тим, що всі відповідні сайти знаходяться на однаковій відстані від загальної вершини, порівнювати кути векторів (або кути серединних перпендикулярів пар сайтів, використовуючи перетворення координат для повороту
+pi / 2 <=> x <- y, y <- -x
), для чого зручна C-функція
std::atan2
, правильно зіставляє знаки і квадранти.
Візуалізація
Для візуалізації я використовую gnuplot. Він цілком швидко малює референтну для мене діаграму Воронного для 100000 сайтів.
Сам алгоритм обрізки ребер за розмірами прямокутника, що обмежує (англ. bounding box) також є в репозиторії.
Опціонально можна включити відображення описаних кіл (виходить така кучерява кольорова картинка) та номерів сайтів. Обидві опції мають сенс тільки на маленьких розмірах введення.
Складний введення
Складними для обробки поганими реалізаціями є наступні конфігурації точок на вході алгоритму:
  • прямокутна сітка (англ. rectangle grid) — точки розташовані через рівні кроки по обом координатам без зміщення
  • діагональна сітка (англ. diagonal grid) — повернена на
    pi / 4
    прямокутна сітка
  • гексагональна сітка (англ. hexagonal grid) — точки розташовані в вершинах правильних (рівносторонній) трикутників, якими вимощено потрібна частина площини
  • трикутна сітка (англ. triangular grid) — точки розташовані в вершинах правильних шестикутників, якими вимощено потрібна частина площини
Міцний горішок (потрібна допомога)
Кожна з перерахованих вище сіток виявила тонни багів на проміжних стадіях моєї імплементації. Остання — це сама жесть, скажу я вам! Вона дуже швидко виявляє неточність в обчисленнях з плаваючою комою. Вже на невеликих разбросах діапазонів абсолютних значень координат різниця в координатах, які обчислюються для центру і радіусу (англ. circumradius) описаної окружності цим формулами (обчислення вершини — ф-я
make_vertex
) та виведеним вручну формулами для обчислення ординати перетину двох горизонтальних парабол (ф-я
intersect
з компаратора для endpoint-ів), сягає зловісно великих значень! Безконтрольне збільшення
eps
— це не вихід, я впевнений.
Не дивлячись на те, що в кожній з цих операцій я домігся використання тільки одного (SIMD :-) ) ділення і одного кореня квадратного, точність їх відрізняється досить сильно. Начебто і
std::sqrt
є точною операцією. Цитата звідси:
std::sqrt is required by the IEEE standard be exact.
Однак істотна розбіжність — встановлений факт. Що робити? Я не знаю. Є тут фахівці з стабільним обчислень з плаваючою комою? Допоможіть — у подяку запишу в співавтори
README
в репозиторії. Мета — перемогти трикутну сітку на нативних типах чисел з плаваючою комою.
Тест продуктивності
Продуктивність я вимірюю для 100'000 точок. Передбачаючи питання "Чому саме для ста тисяч?" я одразу повідомляю: цієї задачі, яка є топ 5 в рейтингу Тимуса за складністю фігурує ця цифра.
Так от, продуктивність для 100'000 точок у загальних позиціях з аллокатором
TCMalloc
Google Tools Performance на досить стародавньому ноутбуці з
Intel® Core(TM) i7-2670QM CPU @ 2.20 GHz
виражається у величині 205 мілісекунд. На серверному залозі без цього хитромудрого аллокатора вона буде приблизно такою ж (перевірено).
до Речі, arena-аллокатор дає таку ж надбавку як і згаданий
TCMalloc
.
Ліцензія
Вона все одно буде нікчемною в Росії. Так що нічого не буду писати тут поки що. Проконсультуйте, якщо є якісь стоять міркування.
Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.