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 — k47

Ne­dávno jsem psal o divoké hyb­ri­di­zaci řa­dí­cích al­go­ritmů, která mě do­vedla k mer­ge­se­lectu. Teď budu v za­po­čaté cestě po­kra­čo­vat.

Metoda sle­do­vání os a hle­dání prů­se­číků uká­zala jedno ne­po­kryté místo: Ra­di­xo­vou obdobu merge sortu. Quicksort má ra­di­xo­vou verzi (ra­di­x­sort) a proto by mělo být možné zkří­žit merge sort s ra­di­xo­vým způ­so­bem života a vy­pro­du­ko­vat řadídí al­go­rit­mus, který stejně jako ra­di­x­sort po­strádá jeden lo­ga­rit­mus z asympto­tické slo­ži­tosti. Vý­sle­dek se v re­trospek­tivě zdá na­to­lik oči­vidný, že se mi o něm skoro ani nechce psát. Wi­ki­pe­die ho přímo ne­zmi­ňuje, ale chybí ji k tomu jen malý krůček.

Začnu však po­po­řadě:

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

Jde o strom slu­čo­vání se­řa­ze­ných polí do vět­ších se­řa­ze­ných polí. Lo­ga­rit­mus je ve slo­ži­tosti proto, že hloubka stromu je lo­ga­rit­mická vzhle­dem k počtu slu­čo­va­ných polí (které budu ozna­čo­vat jako m na rozdíl od cel­ko­vého počtu řa­ze­ných prvků, který do­stane kla­sické jméno n).

Kla­sický merge sort ma­te­ri­a­li­zuje všechna do­časná pole. To však není nutné. Strom slu­čo­va­cích ope­rá­torů můžu na­hra­dit stro­mem ope­rá­torů, které ex­tra­hují mi­ni­mum ze svých pod­stromů, jak uka­zuje další ob­rá­zek.

Vý­sled­kem je líná verze merge sortu, která v čase O(m) se­staví strom MIN ope­rá­torů a potom v čase O(log(m)) ex­tra­huje nejmenší ele­ment bez ja­ké­koli ma­te­ri­a­li­zace me­zi­kroků. Líný merge sort není soubor míst do kte­rých se za­pi­suje, ale série trubek, kte­rými data pro­té­kají.

Ve funk­ci­o­nál­ních ja­zy­cích s líným vy­hod­no­co­vá­ním se takto chová merge sort na­psaný tra­dič­ním stylem. Vstu­pem jsou líné spo­jové se­znamy, které jsou vy­hod­no­co­vány jen tak daleko, jak je nutné. Ze slu­čo­va­cích me­zi­kroků se stane thunk / suspen­ded com­pu­tation (jako vždycky nemám ponětí o české ter­mi­no­lo­gii) a jen jedna cesta stro­mem je při­nu­cena k vy­hod­no­cení s každým ex­tra­ho­va­ným ele­men­tem. Mám za to, že řazení v Haskellu stan­dardně po­u­žívá verzi merge sortu právě proto, jak dobře spo­lu­pra­cuje s líným vy­hod­no­co­vá­ním.2

Asympto­tické časy líné verze vy­pa­dají velice po­vě­domě, skoro jako kdyby šlo o haldu. Tu je taky možné vy­bu­do­vat v li­ne­ár­ním čase a pak ex­tra­ho­vat mi­ni­mum (a ná­sledně ji re­ba­lan­co­vat) v čase lo­ga­rit­mic­kém.

Min-heap je struk­tura, která dělá přesně to, co v tomto pří­padě po­tře­buji. V kla­sic­kém merge sortu vezmu dvě pole a po­rov­nám jejich čela, abych zjis­til, které z nich je menší. Když slu­čuji víc polí, musím také najít nejmenší čelo mezi ně­ko­lika kan­di­dáty – tedy to, co umí min-heap.

Když to do­táhnu do kraj­nosti a zvolím že m = n (všechna vstupní pole jsou sin­gu­lární) pod­půrná datová struk­tuta úplně pře­vezme otěže a do­stanu he­ap­sort. To ale od­bo­čuji.

Min-heap nebo línou kon­strukci je možné na­hra­dit bi­nár­ním vy­hle­dá­va­cím stro­mem, i když asympto­tická ana­lýza uka­zuje, že vý­sle­dek je mo­rálně sporný. Kon­strukce BST trvá O(mlog(m)) a ex­trakce minima O(log(m)).

