Spotlight #60 – Geometria obliczeniowa — upakowanie sfer, krzywe wypełniające przestrzeń i triangulacja

I. Upakowanie sfer — hipoteza Keplera i sieć FCC

W 1611 roku Johannes Kepler stwierdził, że najefektywniejszym sposobem upakowania jednakowych sfer jest układ stosowany przez sprzedawców warzyw przy piętrzeniu pomarańczy: heksagonalnie gęsto upakowane warstwy, z których każda jest przesunięta tak, by osiąść we wgłębieniach warstwy poniżej. Daje to gęstość upakowania:

η = π / (3√2)  ≈  0.7405

Oznacza to, że około 74% przestrzeni jest zajęte przez sfery, a 26% pozostaje pustką. Identyczną gęstość osiąga sieć regularna ściennie centrowana (FCC) oraz układ heksagonalny gęsto upakowany (HCP); różnica polega na sposobie ułożenia kolejnych warstw (ABCABC kontra ABABAB). Hipoteza Keplera pozostawała niedowiedziona przez niemal cztery stulecia. Wyzwaniem nie była sama geometryczna intuicja, lecz rygorystyczne wykluczenie każdego możliwego nieregularnego układu, który mógłby potajemnie pakować gęściej.

Dowód Halesa-Fergusona i projekt Flyspeck

Thomas Hales i Samuel Ferguson ogłosili dowód w 1998 roku, sprowadzając hipotezę do skończonego, lecz ogromnego zagadnienia programowania liniowego, a następnie rozwiązując je obliczeniowo. Złożony manuskrypt, wraz z gigabajtami kodu komputerowego, został przyjęty w 2005 roku przez czasopismo Annals of Mathematics — jednak recenzenci słynnie stwierdzili, że są „pewni w 99%" co do poprawności, nie mogąc zweryfikować ręcznie każdego kroku obliczeniowego. To zapoczątkowało projekt Flyspeck: trwające dekadę przedsięwzięcie mające na celu stworzenie w pełni zweryfikowanego komputerowo dowodu formalnego, ukończone w 2014 roku przy użyciu asystentów dowodzenia HOL Light i Isabelle. Flyspeck stał się kamieniem milowym matematyki formalnej, pokazując, że poważny otwarty problem można rozstrzygnąć z absolutną logiczną pewnością za pomocą maszyny.

