ebvalaim 2015-06-01 0 Comments

Struktura przestrzeni

Tworzenie wszechświata trzeba od czegoś zacząć. Dobrym początkiem będzie określenie jego kształtu.

Najwygodniejszym i chyba najbardziej oczywistym kształtem jest sześcian. Wówczas każda z trzech współrzędnych będzie liczbą z tego samego zakresu, a dodatkowo bardzo łatwo jest sprawić, aby wszechświat nie miał granic - wystarczy dodać warunek, że opuszczenie sześcianu z jednej strony jest równoznaczne z wejściem do sześcianu z drugiej strony. Przypomina to sytuację znaną z gry w węża, w której wąż opuszczający ekran z prawej strony wracał znowu z lewej - tylko w trzech wymiarach.

Mamy zatem sześcienną przestrzeń, w której każdy punkt może być opisany trzema liczbami: $x,\; y,\; z\; \in (a,b)$. Pojawia się pytanie: jakiego typu danych użyć do reprezentowania współrzędnych punktu? Nie uda się na nie odpowiedzieć bez wcześniejszego sprawdzenia, z jak dużymi liczbami mamy do czynienia.

Chcielibyśmy mieć wszechświat o realistycznej wielkości - to daje bok sześcianu rzędu $10^{11}$ lat świetlnych. Rok to ok. 30 000 000 ($3 \times 10^7$) sekund, sekunda świetlna to ok. 300 000 000 ($3 \times 10^8$) metrów, co daje nam bok sześcianu rzędu $10^{27}$ metrów. Dobrze by było, żeby elementarne kawałki wszechświata były nieco mniejsze niż metr (powiedzmy, milimetr) - co oznacza, że musimy radzić sobie z liczbami mającymi 30 cyfr znaczących. Dużo.

Do obsługi tak dużych liczb często wykorzystuje się liczby zmiennoprzecinkowe, ale one nie będą tutaj odpowiednie. Ponieważ takie liczby przechowują stałą liczbę cyfr znaczących i pozycję przecinka, ich precyzja zależy od wielkości. Np. 64-bitowa liczba zmiennoprzecinkowa może przechowywać 15 cyfr znaczących - co oznacza że dla liczb rzędu miliarda ($10^9$) mają dokładność rzędu jednej milionowej ($10^{-6}$), ale przy liczbach wielkości $10^{15}$ mamy wynik już tylko z dokładnością do jedności.

Jasne jest, że konieczne będzie wykorzystanie jakiejś biblioteki liczb o wielokrotnej precyzji. Zdecydowałem się użyć GMP i opisywać współrzędne liczbami całkowitymi (które, stawiając odpowiednio przecinek, można interpretować jako stałoprzecinkowe liczby rzeczywiste). Czemu nie zmiennoprzecinkowymi? Nawet stosując precyzję odpowiednią do największych wartości współrzędnych, otrzymałbym różną precyzję w różnych punktach przestrzeni. Wolę zachować jednorodną precyzję i nie ryzykować różnic między różnymi miejscami, nawet za cenę nieco mniejszej precyzji w pobliżu 0. Oczywiście w tym przypadku współrzędne nie będą opisywały liczby metrów od początku układu współrzędnych, a raczej liczbę jakichś ustalonych jednostek. Przyjąłem, że metr będzie równy 65536 jednostkom, dzięki czemu będę mógł traktować ostatnie 16 bitów liczby jako miejsca po przecinku.

Generowanie galaktyk

Mając określoną przestrzeń, możemy zacząć się zastanawiać, gdzie w tej przestrzeni znajdują się galaktyki.

Moduł generowania galaktyk będzie musiał być w stanie odpowiedzieć na dwa główne pytania:

• Jakie są położenia galaktyk w zadanym obszarze przestrzeni?
• Czy podane współrzędne galaktyki opisują prawidłową galaktykę?

Drugie pytanie będzie potrzebne raczej do upewniania się, czy zewnętrzne dane są poprawne, ale też się przyda.

