15509. Wycieczka 2 [AL_09_06]

Zadanie:
https://pl.spoj.com/problems/AL_09_06

Skrócony opis problemu:

Otrzymujesz na wejściu macierz sąsiedztwa $M$ o rozmiarze $n$ oraz macierz $N$, w której w $i$-tym wierszu i $j$-tej kolumnie powinna być ilość dróg z wierzchołka $i$ do wierzchołka $j$ o długości 2 (a więc z dokładnie jednym pośrednikiem; 2 oznacza ilość krawędzi, a nie wierzchołków). Twoim zadaniem jest zweryfikowanie czy faktycznie dla każdej pary wierzchołków $(i;j)$ na wejściu podano prawidłową liczbę $N_{i,j}$ dróg o długości 2 między tymi wierzchołkami. Jeśli w macierzy $N$ wszystkie liczby się zgadzają wypisz TAK. W przeciwnym wypadku wypisz NIE.

Np. dla macierzy sąsiedztwa:

0 1

1 0

Od wierzchołka 1 można przejść do wierzchołka 1 na 1 sposób (przez wierzchołek 2), od 2 do 2 też na 1 (przez wierzchołek 1), ale już od 1 do 2 i od 2 do 1 nie można (więc na 0 sposobów). Macierz $N$ powinna zatem wyglądać tak:

1 0

0 1

Rozwiązanie naiwne:

Dla każdej pary wierzchołków $a$, $b$ można zliczyć ile jest wierzchołków, do których jest krawędź z $a$ i jednocześnie, z których jest krawędź do $b$.

Złożoność: $O(n^3)$, bo par jest $n^2$ i dla każdej pary sprawdzamy wszystkie $n$ wierzchołków.

Rozwiązanie sprytniejsze:

Autorem tego algorytmu jest Maciej Boniecki.

Algorytm ten jest rozwinięciem poprzedniego z małym przyspieszeniem. Mianowicie macierzy sąsiedztwa nie będziemy przechowywać w macierzy booli czy intów, ale w macierzy unsigned long longów. Konkretnie, podzielimy każdy wiersz i każdą kolumnę na zestawy po 64 boole, które możemy zapisać w jednej liczbie typu unsigned long long, jako że jest ona 64-bitowa. Stworzymy zatem 2 macierze. Jedną o rozmiarze $n$ na $\frac{n}{64}$, a drugą o rozmiarze $\frac{n}{64}$ na $n$. W pierwszej podzielimy wiersze na części (stąd ilość kolumn jest mniejsza), a w drugiej kolumny (stąd ilość wierszy będzie mniejsza).

Np. dla macierzy:

$$\begin{bmatrix} 0&1&1&0\\1&0&0&0\\1&1&0&1\\1&0&1&0 \end{bmatrix}$$

Stworzymy takie macierze:

$$\begin{bmatrix} 6\\8\\13\\10 \end{bmatrix}, \begin{bmatrix} 7&10&9&2 \end{bmatrix}$$

Owe dziwne liczby wzięły się z zamiany 0 i 1 z wierszy i kolumn na system dziesiętny. Np. 3-ci wiersz to $1101_2 = 13_{10}$. Mają one tylko po 1 wierszu/kolumnie, bo 4 bity na luzie zmieszczą się w jednej 64-bitowej liczbie. Dopiero od $n=65$ mielibyśmy 2 wiersze/kolumny.

Dla każdej komórki macierzy $N_{i,j}$ musimy zliczyć na ilu pozycjach w wierszu $i$-tym i kolumnie $j$-tej występuje 1nka. Każda taka 1nka na tej samej pozycji (powiedzmy na pozycji $k$) oznacza, że istnieje droga z $i$ do $k$ oraz z $k$ do $j$, a o to nam właśnie chodzi. Jeśli więc na pozycji $k$ jest 1nka w obu miejscach, to znaczy, że wierzchołek o numerze $k$ jest jednym z pośredników między $i$ a $j$ (czyli mamy ścieżkę $i \rightarrow k \rightarrow j$). Musimy więc zliczyć tych pośredników dla każdego $k$. No i przy tym zliczeniu właśnie przyda się to nasze segmentowanie tych 0 i 1 w 64-bitowe liczby. Możemy bowiem wykonać iloczyn logiczny na dwóch dużych liczbach. Działa on w czasie stałym. Dla każdego zestawu 64 zer i jedynek możemy zatem w czasie stałym "połączyć" wiersz z kolumną. Pozostaje nam jeszcze zliczyć w tym połączeniu 1nki. Robimy to też w czasie stałym za pomocą funkcji __builtin_popcountll. Podsumowując, dla każdej komórki $N_{i,j}$ wykonujemy iloczyn logiczny na każdym z $\frac{n}{64}$ segmentów wierszy i kolumn, a następnie na tym iloczynie __builtin_popcountll i porównujemy czy mają takie same wartości. Jeśli dla którejś komórki się różnią, to wypisujemy NIE. W przeciwnym wypadku TAK.

