Nakon nešto duže pauze uglavnom prouzrokovane ispitnim rokom, vraćamo se programiranju i matematici…

Rekurentne formule predstavljaju poseban oblik zapisivanja vrednosti nekog niza. Jedna od najpoznatijih ovakvih relacija je svakako Fibonačijev niz, u kome vrednost nekog člana niza zavisi od vrednosti prethodna dva.

F_n = F_{n-1} + F_{n-2}, ~ F_0 = 0, ~ F_1 = 1

Ukoliko u programu želimo da izračunamo neki član ovog (ili bilo kog drugog) niza, potrebno je da krenemo od vrednosti koje su nam date i iterativno dođemo do željenog člana, pritom prolazeći kroz sve njegove prethodnike. Naravno, moguće je i rešavanje relacije tako da formula zavisi samo od datih vrednosti i N, ali to bi bilo previše komplikovano za implementaciju.

Iako je ovaj pristup krajnje jednostavan, u nekim slučajevima nije moguće primeniti ga usled problema sa efikasnošću. Linearne homogene rekurentne formule, koje ću obraditi u ovom tekstu, mogu biti vremenski zahtevne ukoliko se traži član se velikim indeksom. Ukoliko svaki član zavisi od 10 prethodnih, a njegov indeks je 100.000.000, ukupno će biti potrebno milijardu sabiranja, što je apsolutno nepotrebno. Postoji i komplikovaniji, ali daleko brži pristup.

Inspiracija za ovaj tekst je nastala tokom rešavanja jednog od zadataka sa ovogodišnjeg BubbleCup-a koji je zahtevao efikasan pristup rekurentnim formulama — Archer’s Travel. Kako indeks člana ide čak do jedne milijarde, jednostavno rešenje nije moglo da se izvrši u okviru vremenskog ograničenja. Srećom, pa postoji rešenje sa logaritamskom složenošću.

Neka nam je data sledeća rekurentna formula:

x_n = a_1 x_{n-1} + a_2 x_{n-2} + \dots + a_k x_{n-k}

Xn zavisi od prethodnih k članova, a vrednosti ai su poznati koeficijenti rekurentne formule. Ovu jednačinu možemo zapisati i na sledeći način upotrebom matrica:

x_n = \begin{bmatrix} a_1 & a_2 & \dots & a_k \end{bmatrix} \begin{bmatrix} x_{n-1} \\ x_{n-2} \\ \dots \\ x_{n-k} \end{bmatrix}

Naravno, kada se ove dve matrice pomnože, rezultat je isti kao u prethodnoj jednačini. Za sledeći korak nama nije dovoljno samo xn, već nam treba i njegovih k-1 prethodnika. Ali zašto ne bismo i njih uvrstili u levu stranu jednakosti i dobili opšti iterativni korak?

\begin{bmatrix} x_{n} \\ x_{n-1} \\ \dots \\ x_{n-k-1} \end{bmatrix} = A \begin{bmatrix} x_{n-1} \\ x_{n-2} \\ \dots \\ x_{n-k} \end{bmatrix}

Sada nam ostaje samo da formiramo kvadratnu matricu A tako da osim xn dobijemo i ostale potrebne članove. Tačnije, prvi red matrice A ćemo popuniti koeficijentima rekurentne formule, a ostatak ćemo ispuniti nulama, uz jedinice po dijagonali (ali ne dijagonali matrice), radi dobijanja željenog proizvoda.

\begin{bmatrix} x_{n} \\ x_{n-1} \\ \dots \\ x_{n-k+2} \\ x_{n-k+1} \end{bmatrix} = \begin{bmatrix} a_1 & a_2 & \dots & a_{k-1} & a_k \\ 1 & 0 & \dots & 0 & 0 \\ 0 & 1 & \dots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \dots & 1 & 0 \end{bmatrix} \begin{bmatrix} x_{n-1} \\ x_{n-2} \\ \dots \\ x_{n-k+1} \\ x_{n-k} \end{bmatrix}

Pošto smo dobili opšti oblik za računanje matrice, možemo ga rekurzivno primeniti da bismo došli do poznatih vrednosti članova niza.

\begin{bmatrix} x_{n} \\ x_{n-1} \\ \dots \\ x_{n-k+1} \end{bmatrix} = A \begin{bmatrix} x_{n-1} \\ x_{n-2} \\ \dots \\ x_{n-k} \end{bmatrix} = A^2 \begin{bmatrix} x_{n-2} \\ x_{n-3} \\ \dots \\ x_{n-k-1} \end{bmatrix} = A^3 \begin{bmatrix} x_{n-3} \\ x_{n-4} \\ \dots \\ x_{n-k-2} \end{bmatrix} = \dots = A^{n-k+1} \begin{bmatrix} x_{k-1} \\ x_{k-2} \\ \dots \\ x_{0} \end{bmatrix}

