Omówienie zadania Komunikacja międzyplanetarna (III etap XXVIII OI)

Dane jest \(n\) punktów na płaszczyźnie i dla każdego z nich chcemy obliczyć sumę odległości do wszystkich pozostałych punktów. Przy czym w treści wyraźnie zaznaczone jest, że nasze obliczenia nie muszą być bardzo dokładne: mogą różnić się od rzeczywistych aż o \(0{,}1\,\%\), to znaczy wystarczy podać odpowiedzi z błędem względnym nie większym niż \(10^{-3}\).

Podzadanie 1

Nie trzeba dość długo się namyślać, żeby zaproponować rozwiązanie działające w czasie \(O(n^2)\). Wystarczy po prostu obliczyć odległość niezależnie dla każdej pary punktów, a następnie dodać ją do sum dla tych punktów. Co więcej, jeżeli będziemy liczyć te sumy odległości na zmiennych zmiennoprzecinkowych double lub long double, to uzyskana przez nas dokładność będzie o wiele większa niż ta wymagana w zadaniu.

Podzadanie 2

Wszystkie punkty leżą na jednej prostej, zatem zadanie tak naprawdę staje się jednowymiarowe. A w jednym wymiarze dość prosto jest obliczyć szukane sumy odległości, stosując proste zamiatanie.

Na początek przenumerujmy nasze punkty tak, by nowa numeracja odzwierciedlała ich kolejność na prostej. Zauważmy, że jeżeli prosta nie jest pionowa, to wystarczy posortować punkty po współrzędnych \(x\). Ta sztuczka nie będzie oczywiście działać dla prostej pionowej, w której współrzędne \(x\) dla wszystkich punktów są równe, ale wtedy możemy posortować punkty po współrzędnych \(y\). Ale nie będziemy musieli przejmować się tym szczególnym przypadkiem, jeżeli posortujemy punkty leksykograficznie po parach \((x, y)\).

Następnie obliczmy odległości między parami kolejnych punktów: \(d_1, d_2, \ldots, d_{n-1}\). Sumą odległości z punktu \(n\)-tego do wszystkich pozostałych jest \(\sum_{i} d_i \cdot i.\) Jest tak dlatego, że długość odcinka \(d_1\) jest uwzględniona w sumie tylko raz (przy odległości od punktu \(1\) do \(n\)-tego), długość odcinka \(d_2\) jest uwzględniona dwa razy (przy odległości od punktów \(1\) i \(2\) do punktu \(n\)-tego) itd.

Teraz mając już policzoną sumę dla wierzchołka \(n\)-tego, bardzo prosto jest ją uaktualnić, żeby stała się sumą dla wierzchołka \(n-1\). Wystarczy bowiem zmienić współczynnik przy odległości \(d_{n-1}\) z \(n-1\) na \(1\), ponieważ wkłady wszystkich odcinków na lewo od punktu \(n-1\) pozostają takie same, a odcinek \(d_{n-1}\), który wcześniej był sumowany w \(n-1\) odległościach, teraz jest sumowany tylko w jednej odległości pomiędzy punktem \(n-1\) a punktem \(n\). Analogicznie przy uaktualnianiu sumy dla punktu \(n-1\) do sumy dla punktu \(n-2\) musimy zmienić współczynnik przy \(d_{n-2}\) z \(n-2\) na \(2\). Postępując tak samo dla kolejnych punktów, każde uaktualnienie zrobimy w czasie stałym, więc wszystkie sumy obliczymy w czasie \(O(n)\).

Z uwagi na konieczność wykonania początkowego posortowania, cały algorytm dla podzadania 2 będzie działał w czasie \(O(n \log n)\).

Pełne rozwiązanie

Kluczowy pomysł jest następujący. Załóżmy, że mamy dwa punkty \(A=(x_A,y_A)\) i \(B=(x_B,y_B)\), których różnice we współrzędnych to odpowiednio \(x = x_B-x_A\) oraz \(y = y_B-y_A\). Załóżmy również, że nachylenie odcinka łączącego punkty \(A\), \(B\) w stosunku do osi OX jest dosyć małe, wobec tego \(x\) jest dużo większe niż \(y\). Prawdziwa odległość euklidesowa pomiędzy tymi dwoma punktami wyraża się oczywiście wzorem \(d = dist(A,B) = \sqrt{x^2 + y^2}\), ale ponieważ \(x \gg y\), ta odległość powinna być dość bliska wartości samego \(x\). Zastanówmy się zatem, kiedy możemy przybliżyć wartość \(d\) wartością \(x\), aby uzyskać wystarczającą w zadaniu dokładność.