Upakowanie losowe — czyli zwyczajne wsypanie sfer do pojemnika — osiąga maksymalną gęstość rzędu 0,64 (granica „losowego gęstego upakowania"). Różnica między 0,64 a 0,7405 ma ogromne praktyczne konsekwencje: decyduje o porowatości tabletek farmaceutycznych, przepuszczalności gleby oraz przewodności materiałów kompozytowych. Symulacja pozwala bezpośrednio porównać te reżimy.

II. Krzywa Hilberta — linia, która wypełnia kwadrat

W 1890 roku Giuseppe Peano wstrząsnął światem matematyki, konstruując ciągłą krzywą przechodzącą przez każdy punkt kwadratu jednostkowego — krzywą o dwuwymiarowym obrazie wywodzącym się z jednowymiarowego parametru. Rok później David Hilbert podał czystszą, geometrycznie przejrzystą konstrukcję. Każda iteracja krzywej Hilberta dzieli kwadrat na cztery ćwiartki, łączy je ścieżką w kształcie litery U i rekurencyjnie stosuje tę samą regułę wewnątrz każdej ćwiartki, obracając i odbijając ją w razie potrzeby.

Krzywa graniczna ma wymiar Hausdorffa dokładnie równy 2, dlatego wypełnia płaszczyznę: nie jest krzywą w intuicyjnym sensie, lecz obiektem fraktalnym, który zaciera granicę między linią a powierzchnią. Formalne stwierdzenie mówi, że istnieje ciągła surjekcja z [0,1] na [0,1]², obalając intuicję, że krzywe i kwadraty muszą mieć różne „rozmiary" nieskończoności punktów.

Lokalność danych i pamięć podręczna

Najcenniejszą właściwością krzywej Hilberta jest jej lokalność: punkty bliskie sobie wzdłuż jednowymiarowego parametru zwykle znajdują się blisko siebie także w przestrzeni dwuwymiarowej. Odwrotność nie jest w pełni prawdziwa (niektóre sąsiadujące punkty 2D leżą daleko od siebie na krzywej), ale korelacja jest na tyle silna, że można ją wykorzystać w informatyce. Systemy informacji geograficznej (GIS) używają indeksów opartych na krzywej Hilberta do sortowania danych przestrzennych, tak by sąsiadujące obiekty trafiały do tego samego bloku dysku lub linii pamięci podręcznej. Gdy zapytanie pobiera wszystkie punkty w obrębie prostokąta ograniczającego, dane uporządkowane przestrzennie powodują znacznie mniej chybień pamięci podręcznej niż naiwne uporządkowanie wierszowe.

Silniki baz danych (w tym Bigtable od Google i jego otwartoźródłowi następcy) wykorzystują indeksy Hilberta do grupowania powiązanych wierszy na dysku. Algorytmy kompresji obrazów korzystają z odwzorowań 2D-na-1D (spokrewnionych z krzywą Hilberta) do budowy słowników kwantyzacji. Nawet kontrolery pamięci DRAM eksperymentowały z porządkiem Hilberta, by usprawnić wielowątkowe wzorce dostępu przestrzennego. Idea, że topologia krzywej może optymalizować wydajność sprzętu, jest pięknym przykładem tego, jak abstrakcyjna matematyka znajduje nieoczekiwane zastosowania inżynierskie.

III. Spirala Ulama — ukryta geometria liczb pierwszych

W 1963 roku Stanisław Ulam siedział na nużącym zebraniu i zaczął bazgrać: zapisał 1 w środku, następnie rozwinął liczby całkowite na zewnątrz w kwadratową spiralę i zakreślił liczby pierwsze. Powstały diagram, rozszerzony do tysięcy liczb całkowitych, ukazuje wyraźne przekątne pasma — linie, wzdłuż których liczby pierwsze skupiają się znacznie gęściej, niż wynikałoby to z losowości. Bazgroły Ulama stały się małą sensacją, gdy opublikował je Scientific American.

Wyjaśnienie tkwi w wielomianach kwadratowych. Każda przekątna spirali odpowiada wielomianowi postaci 4n² + bn + c dla ustalonych b i c. Niektóre wielomiany generują liczby pierwsze z wyjątkowo wysoką częstotliwością. Słynny wzór Eulera n² + n + 41 daje liczby pierwsze dla wszystkich n od 0 do 39 — czterdzieści kolejnych liczb pierwszych. W spirali Ulama wielomian ten tworzy jasne przekątne pasmo. Hipotezy Hardy'ego-Littlewooda o liczbach pierwszych (hipoteza B, 1923) dają heurystyczne wyjaśnienie: gęstość liczb pierwszych generowanych przez wielomian jest powiązana z tym, jak często unika on podzielności przez małe liczby pierwsze, co ujmuje ilościowo iloczyn czynników lokalnych zwany szeregiem osobliwym.

Dlaczego nie ma dowodu

Mimo wyraźnego wzoru wizualnego nie udowodniono, by jakikolwiek wielomian kwadratowy generował nieskończenie wiele liczb pierwszych. Twierdzenie Dirichleta gwarantuje nieskończenie wiele liczb pierwszych w postępach liniowych an + b (gdy NWD(a,b) = 1), lecz analogiczny wynik dla wielomianów kwadratowych — czwarty problem Landaua — pozostaje otwarty. Spirala Ulama wizualizuje więc jedno z najgłębszych otwartych pytań teorii liczb: nieregularny rozkład liczb pierwszych jest jednocześnie lokalnie uporządkowany i globalnie tajemniczy.

IV. Triangulacja Delaunaya — najlepsza siatka trójkątów

Mając zbiór punktów na płaszczyźnie, można je połączyć w nienachodzące na siebie trójkąty na wiele różnych sposobów. Triangulacja Delaunaya to jedyna taka triangulacja, która maksymalizuje minimalny kąt wewnętrzny spośród wszystkich trójkątów — równoważnie, unika ona „wiórowatych" trójkątów w możliwie największym stopniu geometrycznie. Kryterium definiujące to własność pustego okręgu opisanego: dla każdego trójkąta w triangulacji jego okrąg opisany nie zawiera wewnątrz żadnego innego punktu zbioru.

Triangulacja Delaunaya jest dualem geometrycznym diagramu Woronoja: łącząc środki sąsiadujących komórek Woronoja (komórek mających wspólną krawędź), otrzymuje się triangulację Delaunaya. Ta dualność oznacza, że obliczenie jednej struktury daje drugą, i tłumaczy, dlaczego triangulacja Delaunaya dzieli przestrzeń w sposób najbardziej „naturalny" względem danego rozkładu punktów.

Algorytm Bowyera-Watsona i zastosowania

Algorytm Bowyera-Watsona buduje triangulację przyrostowo. Zaczyna się od supertrójkąta zawierającego wszystkie punkty. Przy wstawianiu każdego nowego punktu p: znajduje się wszystkie trójkąty, których okrąg opisany zawiera p („złe" trójkąty); usuwa się je, tworząc wielokątną dziurę; ponownie triangułuje się dziurę, łącząc wszystkie jej krawędzie brzegowe z punktem p. Algorytm działa w oczekiwanym czasie O(n log n) i jest prosty do zaimplementowania, co tłumaczy jego powszechność w oprogramowaniu do generowania siatek.

Zastosowania obejmują imponujący zakres dziedzin. W metodzie elementów skończonych jakość elementu (a więc i dokładność numeryczna) jest bezpośrednio powiązana z minimalnym kątem, co czyni triangulację Delaunaya domyślną siatką startową. Oprogramowanie do modelowania terenu triangułuje nieregularne chmury punktów pomiarowych w modele TIN (triangulowaną nieregularną sieć). Metody interpolacji geoprzestrzennej, takie jak interpolacja naturalnego sąsiedztwa, wykorzystują strukturę Delaunaya do wyznaczania wag. Nawet obliczeniowa mechanika płynów używa siatek Delaunaya jako punktu wyjścia, zanim dalsze algorytmy udoskonalające (algorytm Rupperta, algorytmy Chew) narzucą ograniczenia dotyczące kąta i gęstości.

V. Maszerujące kwadraty — wyznaczanie konturów z pól skalarnych

Algorytm maszerujących kwadratów wyznacza izolinie — linie, wzdłuż których pole skalarne przyjmuje określony próg — na regularnej siatce. Dla każdej komórki siatki każdy z czterech narożników znajduje się powyżej lub poniżej izowartości. Istnieje 2⁴ = 16 możliwych konfiguracji narożników, a każda konfiguracja odpowiada konkretnemu zbiorowi odcinków przecinających krawędzie komórki. Przechowując te odwzorowania w tabeli przeglądowej, algorytm „maszeruje" przez każdą komórkę i składa odpowiednie odcinki w linie konturowe.

Dwa z 16 przypadków (przypadki 5 i 10) są niejednoznaczne: obie pary przekątnych narożników mają ten sam znak, więc istnieją dwie równie poprawne topologie konturu. Standardowe rozwiązanie polega na próbkowaniu wartości skalarnej w środku komórki: jeśli środek leży powyżej izowartości, kontur tworzy dwie oddzielne krzywe (przypadek „rozdzielony"); jeśli poniżej — kontur tworzy przecięcie w punkcie siodłowym. Błędne rozstrzygnięcie niejednoznacznych przypadków prowadzi do błędów topologicznych — konturów, które się samoprzecinają lub pozostawiają luki.

Od 2D do 3D: maszerujące sześciany

Uogólnienie trójwymiarowe, maszerujące sześciany (Lorensen i Cline, 1987), działa na sześciennych wokselach o ośmiu narożnikach, dając 2⁸ = 256 konfiguracji zredukowanych dzięki symetrii do 15 unikalnych przypadków. Jest to koń pociągowy algorytmiczny obrazowania medycznego: skanery TK i MRI produkują trójwymiarowe pola skalarne (jednostki Hounsfielda dla gęstości tkanki), a maszerujące sześciany wyznaczają powierzchnie narządów przy określonych progach gęstości. Powstałe siatki trójkątów można wydrukować w 3D jako modele anatomiczne, zaimportować do oprogramowania planowania zabiegów chirurgicznych lub wykorzystać do planowania radioterapii. Praktycznie każda izopowierzchnia, jaką kiedykolwiek widzieliście w wizualizacji naukowej, została wygenerowana przez ten algorytm lub jednego z jego następców.