Odlično! Sveli smo računanje člana xn na samo jedno množenje. Avaj, An-k+1 nam nije poznato, i računanje ove matrice se svodi na još n-k+1 množenja…

Da li je to slučaj? Naravno da nije. Sada nailazimo na korak koji predstavlja najbitniji deo ove “optimizacije” — računanje stepena nekog broja (odnosno matrice, u našem slučaju) se može izvršiti u logaritamskom vremenu.

A^n = \left\{ \begin{array}{l l} A^{n/2}A^{n/2} & \quad \text{ako je n parno} \\ A^{\lfloor n/2\rfloor}A^{\lfloor n/2 \rfloor}A & \quad \text{ako je n neparno} \\ \end{array} \right.

Iako je množenje matrica nešto kompleksniji posao, ušteda je i dalje jako velika – najjednostavniji algoritam za množenje matrica je složenosti k3, gde je k dimenzija kvadratne matrice.

Za kraj, sledi najjednostavniji algoritam za računanje vrednosti nekog člana rekuretne formule implementiran u Javi. Dalje optimizacije su moguće (izbacivanje rekurzije, efikasnije množenje matrica), ali nemaju ni približno veliki uticaj kao što ima uvođenje matrica i njihovog stepenovanja u algoritam.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
// Brojevi se cuvaju u BigDecimal jer mogu postati jako veliki
// Smatra se da je indeks prvog clana niza nula
// Pocetni clanovi se prosledjuju redom x(0), x(1), ...
// Koeficijenti se prosledjuju obrnutim redosledom - prvo ide
// koeficijent uz x(i), pa uz x(i-1) itd.
 
static BigDecimal izracunaj(int indeksClana, double[] pocetniClanovi, double[] koeficijenti)
{
    int k = pocetniClanovi.length;
 
    // Ukoliko prosledjeni nizovi nisu validni
 
    if (k != koeficijenti.length || k < 1)
        return null;
 
    // Ukoliko se trazi clan koji je vec poznat
 
    if (indeksClana < k)
        return new BigDecimal(pocetniClanovi[indeksClana]);
 
    // Sve je u redu, formiramo matricu
 
    BigDecimal[][] matrica = new BigDecimal[k][k];
 
    for (int i = 0; i < k; i++)
        matrica[0][i] = new BigDecimal(koeficijenti[i]);
 
    for (int i = 1; i < k; i++)
        for (int j = 0; j < k; j++)
            if (i - 1 == j)
                matrica[i][j] = new BigDecimal(1);
            else
                matrica[i][j] = new BigDecimal(0);
 
    // Stepenujemo matricu - A^(indeks-k+1)
 
    matrica = stepenujRekurzivno(indeksClana - k + 1, matrica, k);
 
    // Racunamo x(indeks)
 
    BigDecimal rezultat = new BigDecimal(0);
 
    for (int i = 0; i < k; i++)
        rezultat = rezultat.add(matrica[0][i].multiply(new BigDecimal(pocetniClanovi[k - i - 1])));
 
    return rezultat;
}
 
static BigDecimal[][] stepenujRekurzivno(int stepen, BigDecimal[][] matrica, int dimenzije)
{
    if (stepen == 1)
        return matrica;
 
    BigDecimal[][] tempMatrica = stepenujRekurzivno(stepen / 2, matrica, dimenzije);
 
    // Kvadriramo dobijenu matricu
 
    BigDecimal[][] novaMatrica = new BigDecimal[dimenzije][dimenzije];
 
    for (int i = 0; i < dimenzije; i++)
        for (int j = 0; j < dimenzije; j++)
        {
            novaMatrica[i][j] = new BigDecimal(0);
 
            for (int t = 0; t < dimenzije; t++)
            {
                novaMatrica[i][j] = novaMatrica[i][j].add(tempMatrica[i][t].multiply(tempMatrica[t][j]));
            }
        }
 
    // Mnozimo jos jednom ako je stepen neparan
 
    if (stepen % 2 != 0)
    {
        BigDecimal[][] novaTempMatrica = new BigDecimal[dimenzije][dimenzije];
 
        for (int i = 0; i < dimenzije; i++)
            for (int j = 0; j < dimenzije; j++)
            {
                novaTempMatrica[i][j] = new BigDecimal(0);
 
                for (int t = 0; t < dimenzije; t++)
                    novaTempMatrica[i][j] = novaTempMatrica[i][j].add(matrica[i][t].multiply(novaMatrica[t][j]));
            }
        novaMatrica = novaTempMatrica;
    }
 
    // Vracamo matricu
 
    return novaMatrica;
}