Złożoność: $O(\frac{n^3}{64})$, przy założeniu, że __builtin_popcountll działa w czasie stałym ($O(1)$). Może dziwić, że mamy stałą w złożoności. Umieściłem ją tam, bowiem jest to duża stała. Maksymalnie $n=2000$. 64 jest więc większe od pierwiastka z $n$, który jest mniejszy od 45. Dla tego zadania, gdzie $n=100$ daje nam to prawie złożoność kwadratową.

Algorytm:

int main()
{
    int n, i, j, k, x, y, ok = 1;
    unsigned long long int rows[n][n], columns[n][n];

    get(n);
    for(i = 0; i < n; ++i)
        for(j = 0; j < n; ++j)
            get(x),
            rows[i][j/64] = (rows[i][j/64]*2)|x, // przesuwamy bitowo starą wartość wiersza w lewo i doklejamy nowy bit (operacje na bitach są wydajniejsze, zamiast | można użyć +)
            columns[i/64][j] = (columns[i/64][j]*2)|x;
    for(i = 0; i < n; ++i)
        for(j = 0; j < n; ++j)
        {
            get(x);
            if(ok) // jeśli już kiedyś któraś komórka się nie zgadzała, to po co sprawdzać następne? nie możemy jednak zrobić breaka, bo musimy wczytać pozostałe liczby
                for(y = k = 0; k <= (n-1)/64; ++k)
                    y += __builtin_popcountll(rows[i][k]&columns[k][j]); // & to iloczyn logiczny
            if(y != x)
                ok = 0;
        }
    puts(ok ? "TAK" : "NIE");
}

Rozwiązanie sprytniejsze nr 2:

Rozwiązanie to jest bardzo podobne do poprzedniego i również jest optymalizacją pierwszego algorytmu.

Różnica polega na zastosowaniu innego kontenera. Mianowicie, skorzystamy z bitseta. Tak jak uprzednio, zrobimy 2 osobne bitsety (macierze składające się z samych 0 i 1) - jeden na wiersze, drugi na kolumny. Nie będziemy jednak dzielić każdego wiersza i kolumny na segmenty, jako że to one już będą segmentami. Bitset zapewnia bowiem szybkie operacje na ciągach bitów - tak samo jak przed chwilą wykonywaliśmy szybkie operacje na samych liczbach. Dla każdej pary: wiersz-kolumna wykonujemy zatem tylko jeden iloczyn logiczny, a następnie na wyniku wykonujemy metodę count, która zlicza nam 1nki w powstałym ciągu bitów (odpowiednik __builtin_popcount).

Złożoność: $O(n^2 \cdot b)$, gdzie $b$ to max ze złożoności operacji iloczynu logicznego na 2 wierszach bitsetu oraz funkcji count. Zakładam, że jest ona szybsza od $O(n)$.

Algorytm:

int main()
{
    int n, i, j, k, x, y, ok;
    bitset<n> rows[n], columns[n];

    get(n);
    ok = 1;
    for(i = 0; i < n; ++i)
        for(j = 0; j < n; ++j)
            get(x),
            columns[j][i] = rows[i][j] = x;
    for(i = 0; i < n; ++i)
        for(j = 0; j < n; ++j)
        {
            get(x);
            if(ok)
                y = (int)(rows[i]&columns[j]).count();
            if(y != x)
                ok = 0;
        }
    puts(ok ? "TAK" : "NIE");
}

Rozwiązanie jeszcze sprytniejsze:

Algorytm ów zdradził mi Damian Straszak.

Aby rozwiązać to szybciej, potrzebujemy trochę teorii.

$n$-ta potęga macierzy sąsiedztwa pokazuje liczbę marszrut o długości $n$ między każdymi dwoma wierzchołkami; np. $M^7[6][9]$ będzie zawierać liczbę marszrut o długości 7 z 6 do 9. Uzasadnienie w akapicie "Ścieżki w grafie".

Marszruta (ścieżka) - ciąg krawędzi, taki że dwie kolejne krawędzie są sąsiednie.
Droga (ścieżka prosta) - marszruta, w której każdy wierzchołek występuje najwyżej raz.
Łańcuch - marszruta, w której każda krawędź występuje najwyżej raz.
Cykl - droga, w której wierzchołek początkowy jest równy końcowemu.

Musimy de facto obliczyć kwadrat macierzy $M$ i porównać z $N$. Jeszcze innymi słowami musimy sprawdzić czy macierz $M^2-N$ jest macierzą zerową.

