Recepty z programátorské kuchařky
Rozděl a panuj

Dnešní díl programátorské kuchařky se bude zabývat algoritmy založenými na metodě Rozděl a panuj. Slušelo by se začít tím, jaká je myšlenka této metody: Často se setkáme s úlohami, které lze snadno rozdělit na nějaké menší úlohy a z jejich výsledků zase snadno složit výsledek původní velké úlohy. Přitom menší úlohy můžeme řešit opět týmž algoritmem (zavoláme si ho rekurzivně), leda by již byly tak maličké, že dokážeme odpovědět triviálně bez jakéhokoliv počítání. Zkrátka jak říkali staří římští císařové: Divide et impera. Uveďme si pro začátek jeden staronový příklad:

Quicksort

QuickSort (alias QS) je algoritmus pro třídění posloupnosti prvků. Už jste si o něm mohli přečíst v kuchařce o třídění. Tentokrát se na něj podíváme trochu podrobněji a navíc nám poslouží jako ingredience pro další algoritmy.

QS v každém svém kroku zvolí nějaký prvek (budeme mu říkat pivot) a přerovná prvky v posloupnosti tak, aby napravo od pivota byly pouze prvky větší než pivot a nalevo pouze menší. Pokud se vyskytnou prvky rovné pivotu, můžeme si dle libosti vybrat jak levou, tak pravou stranu posloupnosti, funkčnost algoritmu to nijak neovlivní. Tento postup pak rekurzivně zopakujeme zvlášť pro prvky nalevo a zvlášť pro prvky napravo od pivota, a tak získáme setříděnou posloupnost.

Implementací QS je mnoho a mimo jiné se liší způsobem volby pivota. My si předvedeme jinou, než jsme ukazovali v třídící kuchařce (hlavně proto, že se nám z ní pak budou snadno odvozovat další algoritmy), a pro jednoduchost budeme jako pivota volit poslední prvek zkoumaného úseku:

pole = [1, 2, 8, 42, 9, 17, -5, 20, 2]

# Přerovnej pole od levého do pravého indexu
def prerovnej(pole, L, P):
    pivot = pole[P - 1]
    # i je nejlevější nepřerovnaný prvek
    i = L
    # j je aktuální probíraný prvek
    for j in range(L, P - 1):
        if (pole[j] <= pivot):
            # Prohodíme prvek s nejlevějším
            pole[i],pole[j] = pole[j],pole[i]
            i += 1
    # Dáme pivota na správné místo
    pole[P - 1],pole[i] = pole[i],pole[P - 1]
    return i

def quicksort(pole, levy_index, pravy_index):
    if (levy_index >= pravy_index):
        return
    # Přerovnáme úsek a najdeme pivota...
    p = prerovnej(pole, levy_index, pravy_index)
    # ... a zavoláme se rekurzivně na podúseky
    quicksort(pole, levy_index, p)
    quicksort(pole, p + 1, pravy_index)

{budeme třídit takováto pole}
type Pole=array[1..MaxN] of Integer;

{přerovnávací procedura pro úsek a[l..r]}
function prer(a:Pole; l,r:Integer):Integer;
var i,j,x,q:Integer;
begin
  {pivotem se stane poslední prvek úseku}
  x:=a[r];                      {hodnota pivota}
  i:=l-1;  {a[i] bude vždy poslední <= pivotovi}

  {samotné přerovnávání}
  for j:=l to r-1 do
    if a[j]<=x then   {právě probíraný prvek     }
    begin             {menší/rovný hodnotě pivota}
      Inc(i);         {pak zvyš ukazatel         }
      q:=a[j];        {a proveď přerovnání prvku }
      a[j]:=a[i];
      a[i]:=q;
    end;

  {nakonec přesuneme pivota za poslední <=}
  q:=a[r];
  a[r]:=a[i+1];
  a[i+1]:=q;
  prer:=i+1;      {vrátíme novou pozici pivota}
end;

{hlavní třídící procedura, třídí a[l..r]}
procedure QuickSort(a:Pole; l,r:Integer);
var m:Integer;
begin
  if l<r then            {máme ještě co dělat?}
  begin
    m:=prer(l,r);   {přerovnej, m pozice pivota}
    QuickSort(l,m-1);     {setřiď prvky napravo}
    QuickSort(m+1,r);      {setřiď prvky nalevo}
  end;
end;

