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

Zobraz historii verzí Skryj historii verzí

Aktuální verze kuchařky: prosinec 2016

Leták 29-3prosinec 2016(Jirka Setnička)
* Revize zdrojových kódů v Pythonu
* Drobné vylepšení a zpřehlednění formulací
Leták 27-3prosinec 2014(Ondra Hlavatý)
* Přepis programů do Pythonu
- Odstraněno několik vět z úvodu
* Rozdělení některých odstavců
Leták 21-2listopad 2008
* Aktualizovaný slovní odkaz na třídicí kuchařku
* Přeformulována poznámka u recyklace pravidla pro hledání mediánu
Leták 19-2říjen 2006
* Aktualizovaný slovní odkaz na třídicí kuchařku
* Přeformulovaný slovní odkaz na 16-S
Leták 17-4leden 2005(David Matoušek, Martin Mareš)
* První verze

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(N2).

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/2k, 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 log2 N už jsou všechny části nejvýše jednoprvkové, takže se rekurze zastaví. Celkem tedy máme log2 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?

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;
Ilustrace: Zástup hrochů

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í:

t(N) ≤ cN + t(N/5) + t(7/10 ·N).

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:

dN ≤ (c + 1/5·d + 7/10 ·d) ·N.

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

Ilustrace: Hroch s počítadlem

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 10NA+B, kde A a B jsou N-ciferná. Součin dvou takových čísel pak bude (10NA+B)·(10NC+D) = (102NAC + 10N(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)≈ N2, č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:

t(N) = N + 3(N/2 + 3t(N/4)) =
= N + 3/2·N + 9t(N/4) =
= N + 3/2·N + 9/4·N + 27t(N/8) = … =
= N + 3/2·N + … + 3k-1/2k-1·N + 3kt(N/2k).

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

t(N) = N·[ 1 + 3/2 + 9/4 + … + (3/2)k-1] + 3kd.

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

(
3
2
)k-1
3
2
-1
= O((
3
2
)k)

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

3k=2k log2 3 =2 log2log2 3 =(2 log2 n) log2 3 = n log2 3 ≈ n1.58

Konstanta d se nám „schová do O-čka“, takže algoritmus má časovou složitost přibližně O(n1.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í

David Matoušek & Martin Mareš