Mamy podać wynik z dokładnością względną \(\epsilon = \frac 1{1000}\), co oznacza, że \[|d-x| \leq d \cdot\epsilon.\] Oczywiście w treści zadania mowa jest o błędzie dla sumy odległości z danego punktu do wszystkich pozostałych, ale jeżeli uda nam się ograniczyć błąd dla każdej odległości z punktu \(A\) do innego punktu, to będzie to implikowało ograniczenie błędu dla sumy tych odległości. Ponieważ \(d > x\), więc możemy opuścić znaki wartości bezwzględnej i wykonać kilka prostych przekształceń algebraicznych: \[\frac x {1-\epsilon} \geq d.\] Następnie podstawiamy definicję odległości i podnosimy obie strony do kwadratu: \[\frac{x^2} {(1-\epsilon)^2} \geq x^2 + y^2.\] Po jeszcze kilku przekształceniach algebraicznych dostajemy ograniczenie na iloraz \(\frac{y}{x}\) w zależności od dokładności \(\epsilon\): \[\frac{y}{x} \leq \sqrt{ 1/{(1-\epsilon)^2} - 1} = 0{,}04475\ldots\] Zatem różnica na współrzędnych \(y\) nie może być większa niż około \(4{,}5\,\%\) różnicy na współrzędnych \(x\).

Ustalmy sobie zatem punkt \(A\) i zobaczmy, gdzie znajdują się takie punkty \(B\), że ich odległość od \(A\) może być aproksymowana rzutem na oś OX. Muszą one znajdować się w kącie zaczepionym w punkcie \(A\), którego dwusieczna jest równoległa do osi OX, a rozmiar tego kąta (który oznaczymy przez \(\alpha\)), jest ograniczony ograniczeniem na iloraz \(\frac{y}{x}\), a konkretnie wyraża się wzorem: \[\alpha \leq 2\arctan\sqrt{ 1/{(1-\epsilon)^2}-1} = 5{,}1251\ldots^{\circ}\]

Zatem odległości wszystkich punktów leżących w kącie \(\alpha\) od jego wierzchołka możemy aproksymować rzutami tych odległości na oś OX. Jeśli zatem oznaczymy te punkty przez \(B_1, B_2, \ldots, B_k\), to ich wkład do sumy dla punktu \(A\) możemy przybliżyć przez \[\sum_{1\leq i\leq k} dist(A,B_i) \approx \sum_{1\leq i\leq k} |x_A - x_{B_i}| = -k\cdot x_A + \sum_{1\leq i\leq k} x_{B_i}.\] Zatem wystarczy nam znać liczbę punktów \(B_i\) znajdujących się w kącie oraz sumę ich współrzędnych \(x_{B_i}\).

Oczywiście takie samo rozumowanie z kątem możemy zastosować dla dowolnego wyboru punktu \(A\), a ponieważ dla uzyskania lepszej złożoności niż \(O(n^2)\) dla całego zadania musimy jakoś uwspólnić pracę wykonywaną przy obliczeniach sum dla różnych punktów, więc zastosujemy tutaj drugi kluczowy pomysł: spróbujemy policzyć szukane wartości jednocześnie dla wszystkich kątów zaczepionych we wszystkich punktach z wejścia. Pomysł jak to zrobić nie jest nowy i pojawił się już wcześniej m.in. w zadaniu Lampy słoneczne z finału XXI Olimpiady Informatycznej.

Na początek przekształcamy naszą płaszczyznę afinicznie tak, aby każdy kąt przeszedł na kąt rozmiaru \(90^\circ\), którego ramiona są równoległe do osi układów współrzędnych, a wnętrze jest górną prawą ćwiartką. Teraz będziemy zamiatali punkty od prawej do lewej pionową prostą. Każdy napotkany punkt będzie wrzucony na prostą ze swoją nową współrzędną \(y\) oraz wagą, będącą starą współrzędną \(x\), czyli tę sprzed przekształcenia płaszczyzny. Jednak przed wrzuceniem punktu na miotłę odczytamy z niej potrzebne informacje do obliczeń dla kąta zaczepionego w tym punkcie. A te informacje to liczba punktów na miotle o współrzędnej \(y\) większej od współrzędnej aktualnego punktu oraz suma wag dla tych wszystkich punktów z miotły.

Tak więc nasza miotła będzie potrzebowała udostępnić dwie operacje:

• dodaj punkt na miotłę o współrzędnej \(y\) i wadze \(w\),

• zwróć liczbę punktów o współrzędnej większych niż \(y\) oraz sumę ich wag.

