Recepty z programátorské kuchařky
Rekurze a dynamika

V poslední kuchařce tohoto ročníku se budeme zabývat převážně rekurzí a dynamickým programováním. O čem že je řeč? Rekurzivní funkce je taková funkce, která při svém běhu volá sama sebe. Dynamické programování pak bude technika, kterou často půjde z exponenciálně pomalého rekurzivního algoritmu vyrobit pěkný polynomiální. Ale nepředbíhejme, nejdříve se podíváme na jednoduchý příklad rekurze:

Fibonacciho čísla

Budeme počítat n-té číslo Fibonacciho posloupnosti. To je posloupnost, jejíž první dva členy jsou jedničky a každý další člen je součtem dvou předchozích. Začíná takto:

1123581321345589

Pro nalezení n-tého členu (ten budeme značit Fn) si napíšeme rekurzivní funkci Fibonacci(n), která bude postupovat přesně podle definice: zeptá se sama sebe rekurzivně, jaká jsou dvě předchozí čísla, a pak je sečte. Možná více řekne program:

function Fibonacci(n: Integer): Integer;
begin
  if n <= 2 then
    Fibonacci := 1
  else
    Fibonacci := Fibonacci(n-1) + Fibonacci(n-2)
end;

To, jak funkce volá sama sebe, si můžeme snadno nakreslit třeba pro výpočet čísla F5:

Schéma

Vidíme, že program se rozvětvuje a tvoří strom volání. Všimněme si také, že některé podstromy jsou shodné. Zřejmě to budou ty části, které reprezentují výpočet stejného Fibonacciho čísla – v našem příkladě třeba třetího.

Pokusme se odhadnout časovou složitost Tn naší funkce. Pro n=1 a n=2 funkce skončí hned, tedy v konstantním (řekněme jednotkovém) čase. Pro vyšší n zavolá sama sebe pro dva předchozí členy plus ještě spotřebuje konstantní čas na sčítání:

Tn ≥ Tn-1 + Tn-2 + const,  a proto Tn ≥ Fn.

Tedy na spočítání n-tého Fibonacciho čísla spotřebujeme čas alespoň takový, kolik je ono číslo samo. Ale jak velké takové Fn vlastně je? Můžeme třeba využít toho, že:

Fn = Fn - 1 + Fn - 2 ≥ 2 ·Fn - 2,

z čehož plyne:

Fn ≥ 2n/2.

Funkce Fibonacci má tedy exponenciální časovou složitost, což není nic vítaného. Ovšem jak jsme už řekli, některé výpočty opakujeme stále dokola. Nenabízí se proto nic snazšího, než si tyto mezivýsledky uložit a pak je vytáhnout jako pověstného králíka z klobouku s minimem námahy.

Dynamitika

Bude nám k tomu stačit jednoduché pole Pn prvcích, na počátku inicializované nulami. Kdykoliv budeme chtít spočítat některý člen, nejdříve se podíváme do pole, zda jsme ho již jednou nespočetli. A naopak jakmile hodnotu spočítáme, hned si ji do pole poznamenáme:

var P: array[1..MaxN] of Integer;
function Fibonacci(n: Integer): Integer;
begin
  if P[n] = 0 then
  begin
    if n <= 2 then
      P[n] := 1
    else
      P[n] := Fibonacci(n-1) + Fibonacci(n-2)
  end;
  Fibonacci := P[n]
end;

Podívejme se, jak nyní strom volání nyní:

Schéma

Na každý člen posloupnosti se tentokrát ptáme maximálně dvakrát – k výpočtu ho potřebují dva následující členy. To ale znamená, že funkci Fibonacci zavoláme maximálně 2n-krát, čili jsme touto jednoduchou úpravou zlepšili exponenciální složitost na lineární.

Zdálo by se, že abychom získali čas, museli jsme obětovat paměť, ale to není tak úplně pravda. V prvním příkladu sice nepoužíváme žádné pole, ale při volání funkce si musíme zapamatovat některé údaje, jako je třeba návratová adresa, parametry funkce a její lokální proměnné, a na to samotné potřebujeme určitě paměť lineární s hloubkou vnoření, v našem případě tedy lineární s n.

Určitě vás už také napadlo, že n-té Fibonacciho číslo se dá snadno spočítat i bez rekurze. Stačí prvky našeho pole P plnit od začátku – kdykoliv známe P[1]=F1,… ,Fk=P[k], dokážeme snadno spočítat i P[k+1]=Fk+1:

function Fibonacci(n: Integer): Integer;
var
  P: array[1..MaxN] of Integer;
  I: Integer;
begin
  P[1] := 1;
  P[2] := 1;
  for I := 3 to n do
    P[I] := P[I-1] + P[I-2];
  Fibonacci := P[n]
end;

Zopakujme si, co jsme postupně udělali: nejprve jsme vymysleli pomalou rekurzivní funkci, tu jsme zrychlili zapamatováváním si mezivýsledků a nakonec jsme celou rekurzi „obrátili naruby“ a mezivýsledky počítali od nejmenšího k největšímu, aniž bychom se starali o to, jak se na ně původní rekurze ptala.