Złożoność: $O(n^{2.3727})$ ze względu na mnożenie macierzy. Macierze można bowiem wymnożyć naiwnie w czasie $O(n^3)$, za pomocą algorytmu Strassena w czasie $O(n^log_2{7}) \approx O(n^{2.807})$ lub najszybciej właśnie w czasie $O(n^{2.3727})$ dzięki Virginii Williams, która usprawniła algorytm Coppersmitha-Winograda.

Rozwiązanie najsprytniejsze:

Algorytm ów zdradził mi Damian Straszak.

Ten algorytm jest rozwinięciem poprzedniego. Polega na modyfikacji mnożenia macierzy tak, by działała w czasie kwadratowym.

Czy to znaczy, że Damian jako pierwszy na świecie dokonał rewolucyjnego odkrycia i znalazł kwadratowy algorytm mnożenia macierzy? Niestety nie, ale jest myk, dzięki któremu możemy z dużym prawdopodobieństwem sprawdzić czy macierz jest zerowa.  Mianowicie, jeśli przemnożymy macierz zerową przez wektor to otrzymamy wektor zerowy. Jeśli natomiast macierz nie jest zerowa i przemnożymy macierz przez wektor $v$ to prawdopodobieństwo, że otrzymamy wektor zerowy jest mniejsze od $\frac{1}{n}$ (Damian udowodnił to w komentarzu). Co więcej, powyższe twierdzenie działa dla wektorów o dowolnych składnikach (również ujemnych). W naszym wypadku macierz ma elementy nieujemne, więc jeśli przemnożymy macierz $M^2-N$ przez losowy wektor $v$ o współrzędnych z przedziału $\left<1;n\right>$ i otrzymamy wektor zerowy to wypisujemy TAK, a jeśli niezerowy, to NIE. Nie ma liczb ujemnych, więc nic nie zredukuje się do 0 bez przemnożenia przez 0.

Jak jednak wymnożyć $M^2-N$ przez wektor $v$ w czasie kwadratowym? Możemy na przykład osobno obliczyć $R = N \cdot v$ w czasie kwadratowym (bo mnożymy macierz przez wektor, a to ma złożoność $O(n^2)$), osobno $O = M \cdot v$ w czasie kwadratowym i osobno $P = O \cdot v$, a na końcu $P-R$ w czasie liniowym.

Możemy zrobić to też inaczej (korzystam z tego w poniższym programie). Uogólnijmy mnożenie macierzy $M \cdot M \cdot v$ na mnożenie dwóch różnych macierzy $X \cdot Y \cdot$ (o rozmiarze 3).

$$X = \begin{bmatrix} A&B&C\\D&E&F\\G&H&I \end{bmatrix}, Y = \begin{bmatrix} J&K&L\\M&N&O\\P&Q&R \end{bmatrix}, v = \begin{bmatrix} S\\T\\U \end{bmatrix}$$

Po wymnożeniu okazuje się, że wektor wynikowy wygląda następująco:

$$\begin{bmatrix} A(JS+KT+LU)+B(MS+NT+OU)+C(PS+QT+RU)\\D(JS+KT+LU)+E(MS+NT+OU)+F(PS+QT+RU)\\G(JS+KT+LU)+H(MS+NT+OU)+I(PS+QT+RU) \end{bmatrix}$$

Jak widać nawiasy się powtarzają. Możemy je więc obliczyć jednorazowo w czasie $O(n^2)$ - bo jest $n$ nawiasów, a w każdym z nich $n$ elementów (zapiszmy je w wektorze tymczasowym $tmp$). Na koniec dla każdego z $n$ elementów wektora końcowego sumujemy iloczyny elementu z wektora tymczasowego i elementu z wektora $v$ (których rozmiary to $n$).

Złożoność: $O(n^2)$.

Ostateczny algorytm:

int main()
{
    int n, i, j, ok, x;
    int M[n][n], N[n][n], v[n], tmp[n];

    get(n, M, N);
    for(i = 0; i < n; ++i)
        v[i] = i+1, // tutaj wybrałem takie wartości wektora; też dają AC
        tmp[i] = 0; // zerujemy naszą tymczasową tablicę
    for(i = 0; i < n; ++i)
        for(j = 0; j < n; ++j)
            tmp[i] += v[j]*tab[i][j];
    for(i = 0; i < n; ++i)
    {
        for(j = 0; j < n; ++j)
            x += tab[i][j]*tmp[j]-v[j]*tab2[i][j]; // x to element wektora wynikowego (nie potrzebujemy go nigdzie zapisywać)
        if(x != 0) // jeśli któryś element wektora wynikowego nie jest zerowy, to wektor też nie będzie zerowy, więc nie trzeba sprawdzać pozostałych
        {
            ok = 0;
            break;
        }
    }
    puts(ok ? "TAK" : "NIE");
}