Bohužel volit pivota právě takto je docela nešikovné, protože se nám snadno může stát, že si vybereme nejmenší nebo největší prvek v úseku (rozmyslete si, jak by vypadala posloupnost, ve které to nastane pokaždé), takže dostaneme posloupnost délky N, rozdělíme ji na úseky délek N-1 a 1, načež pokračujeme s úsekem délky N-1, ten rozdělíme na N-2 a 1, atd. Přitom pokaždé na přerovnání spotřebujeme čas lineární s velikostí úseku, celkem tedy O(N+(N-1)+(N-2)+… +1) = O(N²).

Na druhou stranu pokud bychom si za pivota vybrali vždy medián z právě probíraných prvků (tj. prvek, který by se v setříděné posloupnosti nacházel uprostřed; pro sudý počet prvků zvolíme libovolný z obou prostředních prvků), dosáhneme daleko lepší složitosti O(N log N). Ale raději si to dokažme pořádně:

Přerovnávací část algoritmu běží v čase lineárním vůči počtu prvků, které máme přerovnat. V prvním kroku QS pracujeme s celou posloupností, čili přerovnáme celkem N prvků. Následuje rekurzivní volání pro levou a pravou část posloupnosti (obě dlouhé (N-1)/2 ±1); přerovnávání v obou částech dohromady trvá opět O(N) a vzniknou tím části dlouhé nejvýše N/4. Zanoříme-li se v rekurzi do hloubky k, pracujeme s částmi dlouhými nejvýše N/2^k, které dohromady dají nejvýše N (všechny části dohromady dají prvky vstupní posloupnosti bez těch, které jsme si už zvolili jako pivoty). V hloubce ⌈ log₂ N⌉ už jsou všechny části nejvýše jednoprvkové, takže se rekurze zastaví. Celkem tedy máme ⌈ log₂ N⌉ hladin (hloubek) a na každé z nich trávíme lineární čas, dohromady O(N log N).

V tomto důkazu jsme se ale dopustili jednoho podvodu: Zapomněli jsme na to, že také musíme medián umět najít. Jak z této nepříjemné situace ven?

• Naučit se počítat medián. Ale jak?
• Spokojit se se „lžimediánem“: Kdybychom si místo mediánu vybrali libovolný prvek, který bude v setříděné posloupnosti „v prostřední polovině“ (čili alespoň čtvrtina prvků bude větší a alespoň čtvrtina menší než on), získáme také složitost O(N log N), neboť úsek délky N rozložíme na úseky, které budou mít délky nejvýše (1-1/4)·N, takže na k-té hladině budou úseky délek nejvýše (1-1/4)^k ·N, čili hladin bude maximálně log_1-1/4 N = O(log N). Místo 1/4 by fungovala i libovolná jiná konstanta mezi nulou a jedničkou, ale ani to nám nepomůže k tomu, abychom uměli lžimedián najít.
• Recyklovat pravidlo typu „vezmi poslední prvek“ a jen ho trochu vylepšit. To bohužel nebude fungovat, protože pokud budeme při výběru pivota hledět jenom na konstantní počet prvků, bude poměrně snadné přijít na vstup, pro který toto pravidlo bude dávat kvadratickou složitost, i když obvykle půjde dokázat, že takových vstupů je „málo“. (Proto se také QS často implementuje právě s náhodnou volbou pivota.)
• Volit pivota náhodně ze všech prvků zkoumaného úseku. K náhodné volbě samozřejmě potřebujeme náhodný generátor a s těmi je to svízelné, ale zkusme na chvíli věřit, že jeden takový máme nebo alespoň že máme něco s podobnými vlastnostmi. Jak nám to pomůže? Náhodně zvolený pivot nebude sice přesně uprostřed, ale s pravděpodobností 1/2 to bude lžimedián, takže po průměrně dvou hladinách se ke lžimediánu dopracujeme. Proto časová složitost takovéhoto randomizovaného QS bude v průměru 2-krát větší, než lžimediánového QS, čili v průměru také O(N log N). Jednoduše řečeno, zatímco fixní pravidlo nám dalo dobrý čas pro průměrný vstup, ale existovaly vstupy, na kterých bylo pomalé, randomizování nám dává dobrý průměrný čas pro všechny možné vstupy.

Hledání k-tého nejmenšího prvku

Nad QuickSortem jsme zvítězili, ale současně jsme při tom zjistili, že neumíme rychle najít medián. To tak nemůžeme nechat, a proto rovnou zkusíme vyřešit obecnější problém: najít k-tý nejmenší prvek (medián dostáváme pro k=⌊N/2 ⌋).

