funkcionálně.cz

Přední český blog o funkcionálním programování, kde se o funkcionálním programování nepíše
««« »»»

Radix merge sort

25. 5. 2017

Nedávno jsem psal o divoké hybridizaci řadících algoritmů, která mě dovedla k mergeselectu. Teď budu v započaté cestě pokračovat.

Metoda sledování os a hledání průsečíků ukázala jedno nepokryté místo: Radixovou obdobu merge sortu. Quicksort má radixovou verzi (radixsort) a proto by mělo být možné zkřížit merge sort s radixovým způsobem života a vyprodukovat řadídí algoritmus, který stejně jako radixsort postrádá jeden logaritmus z asymptotické složitosti. Výsledek se v retrospektivě zdá natolik očividný, že se mi o něm skoro ani nechce psát. Wikipedie ho přímo nezmiňuje, ale chybí ji k tomu jen malý krůček.

Začnu však popořadě:

Merge sort za běhu vypadá nějak takhle:

Jde o strom slučování seřazených polí do větších seřazených polí. Logaritmus je ve složitosti proto, že hloubka stromu je logaritmická vzhledem k počtu slučovaných polí (které budu označovat jako m na rozdíl od celkového počtu řazených prvků, který dostane klasické jméno n).

Klasický merge sort materializuje všechna dočasná pole. To však není nutné. Strom slučovacích operátorů můžu nahradit stromem operátorů, které extrahují minimum ze svých podstromů, jak ukazuje další obrázek.

Výsledkem je líná verze merge sortu, která v čase O(m) sestaví strom MIN operátorů a potom v čase O(log(m)) extrahuje nejmenší element bez jakékoli materializace mezikroků. Líný merge sort není soubor míst do kterých se zapisuje, ale série trubek, kterými data protékají.

Ve funkcionálních jazycích s líným vyhodnocováním se takto chová merge sort napsaný tradičním stylem. Vstupem jsou líné spojové seznamy, které jsou vyhodnocovány jen tak daleko, jak je nutné. Ze slučovacích mezikroků se stane thunk / suspended computation (jako vždycky nemám ponětí o české terminologii) a jen jedna cesta stromem je přinucena k vyhodnocení s každým extrahovaným elementem. Mám za to, že řazení v Haskellu standardně používá verzi merge sortu právě proto, jak dobře spolupracuje s líným vyhodnocováním.2

Asymptotické časy líné verze vypadají velice povědomě, skoro jako kdyby šlo o haldu. Tu je taky možné vybudovat v lineárním čase a pak extrahovat minimum (a následně ji rebalancovat) v čase logaritmickém.

Min-heap je struktura, která dělá přesně to, co v tomto případě potřebuji. V klasickém merge sortu vezmu dvě pole a porovnám jejich čela, abych zjistil, které z nich je menší. Když slučuji víc polí, musím také najít nejmenší čelo mezi několika kandidáty - tedy to, co umí min-heap.

Když to dotáhnu do krajnosti a zvolím že m = n (všechna vstupní pole jsou singulární) podpůrná datová struktuta úplně převezme otěže a dostanu heapsort. To ale odbočuji.

Min-heap nebo línou konstrukci je možné nahradit binárním vyhledávacím stromem, i když asymptotická analýza ukazuje, že výsledek je morálně sporný. Konstrukce BST trvá O(mlog(m)) a extrakce minima O(log(m)).

V tomto okamžiku už všem musí být jasné, jak vypadá hledaný radix merge sort, žeano? Když vyměním min-heap za trie (neboli prefixový strom) dostanu kýžený algoritmus.

Pro klíče konstantní délky b má trie hloubku O(1). Přesné detaily záleží, na jaké kusy klíče nasekám: zdali na bajty, bity nebo jiné fragmenty.

Algoritmus pracuje tak, že vezme čela seřazených polí a vloží je do trie. Pak ji začne sekvenčně procházet. První prvek, na který narazí, je ten nejmenší. Zapíše ho do výstupu, nové čelo pole, ze kterého prvek pochází, zas vloží do trie a celý postup opakuje. Na diagramu je hledání nejbližšího prvku vyznačeno tečkovaně.

Když mám klíče délky b bitů a rozdělím je na úseky po f bitech, musím v nejhorším případě iterovat 2^f * (2b/f - 1) pozic stromu, než objevím následující nejmenší prvek. 2^f je velikost uzlu a 2b/f-1 představuje cestu z aktuálního listu ke kořeni a pak k dalšímu listu. Dá se předpokládat, že často bude nejbližší větší element relativně blízko, možná i ve stejném listu.

Teoreticky na přesném vzorci nezáleží. b i f jsou konstanty, singulární operace má tedy složitost O(1) a celé řazení pracuje v lineárním čase.

Když začnu přemýšlet o praxi, existuje několik optimalizací, které můžou všechno značně urychlit. Jedna z nich vychází z poznatku, že iteruji jen vpřed a nikdy se nevracím. V momentě, kdy opustím jeden uzel, už ho nikdy nepoužiji, můžu ho tedy zrecyklovat a uštřit nějaké alokace. Trie může mít v nejhorším případě něco pod b/f * m uzlů (m je počet slučovaných polí).

Další zlepšovák spočívá v použití bitmapy, která označuje neprázdné pozice.

struct node {
  void* children[256];
  uint64_t bitmap[4];
}

Nejbližší obsazený slot pak najdu několika LZCNT instrukcemi (number of leading zeros) zcela bez podmíněných skoků (po odrolování).

min = 256;
for (i = 0; i < 4; i++) {
  m = _lzcnt_u64(bitmap[i]);
  min = MIN(min, m < 64 ? m + i*64 : 256);
}

Masochisté, kteří si libují v instrinsics, můžou udělat to samé pomocí AVX2:

zeros      = _mm256_broadcastb_epi8(0);
emptyBytes = _mm256_cmpeq_epi8(bitmap, zeros);
mask       = _mm256_movemask_epi8(emptyBytes);
pos        = _lzcnt_u32(~mask);
i          = _lzcnt_u32(((uint8_t*)&bitmap)[pos]);
return pos * 8 + i;

(Případné chyby jsem v kódu nechal jako cvičení pro čtenáře.)

Bitmapy mají příjemnou synergii s recyklací uzlů, protože nemusím vynulovat celý uzel, ale jen příslušnou bitmapu. Situaci o něco málo komplikují možné duplikáty, ale to není nic, co by se nedalo vyřešit.

To všechno je moc pěkné, ale k čemu je radix merge sort vlastně užitečný? Rozhodně ne pro sortování neseřazeného pole. Kdybych ho použil tak jak je, algoritmus by zdegeneroval na obyčejnou trie, která, tak jak jsem ji tu prezentoval, není příliš paměťově úsporná. Na druhou stranu mi přijde jako celkem užitečný nástroj pro velice široké slučování mnoha seřazených polí1 , protože složitost nenarůstá s jejich počtem.

Na druhou stranu je třeba stále mít na paměti ošklivé tajemství všech radix-sortů: Pro reprezentaci každé z n rozdílných hodnot potřebuji log2(n) bitů. Radixová složitost n*b tedy není nic víc než n * log(n), kde jsme si řekli, že n neporoste nad praktické meze.

  1. Ta mohla vzniknout obyčejným radix sortem a jsou dostatečně malá, že se vešla do CPU cache pro maximální efektivitu
  2. Na líné řazení se spoléhali autoři Haskellové knihovny pro hledání nejbližších sousedů využívající cover tree.
@kaja47, kaja47@k47.cz, deadbeef.k47.cz, starší články