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.

// 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;
}