Pierwsze pytanie wydaje się proste - powinno wystarczyć wziąć zadany zakres współrzędnych, wygenerować trochę liczb pseudolosowych z tego zakresu wychodząc od zdeterminowanego ziarna i już. Co jednak, kiedy weźmiemy dwa obszary, które na siebie zachodzą? Jak upewnić się, że w części wspólnej otrzymamy te same galaktyki? Z pomocą przychodzą drzewa ósemkowe. Algorytm będzie wyglądał tak:

1. Weź sześcian odpowiadający całej przestrzeni i wygeneruj liczbę galaktyk znajdujących się w nim (jak - to osobny szczegół, może to być np. liczba pseudolosowa wygenerowana z ziarna dla wszechświata)
2. Wygeneruj ziarno dla obecnego sześcianu (np. obliczając hash jego współrzędnych)
3. Podziel sześcian na 8 części (dzieląc każdy bok na 2)
4. Każdej z części przypisz prawdopodobieństwo proporcjonalne do jakiejś funkcji gęstości galaktyk (dla rozkładu jednorodnego - po prostu 1/8)
5. Dla każdego sześcianu wygeneruj liczbę galaktyk, które się w nim znajdą:
   • Oblicz prawdopodobieństwo sukcesu p, dzieląc prawdopodobieństwo dla danego bloku przez sumę prawdopodobieństw dla niego i pozostałych bloków
   • Wygeneruj liczbę galaktyk zgodnie z rozkładem dwumianowym - N = liczba pozostałych galaktyk, p = wygenerowane przed chwilą prawdopodobieństwo
   • Odlicz otrzymaną liczbę galaktyk od całkowitej sumy i powtórz procedurę z pozostałymi blokami
6. Jeśli któryś z mniejszych bloków zawiera tylko 1 galaktykę i ma bok 1, zwróć go jako położenie galaktyki.
7. Powtórz kroki od 2 dla mniejszych sześcianów, które mają niezerową liczbę galaktyk i niepustą część wspólną z zadanym obszarem

Ponieważ w tym algorytmie często dzielimy bloki na 2, wygodnie będzie, jeśli długość boku całego wszechświata będzie potęgą dwójki. Wielkość wszechświata będzie więc opisywana przez liczbę $n$, oznaczającą bok $2^{n-16}$ metrów (przypominam: 1 metr to 65536 = $2^{16}$ jednostek).

Ponieważ zawsze wychodzimy od całego wszechświata i schodzimy do coraz mniejszych części, zawsze w tych samych blokach otrzymamy te same galaktyki, niezależnie od tego, jaki był zadany obszar, którego są one częścią. Ważne jest jedynie, aby zawartość danego bloku była zdeterminowana wyłącznie przez własności jego samego - ale to można osiągnąć, wykorzystując jako ziarno jakąś wartość z nim związaną (np. jakiś hash współrzędnych). Może się wydawać, że to wolny algorytm, bo wymaga zawsze dzielenia wszechświata na 2 aż do najniższego poziomu - ale nawet dla realistycznych wszechświatów daje to ok. 100 podziałów do zejścia do bloków o boku 1. Jedyny problem może się pojawić przy próbie generowania galaktyk w zbyt dużym obszarze, gdy ich liczba w tym obszarze okaże się za duża - ale to jest łatwe do uniknięcia, wystarczy nie zadawać zbyt dużych obszarów.

Warto też zauważyć, że ten sam algorytm będzie się dał później zastosować do innych obiektów, np. gwiazd - będzie trzeba tylko zadać inną funkcję gęstości. To jest jednak sprawa na przyszłość.

Implementacja

Póki co zaimplementowałem opisany wyżej algorytm w C. Od tamtej pory dowiedziałem się jednak o istnieniu ciekawego języka, a mianowicie Rusta. Ponieważ nie pozwala on na wycieki pamięci, a przy tym generuje kod maszynowy o wydajności porównywalnej z C/C++, postanowiłem przepisać powyższy algorytm w tym języku. Obecnie pracuję nad stworzeniem biblioteki, pozwalającej na korzystanie z MPFR z poziomu Rusta (jest potrzebna do rozkładu dwumianowego).

Następny post pojawi się, gdy przepiszę generowanie galaktyk w Ruście :)