V případě Fibonacciho čísel je samozřejmě snadné přijít rovnou na nerekurzivní řešení (a dokonce si všimnout, že si stačí pamatovat jen poslední dvě hodnoty a paměťovou složitost tak zredukovat na konstantní), ale zmíněný obecný postup zrychlování rekurze nebo rovnou řešení úlohy od nejmenších podproblémů k těm největším – obvykle se mu říká dynamické programování – funguje i pro řadu složitějších úloh. Třeba na tuto:

Problém batohu

Je dáno N předmětů o hmotnostech m1,… ,mN a také číslo M (nosnost batohu). Úkolem je vybrat některé z předmětů tak, aby součet jejich hmotností byl co největší, a přitom nepřekročil M. My si popíšeme algoritmus, který tento problém řeší v čase O(MN).

Náš algoritmus bude používat pomocné pole A[0… M] a jeho činnost bude rozdělena do N kroků. Na konci k-tého kroku budou nenulové hodnoty v poli A právě na těch pozicích, které odpovídají součtu hmotností předmětů z nějaké podmnožiny prvních k předmětů. Před prvním krokem (po nultém kroku), jsou všechny hodnoty A[i] pro i>0 nulové a A[0] má nějakou nenulovou hodnotu, řekněme -1. Všimněme si, že kroky algoritmu odpovídají podúlohám, které řešíme: nejdříve vyřešíme podúlohu tvořenou jen prvním předmětem, pak prvními dvěma předměty, prvními třemi předměty, atd.

Popišme si nyní k-tý krok algoritmu. Pole A budeme procházet od konce, tj. od i=M po i=mk. Pokud je hodnota A[i] stále nulová, ale hodnota A[i-mk] je nenulová, změníme hodnotu uloženou v A[i] na k. Rozmysleme si, že po provedení k-tého kroku odpovídají nenulové hodnoty v poli A hmotnostem podmnožin z prvních k předmětů, pokud před jeho provedením nenulové hodnoty odpovídaly hmotnostem podmnožin z prvních k-1 předmětů. Pokud je hodnota A[i] nenulová, pak buď byla nenulová před k-tým krokem (a v tom případě odpovídá hmotnosti nějaké podmnožiny prvních k-1 předmětů) anebo se stala nenulovou v k-tém kroku. Potom ale existuje podmnožina prvních k-1 předmětů, jejíž hmotnost je i-mk, a k té stačí přidat k-tý předmět, abychom našli podmnožinu hmotnosti přesně i. Naopak, pokud lze vytvořit podmnožinu I hmotnosti m, pak I je buď tvořena jen prvními k-1 předměty, a tedy hodnota A[m] je nenulová již před k-tým krokem, anebo k∈I. Potom ale hodnota A[m-mk] je nenulová před k-tým krokem (hmotnost podmnožiny I\{k} je m-mk) a hodnota A[m] se stane nenulovou v k-tém kroku.

Po provedení všech N kroků odpovídají nenulové hodnoty A[i] přesně hmotnostem podmnožin ze všech předmětů, co máme k dispozici. Speciálně největší index i0 takový, že hodnota A[i0] je nenulová, odpovídá hmotnosti nejtěžší podmnožiny předmětů, která nepřekročí hmotnost M. Nalézt jednu množinu této hmotnosti také není obtížné: v A[i0] je uloženo číslo jednoho z předmětů nějaké takové podmnožiny, v A[i0-mA[i0]] číslo dalšího předmětu, atd. Samotný kód našeho algoritmu lze nalézt níže.

Časová složitost algoritmu je O(NM), neboť se skládá z N kroků, z nichž každý vyžaduje čas O(M). Paměťová složitost činí O(N+M), což představuje paměť potřebnou pro uložení pomocného pole A a hmotností daných předmětů.

var N: word; {   počet předmětů   }
    M: word; { hmotnostní omezení }
    hmotnost: array[1..N] of word;
             { hmotnosti daných předmětů }
    A: array[0..M] of integer;
    i, k: word;
begin
  A[0]:=-1;
  for i:=1 to M do A[i]:=0;
  for k:=1 to N do
    for i:=M downto hmotnost[k] do
      if (A[i-hmotnost[k]]<>0) and (A[i]=0) then
        A[i]:=k;
  i:=M;
  while A[i]=0 do i:=i-1;
  writeln('Maximální hmotnost: ',i);
  write('Předměty v množině:');
  while A[i]<>-1 do
    begin
      write(' ',A[i]);
      i:=i-hmotnost[A[i]];
    end;
  writeln;
end.

Na rozmyšlenou: Proč pole A procházíme pozadu a ne popředu?

Nejkratší cesty a Floyd-Warshallův algoritmus

Náš další příklad bude z oblasti grafových algoritmů (o grafech se dočtete například v kuchařce třetí série), ale zkusíme si ho nejdříve říci bez grafů:

Bylo-nebylo-je N měst. Mezi některými dvojicemi měst vedou (obousměrné) silnice, jejichž délky jsou dány na vstupu. Předpokládáme, že silnice se jinde než ve městech nepotkávají (pokud se kříží, tak mimoúrovňově). Úkolem je spočítat nejkratší vzdálenosti mezi všemi dvojicemi měst, tj. délky nejkratší cest mezi všemi dvojicemi měst. Cestou rozumíme posloupnost měst spojených silnicemi a délkou cesty součet délek silnic, které spojují po sobě následující města. [V grafové terminologii tedy máme daný ohodnocený neorientovaný graf a chceme zjistit délky nejkratších cest mezi všemi dvojicemi jeho vrcholů.]

Půjdeme na to následovně: Na začátku si uložíme vzdálenosti mezi městy do dvourozměrného polé D, tj. D[i][j] inicializujeme na vzdálenost z města i do města j. Pokud mezi městy i a j nevede žádná silnice, bude D[i][j]=∞, v praxi tedy nějaké dostatečně velké číslo. V průběhu výpočtu si budeme na pozici D[i][j] udržovat délku nejkratší dosud nalezené cesty.

Samotný algoritmus se skládá z N fází. Na konci k-té fáze bude v D[i][j] uložena délka nejkratší cesty mezi městy ij, která může procházet skrz libovolná z měst 1,… ,k. V průběhu k-té fáze tedy stačí vyzkoušet, zda je mezi městy i a j kratší stávající cesta přes města 1,… ,k-1, jejíž délka je uložena v D[i][j], anebo nová cesta přes město k. Pokud nejkratší cesta prochází přes město k, můžeme si ji rozdělit na nejkratší cestu z i do k a nejkratší cestu z k do j. Délka takové cesty je tedy rovna D[i][k]+D[k][j]. Takže pokud je součet D[i][k]+D[k][j] menší než stávající hodnota D[i][j], nahradíme hodnotu na pozici D[i][j] tímto součtem, jinak ji ponecháme.

Z popisu přímo plyne, že po N-té fázi je na pozici D[i][j] uložena délka nejkratší cesty z města i do města j. Každá z N fází algoritmu vyžaduje čas O(N2), a tím pádem časová složitost celého algoritmu bude O(N3). Vystačíme si s pamětí na uložení pole D, tedy O(N2). Program bude vypadat takto.

var N:word; { počet měst }
    D:array[1..N] of array[1..N] of longint;
      { délky silnic mezi městy, D[i][i]=0,
        místo neexistujících je "nekonečno" }
    i,j,k:word;
begin
  for k:=1 to N do
    for i:=1 to N do
      for j:=1 to N do
        if D[i][k]+D[k][j] < D[i][j] then
	  D[i][j]:=D[i][k] + D[k][j];
end.

Popišme si ještě, jak bychom postupovali, kdybychom kromě vzdáleností mezi městy chtěli nalézt i samotné nejkratší cesty. To lze jednoduše vyřešit například tak, že si navíc budeme udržovat pomocné pole E[i][j] a do něj při změně hodnoty D[i][j] uložíme nejvyšší číslo města na cestě z i do j délky D[i][j] (při změně v k-té fázi je to číslo k). Máme-li pak vypsat nejkratší cestu z i do j, vypíšeme nejprve cestu z i do E[i][j] a pak cestu z E[i][j] do j. Tyto cesty nalezneme stejným (rekurzivním) postupem.

Na rozmyšlenou:

  • Jak by algoritmus fungoval, kdyby silnice byly jednosměrné?
  • Na první pohled nejpřirozenější hodnota, kterou bychom mohli použít pro , je maxint. To ovšem nebude fungovat, protože ∞+∞ přeteče. Stačí maxint div 2?
  • Hodnoty v poli si sice přepisujeme pod rukama, takže by se nám mohly poplést hodnoty z předchozí fáze s těmi z fáze současné. Ale zachrání nás to, že čísla, o která jde, vyjdou v obou fázích stejně. Proč?
  • Popis algoritmu vysloveně svádí k „rejpnutí“: Jak víme, že spojením dvou cest, které provádíme, vznikne zase cesta (tj. že se na ní nemohou nějaké vrcholy opakovat)? Inu, to samozřejmě nevíme, ale všimněte si, že kdykoliv by to cesta nebyla, tak si ji nevybereme, protože původní cesta bez vrcholu k bude vždy kratší nebo alespoň stejně dlouhá … tedy alespoň pokud se v naší zemi nevyskytuje cyklus záporné délky. (Což, pokud bychom chtěli být přesní, musíme přidat do předpokladů našeho algoritmu.)
  • Pozor na pořadí cyklů – program vysloveně svádí k tomu, abychom psali cyklus pro k jako vnitřní … jenže pak samozřejmě nebude fungovat.
Dnešní menu vám servírovali
Martin Mareš a Petr Škoda