První řešení této úlohy se nabízí samo. Načteme posloupnost do pole, prvky pole setřídíme nějakým rychlým algoritmem a kýžený k-tý nejmenší prvek nalezneme na k-té pozici v nyní již setříděném poli. Má to však jeden háček. Pokud prvky, které máme na vstupu, umíme pouze porovnat, pak nedosáhneme lepší časové složitosti (a to ani v průměrném případě) než O(N log N) – rychleji prostě třídit nelze, důkaz můžete najít například v třídící kuchařce.

O něco rychlejší řešení je založeno na výše zmíněném algoritmu QuickSort (často se mu proto říká QuickSelect). Opět si vybereme pivota a posloupnost rozdělíme na prvky menší než pivot, pivota a prvky větší než pivot (pro jednoduchost budeme předpokládat, že žádné dva prvky posloupnosti nejsou stejné).

Pokud se pivot nalézá na k-té pozici, je to hledaný k-tý nejmenší prvek posloupnosti, protože právě k-1 prvků je menších. Zbývají dva případy, kdy tomu tak není. Pakliže je pozice pivota v posloupnosti větší než k, pak se hledaný prvek nalézá nalevo od pivota a postačí rekurzivně najít k-tý nejmenší prvek mezi prvky nalevo. V opačném případě, kdy je pozice pivota menší než k, je hledaný prvek v posloupnosti napravo od pivota. Mezi těmito prvky však nebudeme hledat k-tý nejmenší prvek, ale (k-p)-tý nejmenší prvek, kde p je pozice pivota v posloupnosti.

Časovou složitost rozebereme podobně jako u QuickSortu. Nešikovná volba pivota dává opět v nejhorším případě kvadratickou složitost. Pokud bychom naopak volili za pivota medián, budeme nejprve přerovnávat N prvků, pak jich zbude nejvýše N/2, pak nejvýše N/4 atd., což dohromady dává složitost O(N+N/2+N/4+… +1) = O(N). Pro lžimedián dostaneme rovněž lineární složitost a opět stejně jako u QS můžeme nahlédnout, že náhodnou volbou pivota dostaneme v průměru stejný čas jako se lžimediánem.

Program bude velmi jednoduchý, využijeme-li přerovnávací proceduru od QS:

def kty(pole, k, L, P):
    pivot = prerovnej(pole, L, P)

    if (k == pivot):
        return pole[pivot]
    if (k < pivot):
        return kty(pole, k, L, pivot)
    else:
        return kty(pole, k, pivot + 1, P)

function kty(var a:Pole; l,r,k:Integer):Integer;
var x,z:Integer;
begin
  x:=prer(a,l,r); {přerovnej, x je pozice pivota}
  z:=x-l+1;     {pozice pivota vzhledem k [l..r]}
  if k=z then
    kty:=a[x]            {k-tý nejmenší je pivot}
  else if k<z then
    kty:=kty(a,l,x-1,k) {k-tý nejmenší je nalevo}
  else
    kty:=kty(a,x+1,r,k-z);              {napravo}
end;

k-tý nejmenší podruhé, tentokrát lineárně a bez náhody

Existuje však algoritmus, který řeší naši úlohu lineárně, a to i v nejhorším případě. Je založený na ďábelském triku: zvolit vhodného pivota (jak ukážeme, bude to jeden ze lžimediánů) rekurzivním voláním téhož algoritmu. Zařídíme to takto:

1. Pokud jsme dostali méně než 6 prvků, použijeme nějaký triviální algoritmus, například si posloupnost setřídíme a vrátíme k-tý prvek setříděné posloupnosti.
2. Rozdělíme prvky posloupnosti na pětice; pokud není počet prvků dělitelný pěti, poslední pětici necháme nekompletní.
3. Spočítáme medián každé pětice. To můžeme provést například rekurzivním zavoláním celého našeho algoritmu, čili v důsledku tříděním. (Také bychom si mohli pro 5 prvků zkonstruovat rozhodovací strom s nejmenším možným počtem porovnání, což je rychlejší, ale jednak pouze konstanta-krát, jednak je to daleko pracnější.)
4. Máme tedy N/5 mediánů. V nich rekurzivně najdeme medián m (označíme mediány pětic za novou posloupnost a na ní začneme opět od prvního bodu).
5. Přerovnáme vstupní posloupnost po quicksortovsku a jako pivota použijeme prvek m. Po přerovnání je pivot, podobně jako v předchozím algoritmu, na (z+1)-ní pozici v posloupnosti, kde z je počet prvků s menší hodnotou, než má pivot.
6. Opět, podobně jako u předchozího algoritmu, pokud je k=z+1, pak je právě pivot m k-tým nejmenším prvkem posloupnosti. V případě, že tomu tak není a k<z+1, budeme hledat k-tý nejmenší prvek mezi prvními z členy posloupnosti, v opačném případě, kdy k>z+1, budeme hledat (k-z+1)-ní nejmenší prvek mezi posledními n-z-1 prvky.