Obie te operacje możemy zrealizować w czasie \(O(\log n)\), używając standardowego drzewa przedziałowego. Tak więc całkowita złożoność czasowa tego pomysłu to będzie \(O(n \log n)\), ponieważ musimy najpierw posortować punkty, a następnie wykonać zamiatanie, implementując miotłę przy pomocy drzewa przedziałowego. Tak więc w sumarycznym czasie \(O(n \log n)\) obliczymy dla każdego punktu z wejścia wkład do odpowiedzi z wszystkich punktów zawartych w jego kącie \(\alpha\).

No dobrze, ale co z punktami, które leżą poza tymi kątami? Okazuje się, że możemy do nich zastosować ten sam algorytm. Jeśli bowiem obrócimy całą płaszczyznę przeciwnie do ruchu wskazówek zegara o kąt \(\alpha\), to w rozważanych przez nas kątach znajdą się punkty, które wcześniej znajdowały się pod dolnym ramieniem kąta. Możemy zatem takim samym algorytmem policzyć ich wkład do sumy. Natomiast wkłady z kolejnych punktów uzyskamy obracając płaszczyznę o \(2\alpha\), \(3\alpha\) itd.

Ile razy zatem będziemy musieli wywołać jedną fazę zamiatania, żeby dla każdego punktu policzyć całą sumę odległości? Musimy to zrobić tyle razy, żeby kątem \(\alpha\) pokryć cały kąt pełny, co przy żądanej dokładności daje \(360^\circ/\alpha = 70{,}242\ldots\) Możemy zatem nieznacznie zmniejszyć wartość kąta \(\alpha\), tak by w dzieleniu uzyskać liczbę całkowitą, równą co najmniej \(F = 71\).

Zatem całkowita złożoność czasowa całego programu to będzie \(O(F n \log n)\).

Kwestia dokładności obliczeń

Nasze rozwiązanie koncepcyjnie wygląda bardzo fajnie, ale w praktyce czeka na nas jeszcze jedna pułapka związana z dokładnością obliczeń. A mianowicie, żeby to rozwiązanie działało poprawnie, to dla ustalonego punktu \(A\) każdy z innych punktów \(B\) musi wylądować w dokładnie jednym kącie. Musimy być zatem ostrożni z punktami \(B\), które leżą bardzo blisko ramion kątów. Potencjalnie bowiem błędy zaokrągleń mogą spowodować, że zamiast zaliczyć ten punkt do dokładnie jednego z tych kątów, uwzględnimy go przy obu z nich albo żadnym. A to może spowodować niekontrolowaną zmianę wyniku i w konsekwencji błędną odpowiedź.

Okazuje się, że przy zachowaniu umiarkowanej ostrożności jesteśmy w stanie rozwiązać to zadanie tak, jak zostało to zaprezentowane, to znaczy z obracaniem płaszczyzny, używając zmiennych zmiennoprzecinkowych.

Jednak jeżeli chcemy mieć zupełną pewność, że nic złego tu nas nie spotka, to możemy nieco zmodyfikować nasz algorytm, tak żeby przynależność punktów do kątów obliczać przy użyciu tylko zmiennych całkowitoliczbowych. W tym celu cały kąt pełny podzielimy sobie na \(71\) mniej więcej równych części za pomocą wektorów od \(v_1\) do \(v_{71}\) o współrzędnych całkowitych. W tym celu za wektor \(v_i\) możemy wziąć wektor \[(10^8 \cos(i \tfrac{360^\circ} F),\ 10^8 \sin (i \tfrac{360^\circ} F)),\] a następnie zaokrąglić jego współrzędne do liczb całkowitych. Można sprawdzić, że rozmiar kąta utworzonego przez każde dwa kolejne wektory będzie mniejszy niż \(\alpha\), a to już nam wystarczy do rozwiązania.

Teraz, aby sprawdzić przynależność punktu \(B\) do kąta wyznaczonego przez wektory \(v_i\) i \(v_{i+1}\), wystarczy sprawdzić znak odpowiednich iloczynów wektorowych, które będziemy obliczać na zmiennych całkowitych. Co więcej, możemy również zastąpić całą operację afinicznego przekształcania płaszczyzny równoważną z naszego punktu widzenia operacją wykonywaną jedynie na liczbach całkowitych. Zainteresowanych widzów odsyłam do opisu rozwiązania Lampy słoneczne, w którym szczegółowo opisano tę operację.

Złożoność czasowa tej ulepszonej wersji algorytmu jest taka sama jak oryginalnej i wynosi \(O(F n \log n)\).