V tomto oka­mžiku už všem musí být jasné, jak vypadá hle­daný radix merge sort, žeano? Když vy­mě­ním min-heap za trie (neboli pre­fi­xový strom) do­stanu kýžený al­go­rit­mus.

Pro klíče kon­stantní délky b má trie hloubku O(1). Přesné de­taily záleží, na jaké kusy klíče na­se­kám: zdali na bajty, bity nebo jiné frag­menty.

Al­go­rit­mus pra­cuje tak, že vezme čela se­řa­ze­ných polí a vloží je do trie. Pak ji začne sek­venčně pro­chá­zet. První prvek, na který narazí, je ten nejmenší. Zapíše ho do vý­stupu, nové čelo pole, ze kte­rého prvek po­chází, zas vloží do trie a celý postup opa­kuje. Na di­a­gramu je hle­dání nej­bliž­šího prvku vy­zna­čeno teč­ko­vaně.

Když mám klíče délky b bitů a roz­dě­lím je na úseky po f bitech, musím v nej­hor­ším pří­padě ite­ro­vat 2^f * (2b/f – 1) pozic stromu, než ob­je­vím ná­sle­du­jící nejmenší prvek. 2^f je ve­li­kost uzlu a 2b/f-1 před­sta­vuje cestu z ak­tu­ál­ního listu ke kořeni a pak k dal­šímu listu. Dá se před­po­klá­dat, že často bude nej­bližší větší ele­ment re­la­tivně blízko, možná i ve stej­ném listu.

Te­o­re­ticky na přes­ném vzorci ne­zá­leží. bf jsou kon­stanty, sin­gu­lární ope­race má tedy slo­ži­tost O(1) a celé řazení pra­cuje v li­ne­ár­ním čase.

Když začnu pře­mýš­let o praxi, exis­tuje ně­ko­lik op­ti­ma­li­zací, které můžou všechno značně urych­lit. Jedna z nich vy­chází z po­znatku, že ite­ruji jen vpřed a nikdy se ne­vra­cím. V mo­mentě, kdy opus­tím jeden uzel, už ho nikdy ne­po­u­žiji, můžu ho tedy zrecyklo­vat a uštřit nějaké alo­kace. Trie může mít v nej­hor­ším pří­padě něco pod b/f * m uzlů (m je počet slu­čo­va­ných polí).

Další zlep­šo­vák spo­čívá v po­u­žití bitmapy, která ozna­čuje ne­prázdné pozice.

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

Nej­bližší ob­sa­zený slot pak najdu ně­ko­lika LZCNT in­struk­cemi (number of lea­ding zeros) zcela bez pod­mí­ně­ných skoků (po odro­lo­vání).

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

Ma­so­chisté, kteří si libují v in­strin­sics, 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 čte­náře.)

Bitmapy mají pří­jem­nou sy­ner­gii s recyk­lací uzlů, pro­tože ne­mu­sím vy­nu­lo­vat celý uzel, ale jen pří­sluš­nou bitmapu. Si­tu­aci o něco málo kom­pli­kují možné du­pli­ká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ži­tečný? Roz­hodně ne pro sor­to­vání ne­se­řa­ze­ného pole. Kdy­bych ho použil tak jak je, al­go­rit­mus by zde­ge­ne­ro­val na oby­čej­nou trie, která, tak jak jsem ji tu pre­zen­to­val, není příliš pa­mě­ťově úsporná. Na druhou stranu mi přijde jako celkem uži­tečný ná­stroj pro velice široké slu­čo­vání mnoha se­řa­ze­ných polí1 , pro­tože slo­ži­tost ne­na­růstá s jejich počtem.

Na druhou stranu je třeba stále mít na paměti oš­k­livé ta­jem­ství všech radix-sortů: Pro re­pre­zen­taci každé z n roz­díl­ných hodnot po­tře­buji log2(n) bitů. Ra­di­xová slo­ži­tost n*b tedy není nic víc než n * log(n), kde jsme si řekli, že n ne­po­roste nad prak­tické meze.

  1. Ta mohla vznik­nout oby­čej­ným radix sortem a jsou do­sta­tečně malá, že se vešla do CPU cache pro ma­xi­mální efek­ti­vitu
  2. Na líné řazení se spo­lé­hali autoři Haskellové knihovny pro hle­dání nej­bliž­ších sou­sedů vy­u­ží­va­jící cover tree.
@kaja47, kaja47@k47.cz, deadbeef.k47.cz, starší články