# Příprava pro přerovnávání z QuickSortu
# -> chceme pivota jako poslední prvek
def prerovnej_podle(pole, L, P, podle):
    q = L
    while (pole[q] != podle):
        q += 1
    pole[q], pole[P-1] = pole[P-1], pole[q]
    return prerovnej(pole, L, P)

def kty(pole, k, L, P):
    pocet = P - L
    # Jednoduché případy
    if (pocet <= 1):
        return pole[L]
    if (pocet <= 5):
        quicksort(pole, L, P)
        return pole[k]

    # Rozdělení na pětice
    petic = (pocet + 4) // 5
    mediany = [0] * petic
    for i in range(0, pocet, 5):
        if (i + 5 > pocet):
            # Ignorujeme neúplnou pětici
            break
        quicksort(pole, L + i, L + i + 5)
        mediany[i // 5] = pole[i + 2]

    # Nalezneme medián mediánů pětic
    median = kty(mediany, petic//2, 0, petic)
    pivot = prerovnej_podle(pole, L, P, median)

    if (pivot == k):
        return median
    if (k < pivot):
        return kty(pole, k, L, pivot)
    else:
        return kty(pole, k, pivot + 1, P)

{potřebujeme přerovnávací funkci, která
 dostane hodnotu pivota jako parametr}
function prerp(var a:Pole; l, r, m:Integer):Integer;
var q,p:Integer;
begin
  {nalezneme pozici pivota}
  p:=l;
  while a[p]<>m do
    inc(p);
  {pivota prohodíme s posledním prvkem}
  q:=a[p]; a[p]:=a[r]; a[r]:=q;
  {a zavoláme původní přerovnávací fci}
  prerp := prer(a,l,r);
end;

{hledání k-tého nejmenšího prvku z a[l..r]}
function kth(var a:Pole; l, r, k:Integer):Integer;
var medp:Pole;             {pole pro mediány pětic}
    i,j,q,x,pocet,m,z:Integer;
begin
  pocet:=r-l+1;          {s kolika prvky pracujeme}

  if pocet<=1 then             {pouze jeden prvek?}
    kth:=a[l]            {výsledek nemůže být jiný}
  else if pocet<6 then begin     {méně než 6 prvků}
    QuickSort(a,l,r);
    kth:=a[l+k-1];
    end
  else begin         {mnoho prvků, jde to tuhého}
    {rozdělíme prvky do pětic}
    q:=1;                 {zatím máme jednu pětici}
    i:=l;                 {levý okraj první pětice}
    j:=i+4;              {pravý okraj první pětice}
    while j<=r do begin    {procházíme celé pětice}
      QuickSort(a,i,j);
      medp[q]:=a[i+2];              {medián pětice}
      Inc(q);                    {zvyš počet pětic}
      Inc(i,5);          {nastav levý okraj pětice}
      Inc(j,5);         {nastav pravý okraj pětice}
    end;
    {případnou neúplnou pětici můžeme ignorovat}

    {najdeme medián mediánů pětic}
    m:=kth(medp,1,q-1,q div 2);

    {přerovnej a zjisti, kde skončil pivot}
    x:=prer(a,l,r,m);
    z:=x-l+1;            {pozice vzhledem k [l..r]}
    if k=z then
      kth:=m               {k-tý nejmenší je pivot}
    else if k<z then
      kth:=kth(a,l,x-1,k)    {k-tý nejmenší nalevo}
    else
      kth:=kth(a,x+1,r,k-z);              {napravo}
  end;
end;

Zbývá dokázat, že tato dvojitá rekurze má slíbenou lineární složitost. Zkusme se proto podívat, kolik prvků posloupnosti po přerovnání je větších než prvek m. Všech pětic je N/5 a alespoň polovina z nich (tedy N/10) má medián menší než m. V každé takové pětici pak navíc najdeme dva prvky menší než medián pětice (takže jsou zde tři prvky menší, než m). Celkem tak existuje alespoň 3/10 ·N prvků menších než m. Větších tedy může být maximálně 7/10 ·N. Symetricky ukážeme, že i menších prvků může být nejvýše 7/10 ·N.

Rozdělení na pětice, hledání mediánů pětic a přerovnávání trvá lineárně, tedy nejvýše cN kroků pro nějakou konstantu c>0. Pak už algoritmus pouze dvakrát rekurzivně volá sám sebe: nejprve pro N/5 mediánů pětic, pak pro ≤ 7/10 ·N prvků před/za pivotem. Pro celkovou časovou složitost t(N) našeho algoritmu tedy platí:

Nyní zbývá tuto rekurzivní nerovnici vyřešit, což provedeme drobným úskokem: uhodneme, že výsledkem bude lineární funkce, tedy že t(N)=dN pro nějaké d>0. Dostaneme:

To platí např. pro d=10c, takže opravdu t(N)=O(N).

Násobení dlouhých čísel

Dalším pěkným příkladem na rozdělování a panování je násobení dlouhých čísel – tak dlouhých, že se už nevejdou do integeru, takže s nimi musíme počítat po číslicích (ať už v jakékoliv soustavě – teď zvolíme desítkovou, často se hodí třeba 256-ková). Klasickým „školním“ algoritmem pro násobení na papíře to zvládneme na kvadratický počet operací, zde si předvedeme efektivnější způsob.

Libovolné 2N-ciferné číslo můžeme zapsat jako 10^NA+B, kde A a B jsou N-ciferná. Součin dvou takových čísel pak bude (10^NA+B)·(10^NC+D) = (10^2NAC + 10^N(AD+BC) + BD). Sčítat dokážeme v lineárním čase, násobit mocninou deseti také (dopíšeme příslušný počet nul na konec čísla), N-ciferná čísla budeme násobit rekurzivním zavoláním téhož algoritmu. Pro časovou složitost tedy bude platit t(N)=cN + 4t(N/2). Nyní tuto rovnici můžeme snadno vyřešit, ale ani to dělat nebudeme, neboť nám vyjde, že t(N)≈ N², čili jsme si oproti původnímu algoritmu vůbec nepomohli.

Přijde trik. Místo čtyř násobení čísel poloviční délky nám budou stačit jen tři: spočteme AC, BD a (A+B)·(C+D) = AC+AD+BC+BD, přičemž pokud od posledního součinu odečteme AC a BD, dostaneme přesně AD+BC, které jsme předtím počítali dvěma násobeními. Časová složitost nyní bude t(N)=c' N + 3t(N/2). (Konstanta c' je o něco větší než c, protože přibylo sčítání a odčítání, ale stále je to konstanta. My si ovšem zvolíme jednotku času tak, aby bylo c'=1, a ušetříme si tak spoustu psaní.)

Jak naši rovnici vyřešíme? Zkusíme ji dosadit do sebe samé a pozorovat, co se bude dít:

Pokud zvolíme k= log₂ N, vyjde N/2^k=1, čili t(N/2^k)=t(1)=d, kde d je nějaká konstanta. To znamená, že:

Výraz v hranaté závorce je součet prvních k členů geometrické řady s kvocientem (neboli podílem dvou po sobě jdoucích prvků) 3/2. Tuto geometrickou řadu si můžeme sečíst jako:

Tato funkce však roste pomaleji než zbylý člen 3^kd, takže ji klidně můžeme zanedbat a zabývat se pouze oním posledním členem:

Konstanta d se nám „schová do O-čka“, takže algoritmus má časovou složitost přibližně O(n^1.58). Existují i rychlejší algoritmy se složitostí až O(n log n), ale ty jsou mnohem ďábelštější a pro malá n se to sotva vyplatí.

Program si pro dnešek odpustíme, šetřímeť naše lesy.

K zamyšlení

• Při hledání k-tého nejmenšího prvku jsme předpokládali, že všechny prvky jsou různé. Prohlédněte si algoritmy pozorně a rozmyslete si, že budou fungovat i bez toho. Opravdu?
• Proč jsme zvolili zrovna pětice? Jak by to dopadlo pro trojice? A jak pro sedmice? Fungoval by takový algoritmus? Byl by také lineární?
• Ve výpočtu t(N) jsme si nedali pozor na neúplné pětice a také jsme předpokládali, že pětic je sudý počet. Ono se totiž nic zlého nemůže stát. Jak se to snadno nahlédne? Proč nestačí na začátku doplnit vstup „nekonečny“ na délku, která je mocninou deseti?
• Kdybychom neuhodli, že t(N) je lineární, jak by se na to dalo přijít?
• Ještě jednou QS: Představte si, že budujete binární vyhledávací strom vkládáním prvků v náhodném pořadí. Obecně nemusí být vyvážený, ale v průměru v něm půjde vyhledávat v čase O(log N). Žádný div: Stromy, které nám vzniknou, odpovídají přesně možným průběhům QuickSortu.

David Matoušek & Martin Mareš