How to sort out your life in O(n) time

slides

Flattr this!

Posted in Uncategorized | Leave a comment

Jak řadit v lineárním čase, křísit mrtvé a dosáhnout osvícení

Když se řekne řazení, většina lidí, kteří napsali víc než jednu smyčku a dva ify se zeptá: „A je to nutné?“ Řazení nemá pověst nejrychlejší operace pod sluncem a je preferováno, když se tomu dá vyhnout např. hashováním. Někdy je to ale nutné a v tom případ je třeba zatnout zuby a připravit se na nějaký ten rozpálený křemík.

Jedinci, kterým na zdi visí portrét Donalda Knutha si okamžitě vybaví, že řazení má v nejlepším případě časovou složitost O(n log(n)) a vzpomenou si na všechny ty quicksorty, merge sorty, heap sorty a různé monstrózní hybridy.

Pravda je ale o něco komplikovanější: Řadící algoritmy založené na porovnání potřebují provést řádově n * log(n) porovnání. Existují však i jiné způsoby řazení, které jsou založené na frekvenci, rozložení a histogramech, a ty nepotřebují porovnávat nic s ničím. Stačí jim lineární čas, a co víc: Jsou zatraceně rychlé.

Jedním z nejjednodušších je counting sort, který funguje tak, že vstupní pole jednou projdu, zaznamenávám kolikrát jsem na který element narazil a na konci položky vyplivnu ve požadovaném pořadí.

Může to vypadat například takhle (psáno v jazyce fast-and-loose Scala):

def countingSort(arrayToBeSorted: Array[Int]): Array[Int] = {
  val buckets = new ArrayBuffer[Int]()

  for (el <- arrayToBeSorted) {
    buckets(el) += 1 // increment counter
  }

  val result = new ArrayBuffer[Int]()

  for {
    i <- 0 until buckets.size // element
    _ <- 0 until buckets(i)   // count
  } result += i // append element

  result.toArray
}

Bystří si jistě už všimli, že slibovaná lineární složitost má pár nedostatků. Jak je vidět, je třeba mít po ruce pole, které je velké jako největší element ve vstupní kolekci. Co se stane, když bych se pokusil seřadit čísla 1 a 1844674407370­9551616? Musel bych přikoupit hodně paměti.

Counting sort má lineární časovou složitost, ale potřebuje prostor přímo úměrný rozpětí možných hodnot. To je ideální, když řadím například pole bajtů, ale v ostatních případech to má k ideálu velice daleko.

Řešením je radix sort. Ten funguje velice podobně, ale k řazení nepoužívá celou hodnotu najednou (tj. všechny dostupné bity), ale vždy jen její část. Tak například, když bych chtěl seřadit pole čtyřbajtových intů, radix sort může v první iteraci data rozdělit podle nejdůležitějšího bajtu. Vznikne tak 256 skupin (partition), kdy každá z nich obsahuje jen elementy shodující se v nejvyšším bajtu. V další iteraci provedu stejný krok pro každou skupinu zvlášť s druhým nejdůležitějším bajtem. Vznikne tak 256×256 podskupin, z nichž každá sdílí horních 16 bitů. Když udělám 4 rekurzivní iterace a spojím vzniklé podoblasti, dostanu seřazenou sekvenci. Stačí mi k tomu jen 4 iterace vstupním polem a pokud se nepletu, 4 lze považovat za konstantu.

Časová složitost radix sortu je O(kn), kde n je délka vstupu a k je velikost jedné položky. Když řadím pole primitivních typů nebo malých structů, k je malé a všechno běží velice rychle. Když se však rozhodnu radixovat dlouhé stringy nebo velké a komplikované (nedej bože nekonečné) objekty, k začne dominovat a situace je o poznaní méně růžová. Zvlášť s přihlédnutím k faktu, že log2(232), což je maximální rozsah 32 bitového intu, je jen 32, k tedy nemá moc místa pro manévrování. Stačí když radixovaný objekt trochu nakyne a najednou n log(n) představuje správnou cestu vpřed.


Majíce tohle všechno na paměti jsem napsal vlastní implementaci radix sortu ve Scale, abych na vlastní oči viděl jak rychlé je rychlé řazení. Můj kód je založený na moudrosti tohoto článku a jde o LSD radix sort. Neřadí tedy od nejvýznamnějšího bajtu, ale od toho nejméně významného. To se může zdát neintuitivní, ale funguje to, protože radix sort je stabilní algoritmus – více detailů ve zmiňovaném článku. Pro čtyřbajtové inty potřebuje provést 5 iterací nad daty – v jedné spočítá histogram frekvencí jednotlivých bajtů a ve čtyřech zbývajících přesouvá položky podle histogramů.

To však není jediný způsob. Existuje spousta dalších odrůd radix sortu: je možné řadit od nejdůležitějších bitů, od těch nejméně důležitých, může být in-place, může potřebovat extra O(n) místa, může být hybridní (použije se jiný algoritmus, když je jasné, že se nevyplatí radixovat), může být paralelní, může brát v potaz cache, může používat software managed buffers, atd. atd.

Výsledek? Jak řekl klasik: ho-ly-shit. Přesvědčte se sami:

Graf ukazuje rychlost řazení polí obsahujících náhodné inty mým radix sortem a funkcí java.util.Arra­ys.sort, což je hybrid mezi dual-pivot quicksortem, mergesortem a insertion sortem. Jak je vidět, radix sort běží v lineárním čase a navíc je velice rychlý. A to jde o neaktuální graf, který měří rychlost verze bez několika optimalizací3, a navíc do průměru započítává i počáteční pomalé běhy, kdy JVM nebyla zahřátá. Malé zpomalení s rostoucí velikostí dat není způsobeno asymptotami algoritmu, ale hardwarovými efekty. Všechny operace jsou rychlé, čtení je vždycky sekvenční, zápis je skoro sekvenční do 256 míst v paměti, které se pohodlně vejdou do cache5, neobsahuje žádné nepředvídatelné skoky, které by vedly k brach-miss – situace je zcela odlišná od quicksortu, který neobsahuje nic jiného než nepředvídatelné skoky.1

1GB dat (250M zcela náhodných intů) je na mém stroji možné seřadit za 4 vteřiny, a 1GB floatů za 2.9 vteřiny2. To odpovídá 14.9 ns, resp 10.8 ns na každý element pole. Když se celé pole vejde do cache, řazení je o něco rychlejší. Každému intu i floatu stačí 9.7 ns, což odpovídá rychlosti řazení přes 100M/s v jednom vlákně. Radix sort navíc umí poznat, když není třeba řadit podle určitého bajtu a přeskočí celý jeden řadící krok, což vede ke značnému zrychlení. Například když řadím pole obsahující jen čísla mezi 0 a 64k, dvě iterace z pěti jsou přeskočeny a je třeba jen 8 ns místo 14.9 ns na element.

Řazení floatů radix sortem je poměrně delikátní záležitost, která se spoléhá na to, že čísla mají bitovou reprezentaci podle specifikace IEEE754. Nejjednodušší možností je data seřadit jako integery a pak opravit pořadí záporných hodnot. O něco rafinovanější je floaty na vstupu transformovat, tak aby by je bylo možné řadit jako 32-bitová čísla bez znaménka a na konci je změnit zpátky na float. Nejlepší řešení je však tomu předejít a upravit offsety záporných radixů.


Ok, takže když je radix sort tak rychlý, k čemu se dá použít? K řazení, pochopitelně. Ale najde i uplatnění tam, kde se řazení zdá, jako jít s kanónem na vrabce.

Jedno takové místo je seskupování.

Když bych chtěl seskupit moře hodnot podle jejich klíče z páru klíč-hodnota (kde klíč i hodnota jsou čtyřbajtové inty zabalené do osmibajtového longu), mohl bych to s do očí bijící neefektivitou napsat v jednom řádku Scaly.

keyvals.groupBy(getKey).map { case (k, vs) => vs.map(getKey) }

To je pěkné, ale nešlo by to rychleji, nejlépe s minimem alokací? Šlo. S pomocí mapy, která obsahuje měnitelné buffery.

val buffers = mutable.Map[Int, mutable.ArrayBuilder.ofInt]()

for (kv <- keyvals) {
  val k = getKey(kv)
  val v = getVal(kv)
  buffers.getOrElseUpdate(k, new mutable.ArrayBuilder.ofInt) += v
}

buffers.toSeq.sortBy(_._1).values

Když zanedbám režii mapy (která může být zanedbatelná, ale také nemusí), za pozornost stojí poslední řádek: sortBy. Najednou jsem tam, kde jsem začal. Když je velikost každé skupiny malá, je jejich počet úměrný velikosti vstupních dat a konečné řazení bude velice drané.

Když vím, že klíče náleží do určitého rozpětí, můžou použít pole bufferů.

val buffers = new Array[mutable.ArrayBuilder.ofInt](length)

for (kv <- keyvals) {
  val k = getKey(kv)
  val v = getVal(kv)
  if (buffers(k) == null) {
    buffers(k) = new mutable.ArrayBuilder.ofInt
  }

  buffers(k) += v
}

buffers.filter(_ != null).map(_.result)

Jak je vidět, na konci není třeba nic řadit, protože výsledek je přirozeně seřazen podle klíče ve stylu counting sortu. Pokud mě trápí opakované alokace a kopírování dat v rostoucích bufferech a můžu si dovolit dva průchody daty, je možné v jenom průchodu spočítat velikosti jednotlivých skupin, alokovat pro ně pole velká tak jak je potřeba, a v druhém průchodu je naplnit. To může být v některých případech rychlejší a v jiných, kdy nechci přijít s OutOfMemoryError po funuse, nezbytné.

Tohle je lepší, ale nešlo by to… však víte… rychleji?

Tady do hry vstupuje radix sort. Co kdybych pole seřadil podle prvních 4 bajtů (tj. té části, která obsahuje klíč) a pak výsledek v jedné iteraci rozřezal na kousíčky. Nějak takhle:

val scratch = new Array[Long](keyvals.length)
val (sorted, _) = RadixSort.sort(keyvals, scratch, 0, keyvals.length, 4, 8, false)

var i = 0
(0 until numberOfKeys iterator) map { key =>
  var start = i

  // find end of the current group
  while (i < keyvals.length && getKey(sorted(i)) == key) { i += 1 }
  val end = i

  // empty group
  if (start == end) {
    null

  // copy values  into a separate array
  } else {
    val res = new Array[Int](end - start)
    var j = 0
    while (start < end) {
      res(j) = getVal(sorted(start))
      start += 1
      j += 1
    }

    (key, res)
  }
} filter (_ != null )

Na rozdíl od předchozích případů je kód o něco složitější, na druhou stranu nejde o pseudokód, ale skutečně fungující fragment. A co víc: Je skutečně rychlejší než všechny předchozí ukázky.

Proč?

Jedno vysvětlení je nasnadě: Pokud je skupin hodně a data jsou více méně náhodná, v předchozím případě každý nový klíč může vést ke několika drahým cache-miss. Jeden načte více méně náhodný pointer na buffer a druhý načte objekt bufferu a třetí konec pole v bufferu na které se bude vkládat. A už jsem zmiňoval, že jít do DRAM je drahé? Myslím, že ano. Radix sort tyhle problémy nemá, párkrát projede data, seřadí je a výsledek se nakonec ještě jednou proletí a rozseká na kousíčky. A všechno udělá velice rychle.


Dále k tématu:


  1. Špatně odhadnuté skoky představují tak velký problém, že někdy je quicksort rychlejší, pokud je záměrně zvolen špatný pivot. Když pivot data rozdělí fifty-fifty, je skok nepředvídatelný. Když je ale rozdělí třeba 25–75, procesor bude odhadovat, že skok půjde doprava a většinou se trefí. viz: An Experimental Study of Sorting and Branch Prediction. V tomto kontextu za pozornost také stojí: Branch Prediction and the Performance of Interpreters – Don't Trust Folklore a Branchless code sequences.
  2. Předpokládám, že v případě floatů JVM JIT zkompiluje kód s využitím SIMD instrukcí a proto je výsledek rychlejší. Nebo také může jít o korelace mezi bitovými vzory v reprezentaci floatů. Podrobnosti jsem nezkoumal.
  3. Jedna z optimalizací, která může za velice specifických situací vést k dvojnásobnému zrychlení je zpětná iterace jedné smyčky. V běžných situacích, kdy jsou vstupní data větší než cache přidá k dobru ~10%. Když je vstupní pole velké, nemá žádný vliv. Další optimalizaci, kterou jsem zkoušel, byly software managed buffers zmiňované zde a zde. Ty fungují tak, že data nezapisuji přímo do paměti, ale do malého bufferu, který se vejde do cache a do paměti ho zkopíruji jen, když se naplní. Tohle přidalo pár procent výkonu pro řazení intů, ale vedlo to k drastickému zpomalení v případě floatů. Výsledek nebyl nijak přesvědčivý a proto jsem tuhle optimalizaci nakonec nepoužil.
  4. To s přihládnutím k tomu, že radix sort musí udělat 5 průchodů celým polem znamená 2.2ns/7taktů CPU na element v každé iteraci. Z paměti je přitom třeba číst/zapisovat 4.8 GB/s (vysvětlení v tomto paperu, který hlásí, že jeden průchod jejich radix sortu, který dává 240M/s na čtyřech jádrech, trvá 13.2 taktu na element – pozor jejich čísla nesedí).
  5. V literatuře se uvádí, že radix sort nedělá sekvenční zápis v tom smyslu, že nezapisuje do jednoho místa, ale do mnoha, v tomto případě 256. To vede k neoptimálnímu využití propustnosti pamětí a může vést k neustálým výpadkům TLB. Moje experimenty na moderním hardwaru však nejsou přesvědčivé aspoň, co se týká TLB. Perf hlásí skoro nulu dTLB-load-misses a dTLB-store-misses, ale na druhou stranu cache miss se stále nějaké vyskytují (jako load-miss tak i store-miss).
  6. Se složitostí radix sortu je to o něco komplikovanější. Je pravda, že běží v čase O(kn), ale povaha onoho k je trochu matoucí. Když bych měl n různých čísel, nejmenší počet bitů potřebných k reprezentaci každého z nich je log<sub>2</sub>(n). Ono O(kn) je tedy vlastně jen O(n log(n)) v přestrojení.

Flattr this!

Posted in Algo, CPU, Paměť | 2 Comments

Dualismus hardwaru a softwaru, strojů a virtuálních strojů

René Descartes věřil v myšlenku dualismu, který se dá diletantsky vysvětlit tak, že tělo a duše jsou dvě různé vzájemně neslučitelné kategorie.

V podobném duchu se nese jiný dualismus, který je na rozdíl od Descarterových tvrzení méně metafyzický, ale zato je přítomný a hmatatelný. Oba pak mají jedno společné: Jde o falešnou dichotomii i když všichni v hloubi duši toužíme po opaku.

Pochopitelně mluvím o dualismu hardware a software, strojů a virtuálních strojů, kdy jeden dělá to, co to druhý nemůže, kdy „myslící“ program ovládá „hmotu“ stroje. Tato demarkace na první pohled působí správně a uspokojivě, jako kdyby to tak skutečně bylo. Když ale začnu problematiku zkoumat podrobněji a vezmu to oklikou přes historii CPU architektur, toto rozdělení, před chvílí ještě krystalicky jasné, najednou začíná mizet v mlhách.

Všechno, co vypadá jako neměnné pravdy vytesané do kamene, jsou jen rozhraní a úrovně abstrakce.

Už instrukční sada procesoru (ISA) není ten pomyslný pevný bod, na kterém stojí všechny ostatní želvy, ale také jen abstrakce. Téměř žádný dnešní procesor (možná s výjimkou některých jednoduchých RICSů nebo DSP) nevykonává instrukce tak, jak jsou za sebou vyskládány v binárce, ale dynamicky je v hardwaru překládá na jednodušší (RISCovější, i když tohle označení už moc neznamená) mikro operace (μops). Jedna instrukce se může rozpadnout na několik μops nebo naopak může být běžná sekvence několika instrukcí přeložena na jednu mikro operaci (tzv macro-op fusion)1. Situaci dále komplikuje mikrokód, kdy jsou některé instrukce procesorem přeloženy na malý program (například trigonometrické funkce jsou mikrokódovány, práce se denormálními floaty nebo GATHER instrukce na Haswellech, vylepšená v následujících generacích)14. Instrukce tedy není nutně nejmenší jednotka práce2 a často ani není přímo implementována v hardware.

Tato organizace není použita zbůhdarma, ale má dobrý důvod. Není totiž nutné plýtvat křemíkem na každou funkcionalitu, která je v programech používána jen okrajově. Místo toho je možné použít mikroarchitekturu s menším počet vysoce optimalizovaných funkčních jednotek, které jsou sdílené mezi instrukcemi3. Dále to zjednodušuje návrh pipeline, protože není třeba počítat s operacemi, které mají divoce rozdílné latence, a vede to k efektivnějším out-of-order jádrům, které můžou vykonávat různé μop (které trvají ±stejný počet taktů) různých CISC instrukcí najednou a využít tak dostupný ILP. Jako příklad můžu uvést stařičké VAXy. Ty například přímo podporovaly komplikovanou instrukci INDEX, která byla používána jen velice zřídka, a proto nikdo nestrávil čas optimalizací její hardwarové implementace. To došlo tak daleko, že když někdo napal stejnou funkcionalitu pomocí jednodušších operací, softwarový výsledek byl rychlejší než přímá implementace v křemíku. Víc o tom mluví Onur Mutlu v jedné ze svých skvělých přednášek o architektuře procesorů.

Jestliže instrukce jsou atomy, pak mikro operace jsou jejich kvarky.

Jako další příklad na spektru abstrakce může posloužit Transmeta. Ta v dobách své nedlouhé slávy produkovala čipy, které dokázaly rozběhat barokní x86 binárky. Sám hardware přitom neuměl nic z monstrózní x86 ISA, ale místo toho běhal na vlastní proprietární VLIW architektuře a x86 program byl za běhu překládán pomocí takzvaného Code Morphing software (CMS) do interní ISA4. Nepřekládal ale všechno. Zpočátku kód jen naivně interpretoval, když pak objevil opakovaně prováděný blok, rychle ho překompiloval z x86 do nativního VLIW formátu, když pak zjistil, že jde o skutečný hotspot, provedl důkladnou (a pomalou) kompilaci, jejímž výsledkem byl kvalitní a rychlý kód. Procesory od Transmety se tak v mnohém podobají JVM, jen s tím rozdílem, že jejich vstupem není Java bytekód, ale x86 binárka, cíl JITu není x86, ale interní nízkoúrovňový VLIW formát a hardware je specificky navržen pro tento účel.

VLIW instrukce se v mnohém podobají mikro operacím – jde o operace na velice nízké úrovni, které věrně opisují chování a mikroarchitekturu hardware5. Jsou prováděny ve statickém pořadí určeném kompilátorem (statically scheduled) a obnažují všechny interní detaily procesoru, jako jsou počty funkčních jednotek, délku pipeline a latence jednotlivých operací. To vede k problémům s přenositelností, kdy program zkompilovaný pro jeden model VLIW stroje nefunguje na novějším modelu, který má víc nebo rychlejší funkční jednotky.9 Když je však tento interní formát skryt pod jednou úrovní abstrakce (jako v případě x86 a dalších ne-VLIW ISA), může se mikroarchitektura libovolně měnit a vrstva o jednu úroveň výše se postará, aby byly všechny prostředky co nejlépe využity – nezáleží jestli jde o celkem naivní překlad v hardware, nebo dynamickou rekompilaci v software. Dokonce je teoreticky možné takto postavit stroj, který zvládá několik různých ISA najednou (třeba x86 a ARM).

V případě Transmety se interní formát měnil celkem znatelně – první generace procesorů pojmenovaná Crusoe byla široká 128 bitů a každá instrukce obsahovala 4 operace. Další generace byla dvakrát širší s dvakrát větším počtem operací v instrukci. Více informací je v oficiálním whitepaperu nebo ve dvou článcích o snahách reverzně vyinženýrovat mikroarchitekturu Transmety.

Velice podobně jako čipy už dlouho mrtvé Transmety pracují některé ARM čipy z rodiny Tegra od Nvidie založené na mikroarchitektuře Denver (nebo ruský Elbrus, když už jsem u toho). Na rozdíl od Transmety, která začala program v interpretru, Denver umí přímo vykonávat ARM kód na celkem základní úrovni. Hlavní síla čipů však spočívá v proprietární interní VLIW ISA, která je 7-wide (7 ops v instrukci). Program je zpočátku vykonáván v ARM módu s tím, že jakmile software objeví horký kus kódu, začne kompilovat. Rozpoznávání hotspotů nefunguje na úrovni základních bloků, ale na úrovni průchodů (trace), které jsou typicky tvořeny jednou iterací smyčky a můžou procházet mnoho základních bloků a dokonce i volání funkcí. Když je nalezena frekventovaná trace ARM instrukcí, je softwarově překompilována do VLIW ISA jako lineární sekvence operací, bez větvení a bez odboček. V mnohém tedy připomíná trasovací (tracing) JIT kompilátory, jako je třeba ten v PyPy (používaný i v jiných jazycích než Python), jen netrasuje nějakou formu bytekódu nebo operace interpretru (jako v případě meta-tracing JIT), ale přímo ARM instrukce.

Jako další příklad může sloužit společnost Azul, která dříve vyráběla čipy Vega specializované pro běh Java aplikací. Interně šlo o velice jednoduché a nepříliš rychlé in-order RISC procesory, které byly navrženy s ohledem na specifikaci Javy. Nepokoušely se však přímo vykonávat Javovský bytekód – lidé se o to dříve pokoušeli a zjistili, že je to velice špatný nápad. Na rozdíl od Teger a čipů Transmety měly Vegy vlastní operační systém, standardní knihovny a JVM zkompilované do nativního kódu, ale rozhraním pro uživatelské aplikace byl jedině Java bytekód. Šlo tedy o poměrně vysokou úroveň abstrakce, která návrhářům10 dala velkou míru volnosti jak společně navrhovat michroarchitekturu a JIT. Hardware se mohl mezi generacemi měnit (a také se měnil), aniž by to narušilo kompatibilitu pro uživatelské aplikace.

Úroveň abstrakce může být však i jinde, o čemž svědčí chystaný procesor Mill. Jde opět o VLIW staticaly scheduled, open pipeline stroj, který šířkou, kdy v jedné instrukci může být až 33 operací, překonává všechno ostatní. ISA je zamýšlena jako cíl pro kompilátor a bude přímo vykonána procesorem na nejnižší možné úrovni bez dalších vrstev abstrakce. Od ostatních VLIWů se Mill liší v tom, že není navrhován jako jedna architektura, ale rodina architektur s různými schopnostmi a hardwarovými prostředky.12

VLIW má problém v tom, že jde o příliš nízkou úroveň abstrakce – v podstatě představuje obnaženou mikroarchitekturu. Když se hardware změní, software se musí přizpůsobit. Mill se tohle snaží vyřešit tím, že programy nebudou distribuovány jako binárky přímo určené ke spuštění, ale budou na něco abstraktnější úrovni s tím, že se před spuštěním nebo během instalace staticky specializují pro konkrétní hardware. Na rozdíl od příkladů v minulých odstavcích by nemělo jít o nijak komplikovanou rekompilaci, ale celkem přímé mapování požadavků programu na prostředky hardware. Mill tak abstrahuje nejen možnosti daného stroje, ale i binární kódování instrukcí, které může být specifické pro každý model a přesto si zachovává výhody statického schedule, který nepotřebuje komplikovaný dekodér nebo out-of-order hardware.6

O Millu zatím nemůžu prohlašovat nic konkrétního, protože se ještě nedá koupit reálný křemík. Všechny informace pocházejí ze série velice detailních přednášek, které uspořádal Ivan Godard a ukazuje v nich všechny vlastnosti, kterými se Mill liší od běžných strojů a architektur.

Z předešlých odstavců lze vypozorovat jedno společné téma: všechno je jen abstrakce na určité úrovni. Od velice nízké jako v případě Millu, přes o něco vyšší (x86 a ARM ISA) až po relativně vysokou jako je Java bytekód. A to nemluvím o FPGA, které všechno tohle staví na hlavu, protože dovolí přímo programovat hardware. A víte jak Intel nedávno koupil Alteru? To znamená jediné: dočkáme se procesorů s integrovaný­mi FPGA.


Teď když jsem tu napsal přes 1000 slov o tom, že hranice mezi hardwarem a softwarem je mlhavější, než se na první pohled může zdát, nabízí se jedna zřejmá otázka: Neexistují nějaké principy psaní programů, které nějakým zásadním způsobem benefitují hardwarusoftwaru?

Zcela nepřekvapivě: Ano.

Jedna vlastnost, ze které profitují všichni, je předvídatelnost. Pokud se strojový kód chová předvídatelně, může hardware rozpoznat vzory jeho chování a pomoci. Sázka na předvídatelnost vedla k obohacení procesorů o branch prediction jednotky, které se snaží uhodnout zdali program skočí nebo propadne7, a branch target buffery, který předvídá cíle skoků na registr (typické při implementaci virtuálních metod), cache a prefetching.

Předvídatelnost pomáhá softwaru v tom, že dynamická prostředí, jako je JIT kompilátor v JVM, může specializovat kód pro běžné případy a ty, které nastanou jen zřídka nebo vůbec, ignorovat. Takto může devirtualizovat a inlinovat virtuální metody, což otevře dveře celé plejádě dalších optimalizací. Pro dynamické jazyky je předvídatelnost ještě důležitější, protože pokud má program rozumně statickou strukturu, virtuální stroj může objevit skryté třídy a typy argumentů funkcí a specializovat pro ně kód8. Předvídatelnost je také benefitem pro tracing JIT, které nahrávají (trasují) lineární průchod kódem, který je často prováděn, a čím méně odskoků a odboček v něm je, tím jednodušší má práci13. Stejně tak práce Maxime Chevalier-BoisvertBasic Block Versioning, ocení předvídatelnost tím, že bude potřebovat vygenerovat menší množství verzí základních bloků. Z předvídatelnosti také těží přístup založený na specializaci AST (jako je Truffle+Graal backend pro Jruby), protože (očividně) není třeba provádět takové množství specializací, což vede k menším a rychlejším několikaúrovňovým polymorfním inline cache.

Další dobrou vlastností pro všechny zúčastněné je malý a kompaktní kód. Na úrovni hardwaru je to proto, že se kód vejde do cache. Nejde jen o ty hlavní úrovně jako L3/L3/L1I, ale také L0, která na některých strojích obsahuje dekódované mikro operace. To pomůže v případech, kdy hlavním úzkým hrdlem je dekodér instrukcí. Pokud je však smyčka velice těsná (něco jako pár desítek μops na Haswell/Broad­well/Skylake x86) a celá se vejde do loop bufferu, hardware rozpozná, že si vystačí s tím, co je v loop bufferu a vůbec nebude komunikovat s cache a dekodérem. To může vést ke zrychlení programu a úsporám energie. Navíc nepředvídatelnému a nerozumně rozvětvenému kódu hrozí, že skočí na místo, které ještě není v cache a bude muset čekat.11

Virtuální stroje kompilující kód po metodách (jako typické JVM) preferují malé metody, protože jim to dává volnost v rozhodování, zdali danou metodu inlinovat nebo ne. Inlinování zvětšuje výsledný kód a proto, kdyby to JIT začal dělat příliš velkoryse, kód by nakynul a nemusel by se vejít do cache, což by vedlo ke znatelnému zpomalení. Inlinování je důležité, ale jen v rozumné míře. Když se to přežene, může to uškodit. Z toho důvodu JIT typicky vybírá kandidáty pro inlinování heuristicky a jedna z nich je velikost metody s tím, že velké metody mají jen malou šanci být inlinovány. Například JVM miluje malé metody a když se rozjede, začne je spekulativně inlinovat několik úrovní hluboko. Malé metody zároveň vedou k organizaci kódu s malými funkčními celky, které mají jasně definovanou roli, kdy každá funkce dělá jednu věc, což je obecně dobrá věc.


Nakonec to vypadá, že všechno jsou to jen virtuální stroje stojící na ramenou jiných virtuálních strojů. Dokonce i samotné Céčko, kterému se někdy říká přenositelný assembler definuje chování virtuálního stroje, který není nikde implementovaný v hardware, ale shodou náhod se na skutečný hardware dá docela pěkně napasovat. Všechny virtuální stroje se nakonec od sebe liší jen tloušťkou abstrakce, kterou je třeba překonat, aby program dokázal rozpohybovat reálné elektrony na reálném křemíku.


Dále k tématu

  1. Například poslední Intelí x86 kusy překládají sekvence některých instrukcí, jako cmp nebo aritmetické operace, následovaných skokem na jeden μop. Je tedy možné, že všechny instrukce, které řídí smyčku (inc, cmp, jXX) skončí jako jediný μop přímo implementovaný v hardware. Tohle je podle mě dobrý způsob jak utratit dostupný křemík, protože těsné smyčky jsou časté a program v nich stráví hodně času.
  2. I když některé jednoduché instrukce jsou přímo implementované v hardware. Jak je vidět z některých mikroarchitek­totických schémat, x86 CPU mají několik dekodérů z nichž většina dokáže dekódovat jen jednoduché instrukce, které jsou přeloženy na jeden μop a jeden, který překládá komplikovanější operace nebo pracuje s mikrokódem. Dekódování může pro komplikované CISC ISA představovat úzké hrdlo.
  3. Například aritmetické instrukce, které mají cíl nebe zdroj v paměti jsou přeloženy na jeden μop, který provede load/store a je neplánovaný na některou z dostupných load/store jednotek a další μop samotné aritmetické operace. Pro programátora viditelná architektura tak může být register-memory, ale mikroarchitektura je typu load-store.
  4. CMS naběhl ještě před startem operačního systému a byl to jediný kus softwaru napsaný přímo pro interní VLIW architekturu, všechno ostatní byl tradiční x86 od operačního systému až po uživatelské aplikace.
  5. Jenom s tím rozdílem, že VLIW přesně pasuje na dostupné funkční jednotky a celá dlouhá instrukce je vykonána paralelně, zatímco každý μop je na out-of-order strojích plánován a prováděn dynamicky. VLIW zachycuje statický instrukční paralelismus, μop na OOO dynamický.
  6. Mill se v mnohém podobá Itaniu, které částečně zapadá do rodiny VLIW strojů – každá instrukce je „dlouhá“ a obsahuje několik operací (v tomto případě se skupině operací říká bundle) a explicitní označení datových závislostí mezi operacemi. Pokud je mikroarchitektura dané implementace Itania dostatečně široká, stroj může vykonat pár bundlů paralelně a nemusí sám dynamicky zjišťovat závislosti mezi operacemi. Záměrem bylo dosáhnout většího ILP na novějších strojích bez nutnosti out-of-order hardwaru a se zachováním binární kompatibility.
  7. V současnosti branch predictor dosahuje úctyhodné přesnosti. Z paperu Prediction and the Performance of Interpreters je vidět, že procesor dokáže předpovídat skoky ve smyčce interpretru a svým způsobem pracuje na meta úrovni – nepředvídá skoky v kódu, ale v kódu, který interpretuje kód.
  8. Specializace se dá dotáhnout dál než skryté třídy/specializace kódu. Paper Adaptive Just-in-time Value Class Optimization se zabývá specializací kompozitních datových struktur.
  9. I některé ne-VLIW ISA odhalují interní detaily mikroarchitektury. Jde například o delay-slot.
  10. Mezi mě mimo jiné patřil Cliff Click, který měl o procesorech Vega skvělou přednášku, která kdysi, spolu s Performance Anxiety od Joshe Blocha, nasměrovala můj zájem směrem k hardwaru a procesorům.
  11. Ukazuje se, že některým serverovým úlohám nestačí instrukčí cache.
  12. Ale ani toto není nový nápad. Sám Ivan Godard říká, že Mill je rodina architektur ve stejném smyslu jako IBM System/360.
  13. Tady stojí za zmínku trace cache Pentia 4, která pracuje v duchu tracing JITu a ukládá lineární sekvenci dekódovaných mikro-operací.
  14. Podle knihy The Soul of a New Machine Tracy Kiddera stroje v sedmdesátých letech přecházely od přímé implementace instrukcí v hardware na mikrokód.

Flattr this!

Posted in CPU, VM | Leave a comment

Někdy je nejchytřejší nedělat nic chytrého (další kapitola nekonečného příběhu o optimalizaci)

Nedávno se na Twitteru objevil článek Fast Nearest Neighbor Queries in Haskell, který představuje novou knihovnu, která dokáže za pomoci velice chytrých algoritmů dělat nearest neighbor dotazy rychleji než jakýkoli jiný nástroj na planetě.

Nearest neighbor dotaz je takový, který má najít nejbližší bod z dané množiny bodů v určitém metrickém prostoru. Naivně se dá spočítat tak, že jeden po druhém projde všechny body s tím, že si pamatuje ten nejbližší a jeho vzdálenost1. Chytré řešení všechny body zaindexuje do stromové datové struktury, která nějakým rafinovaným způsobem rekurzivně rozdělí prostor do oblastí, a během hledání většinu těchto regionů přeskočí, protože je garantováno, že v nich obsažené body nemohou být blíže než současný kandidát. Zmíněná knihovna k tomuto účelu používá cover tree.

Článek rozhodně stojí za přečtení, už jen proto, že popisuje, jak psát rychlý numerický kód v Haskellu. Zcela nepřekvapivě jde hlavně o jednu věc: efektivní využití CPU cache ať už pomocí inlinování datových struktur, aby nebylo třeba nahánět zdivočelé pointery, nebo cache oblivious datových struktur jako třeba Van Emde Boas. Výsledky těchto snah jsou shrnuty v několika grafech s výsledky benchmarků.

A právě tady do příběhu vstupuji já.

Když jsem poprvé zíral na výsledná čísla, něco mi v nich nehrálo. Něco mi na nich nesedělo. Jeden benchmark pracoval se 100000 body o 384 dimenzích, ke každému hledal jednoho nejbližšího souseda měřeného euklidovskou vzdáleností a běžel 13.6 minuty. To by na mém stroji odpovídalo víc jak 2600 miliardám taktů procesoru. Naivní řešení by muselo projít všech 100000 bodů a každý porovnat se 100000 ostatními body a při každém porovnání pracovat s 384 prvky vektoru. Dohromady by bylo třeba 3840 miliard operací.

Skoro to vypadalo, že naivní přístup by nemusel být o moc pomalejší.

Bruteforcing for fun and profit

Sedl jsem si tedy ke stroji a za chvíli vyšvihl jednoduchý program v Céčku, který hledá nejbližší sousedy nejnaivnější možnou metodou.

double euclidean_distance_sq(double *vecs, int veclen, int a, int b) {
  double d = 0;
  for (int i = 0; i < veclen; i++) {
    double t = vecs[veclen * a + i] - vecs[veclen * b + i];
    d += t * t;
  }
  return d;
}

int nn(double *vecs, int veccnt, int veclen, int idx) {
  int nearest = -1;
  double dist = INFINITY;
  for (int i = 0; i < veccnt; i++) {
    double d = euclidean_distance_sq_vec_fl(vecs, veclen, idx, i);
    if (i != idx && d < dist) {
      nearest = i;
      dist = d;
    }
  }
  return nearest;
}

struct idxdist {
  int idx; // -1
  double dist; // INFINITY
};

int veccnt = atoi(argv[1]); // počet vektorů
int veclen = atoi(argv[2]); // délka jednoho vektoru

double *vecs = loadVectorsOrGenerateThemOrWhatever(veccnt, veclen);
struct idxdist *idxdists = initializeIdxsToMinusOneAndDistToMinusInfinity(veccnt);

for (int i = 0; i < veccnt; i++) {
  idxdists[i].idx = nn(vecs, veccnt, veclen, i);
}

Brnkačka. Program zkompilovaný pomocí gcc -O3 -funroll-loops běžel na 10× menších datech 36 vteřin. Aby byl srovnatelný s chytrým algoritmem, musel by běžet něco mezi 8 a 10 vteřinami.

Bylo na čase přitlačit na pilu, na vektorizovanou pilu. Nepodařilo se mi donutit GCC, aby automaticky vektorizoval smyčky a generoval SIMD instrukce7, a proto jsem se musel uchýlit k drastickým opatřením. Našel jsem seznam intrinsics – speciálních funkcí, které GCC přeloží na odpovídající instrukce9 – a ručně jsem vektorizoval funkci, která počítá euklidovskou vzdálenost. Nová verze provádí operace se čtyřmi doubly najednou2,8 a GCC tuto smyčku ještě odroluje, což přidá dalších pár procent výkonu6.

double euclidean_distance_sq_vec(double *vecs, int veclen, int a, int b) {
  __m256d d = _mm256_set1_pd(0.0);
  for (int i = 0; i < veclen; i+=4) {
    __m256d aa = _mm256_load_pd(vecs + veclen * a + i);
    __m256d bb = _mm256_load_pd(vecs + veclen * b + i);
    __m256d t = _mm256_sub_pd(aa, bb);
    d = _mm256_add_pd(d, _mm256_mul_pd(t, t));
  }

  __m256d sum = _mm256_hadd_pd(d, d);
  return ((double*)&sum)[0] + ((double*)&sum)[2];
}

Výsledek není na pohled nijak pěkný, ale mě nešlo o krásu, ale jen a pouze o surovou rychlost. A ta se zvýšila: Vektorizované verzi na výpočet stačí jen 20 vteřin. To je lepší, ale zdaleka ne ideální.

Jako další krok jsem zlepšil využití CPU cache. I když jsou vektory uložené v paměti jeden za druhým a procesor maskuje latenci pamětí načítáním dat dopředu (prefetch), limitujícím faktorem se stala propustnost paměti. Kód počítá vzdálenost jednoho vektoru se všemi ostatními. Onen jeden vektor si drží v cache, a všechny ostatní streamuje z DRAM. Každého se dotkne jednou a pak se hned přesune na další. Když začne porovnávat všechno s dalším vektorem, první vektory jsou už dávno vyhozené z cache a je třeba je opět streamovat z DRAM. To vede k mizernému využití cache.

Řešením je rozdělit pole vektorů do malých bloků, které se vejdou do cache, a počítat podobnost vždycky v jednom z těchto bloků. Když je velikost bloku 64, můžu porovnat 64*64 kombinací, ale z paměti stačí načíst jen 64+64 vektorů.3

int tile = 64;

for (int ti = 0; ti < veccnt; ti += tile) {
  for (int tj = 0; tj < veccnt; tj += tile) {
    int maxi = MIN(ti+tile, veccnt);
    for (int i = ti; i < maxi; i++) { // current vector index
      int maxj = MIN(tj+tile, veccnt);
      for (int j = tj; j < maxj; j++) {

        float d = euclidean_distance_sq_vec_fl(vecs, veclen, i, j);
        if (j != i && d < idxdists[i].dist) {
          idxdists[i].idx = j;
          idxdists[i].dist = d;
        }

      }
    }
  }
}

Kód obsahuje 4 vnořené smyčky, dvě vnější iterují po blocích (tile) a dvě vnitřní iterují obsahem bloků. Tohle uspořádání nevypadá na pohled příliš pěkně, ale plní svůj účel: program běží už jen 8.4 vteřiny a je tedy stejně rychlý jako chytrý algoritmus.

Ale to zdaleka nebylo všechno. V ten okamžik jsem se teprve dostal do ráže a v rukávu jsem měl připraveny ještě dva triky.

Dalším krokem bylo použití float místo double. Float potřebuje 4 bajty namísto osmi, z paměti tedy procesor dokáže protlačit 2× více dat, má 2× lepší využití cache a vektorové instrukce pracují s 8 floaty namísto 4 doublů.

float euclidean_distance_sq_vec_fl(float *vecs, int veclen, int a, int b) {
  __m256 d = _mm256_set1_ps(0.0);
  for (int i = 0; i < veclen; i+=8) {
    __m256 aa = _mm256_load_ps(vecs + veclen * a + i);
    __m256 bb = _mm256_load_ps(vecs + veclen * b + i);
    __m256 t = _mm256_sub_ps(aa, bb);
    d = _mm256_add_ps(d, _mm256_mul_ps(t, t));
  }

  d = _mm256_hadd_ps(d, d);
  d = _mm256_hadd_ps(d, d);
  return ((float*)&d)[0] + ((float*)&d)[4];
}

Tahle funkce je podobná té minulé, jen o něco málo ošklivější a 2× rychlejší. Program s těmito změnami běží jen 4.2 vteřiny – 2× rychleji než chytrý algoritmus.

Všechny tyto změny a úpravy přinesly skoro desetinásobnou akceleraci programu, který sám o sobě nebyl nijak strašně neefektivní. Ale samy o sobě jen optimalizovaly jedno-vláknový program a nesnažily se o nějakou razantní změnu. Poslední logickou metou byla paralelizace. Jak se ukázalo, tato změna byla nejjednodušší a nejrychlejší – stačilo před vnější smyčku přidat

#pragma omp parallel for

a zkompilovat s přepínačem -fopenmp4. Ze smyčky se rázem stala paralelní smyčka, která dokáže do posledního vytížit všechna jádra v systému. Program běžel 1.1 vteřinu. To je podle mě docela dobrý výsledek.

Pro zajímavost paralelní verze běží na mém netbooku s Atomem (2 jádra + hyperthreading, 1.66GHz) 17 vteřin, ±1000 vteřin na první generaci RaspberryPi (1 jádro, 700Mhz) a 51 vteřin na druhé generaci RasberryPi (4 jádra, 900MHz, NEON SIMD). Atom bohužel nepodporuje AVX a tak bylo třeba použít SSSE3, které pracuje jen se 128-bitovými vektory.

Shrnuto do jedné tabulky:

základní verze 36s
vektorizace 20s
lepší využití cache 8.4s
float namísto double 4.2s
paralelizace, OpenMP 1.1s

Kompletní program je k dispozici v tomto gistu.

Slovo na závěr

Na závěr musím dodat pár věcí, které zchladí všechny, kteří už-už vyrážejí slavit do ulic.

Autoři testovali onen chytrý algoritmus na Xeonu taktovaném na 2.8 GHz, takže jejich výsledky jsou o něco málo lepší, než se můžou zdát. Z prezentovaných čísel a grafů se dají těžko vyčíst nějaké absolutní hodnoty, ale vypadá to, že do testů zahrnují jak konstrukci cover tree, tak i jeho dotazování. Z jednoho grafu je vidět, že konstrukce nedominuje celkovému času, ale zabere něco jako 20%-30%. Z jiného se pak dá odvodit, že 50–90% všech přístupů do paměti skončilo jako cache-miss a kolem poloviny všech taktů procesor čekal na paměť a neměl co na práci, zatímco můj naivní program každý takt vykonal 2 instrukce a mnoho z nich pracovalo s osmi floaty najednou. Benchmarky jsem (věrný ideálům diletantství) nesnažil spustit na vlastním stroji, takže je možné, že mi mohlo uniknout něco zásadního.

Navíc: I kdyby byl onen cover tree v některých případech pomalejší, můžou existovat situace, kde bude kralovat5 – jako třeba menší počet dimenzí, jiná metrika nebo větší počet datových bodů. K tomu má stále výhodu parametričnosti, kdy můžu snadno změnit datové typy, reprezentaci a funkci vzdálenosti a nepotřebuje hledat několik bodů najednou, aby využil výhody rozdělení do bloků.

Závěr je následující:

Moderní procesory jsou neuvěřitelně rychlé a někdy hloupé řešení s dobrými konstantami může být rychlejší než chytrý algoritmus, který má na papíře lepší asymptoty. Neslavte však předčasně, vždycky nejdřív musíte vědět, čeho je možné dosáhnout.

Dále k tématu:


  1. Pokud mě zajímá n-nejbližších sousedů, pak je budu udržovat v haldě, která mi umožňuje rychle přistoupit k maximu.
  2. Vektor musí mít délku dělitelnou osmi a funkce neřeší, co se stane, když to neplatí. Není problém to dodělat přidáním malé smyčky nebo switche na konci. Dále je třeba mít na paměti, že _mm256_load_pd potřebuje, aby byla paměť zarovnaná na 32 bajtů, jinak segfaultne. Existuje také _mm256_loadu_pd, která načítá nezarovnanou paměť, ale je o něco pomalejší. V tomto případě je třeba místo malloc použít posix_memalign nebo podobnou funkci, která alokuje zarovnanou paměť.
  3. Velikost bloku musí být zvolena tak, aby se data vešla do nějaké úrovně cache. Pokud je tento parametr dimenzován pro L1 cache, je výsledek rychlejší než když je dimenzován pro L2 cache, ale může pojmout méně dat. Nastavení parametru pro sdílenou L3 cache je ještě pomalejší a má navíc ten problém, že L3 cache je většinou sdílená všemi jádry a každé z nich má k dispozici jen odpovídající část kapacity. Hodnota tohoto parametru je závislá na mnoha parametrech a proto ji není dobré nastavovat napevno.
  4. Více informací o OpenMP se dá najít v knize Programming on Parallel Machines
  5. Tady přichází ke slovu curse of dimensionality, která říká, že čím víc rozměrů má daný prostor, tím jsou všechny úhly bližší pravému úhlu a všechny vzdálenosti jsou si podobnější. To má za následek, že algoritmy jako cover tree a podobné fungují hůře a hůře, protože nemají jak prostor rozdělit na výrazně odlišné oblasti.
  6. I když to nemusí být tak úplně pravda na procesorech, které mají micro-op cache. Na ostatních strojích není nadměrná duplikace kódu také příliš vítaná. Když se program nevejde di instrukční cache, může to mít neblahé následky.
  7. Auto-vektorizace na první pohled vypadá jako triviální optimalizace, kterou kompilátor musí dát na první dobrou. Ale ve skutečnosti to ale není vůbec snadné a často je třeba kompilátoru pomoci. Viz: Auto-vectorization with gcc 4.7.
  8. Paralelismus vektorových instrukcí není to samé, co paralelismus out-of-order superskalárních CPU. OOO paralelismus se pokouší najít různé nezávislé instrukce, které může vykonat najednou, vektorové instrukce naproti tomu provedou stejnou operaci na několika datových položkách najednou. OOO je MIMD, vektorové instrukce jsou SIMD.
  9. Pro všechny, kdo nechtějí používat intrinsics, tu je vector extension GCC.

Flattr this!

Posted in Algo, CPU, Paměť | 2 Comments

Kolize hashů pro mírně pokročilé

Hash tabulka je jedna ze základních datových struktur, byla vynalezena roku 1953 a všichni, kdo programování viděli aspoň z rychlíku, ji velice dobře znají. Hash tabulky (někdy také hashmapy nebo slovníky/dicti­onary) umí vyhledat, přidat nebo smazat hodnotu asociovanou s určitým klíčem v konstantním čase. Všechno je možné jen díky jednomu velice chytrému triku – použití hashovací funkce, která klíč, ať už je jakéhokoli typu a velikosti, zobrazí na číslo. To je potom využito k adresování do obyčejného plochého pole, které všichni známe a milujeme.

Tohle všechno funguje skoro až zázračně. Tedy aspoň do okamžiku, kdy to přestane. Pak se dostáváme do velkých problémů.

Všechno záleží na použité hashovací funkci. Ta by měla v ideálním případě být pseudo-náhodná – měla by produkovat hashe, které jsou rovnoměrně rozložené v oboru hodnot a které nezachovávají korelace mezi vstupními hodnotami, kdy např. hashe posloupnosti čísel také připomínají posloupnost. Jinými slovy: mělo by jít alespoň o jednu z rodiny univerzálních hashovacích funkcí.

Pokud hashovací funkce produkuje nadprůměrné množství stejných hashů pro různé hodnoty, nebo je předvídatelná a hodnoty s kolidujícími hashi může někdo, jehož záměry jsou méně než upřímné, snadno spočítat, hash-tabulky najednou začnou vykazovat patologické chování, kdy dříve konstantní operace potřebují lineární čas úměrný velikosti tabulky. Může tak jít o velice efektivní DOS útok.

Abych to osvětlil, musím se ponořit do krvavých detailů.


Hash tabulky mají tři základní varianty: chaining, probing a cuckoo hashing. Chaining a probing se také označují jako open hashing a open addressing, ale tohle označení nemám rád, protože si nikdy nemůžu vzpomenout, které je které.

Řetězící (chaining) varianta se zakládá na poli pointerů (buckets). Každý z nich ukazuje na spojový seznam, jehož každý článek obsahuje klíč, odpovídající hodnotu a pointer na následující článek (nebo null). Když chci vložit pár (key, value), spočítám hash klíče hash = ϝ(key) a z něj odpovídající pozici v poli pointerů index = hash % arraySize. Pokud je daná pozice prázdná, vytvořím článek a nastavím pointer na pozici index, aby na něj ukazoval. Pokud už je index obsazený, vytvořený článek připojím na konec seznamu.

Aby tohle fungovalo v konstantním čase, kolize musejí být nepravděpodobné. Pokud to platí, většina spojových seznamů obsahuje jen málo článků (v průměru n/m, kde m je velikost tabulky a n počet vložených elementů) a pole pointerů má délku úměrnou počtu položek v tabulce. V praxi jde o určitý malý násobek a je dovolené jen určité maximální zaplnění tabulky. Když je tato mez překročena, tabulka se zvětší na dvojnásobek, všechny klíče se rehashují na nové pozice a řetězy kolizí se zkrátí. Zvětšování o násobek garantuje amortizovaný konstantní čas operace vkládání.


Probing varianta typicky používá tři interní pole – jedno obsahuje klíče, druhé odpovídající hodnoty a třetí označuje smazané položky.

Pokud chci vložit klíč, stejně jako v minulém případě spočítám jeho hash a z něho odpovídající index v trojici stejně dlouhých polí. Když je odpovídající místo prázdné, klíč a hodnotu tam jednoduše vložím. Pokud je však místo už zabrané, vložím klíč a hodnotu o jedno místo doprava (linear probing), o 2i míst doprava (quadratic probing) nebo o počet míst vypočítaný jinou hashovací funkcí (double hashing). Pokud je plno i tam, zkouším další místa napravo dokud nenajdu volný flek.

Když chci najít určitý klíč, spočítám hash a z něho pozici, na které by se měl nacházet. Pokud se na tomto místě nachází, skončím s úspěchem. Pokud se tam nachází jiný klíč, zkusím to o jedno místo vpravo. Takto pokračuji, dokud nenajdu požadovaný klíč, smazaná místa ignoruji. Když narazím na prázdnou hodnotu, prohlásím, že klíč nebyl nalezen. Šipky v diagramu ukazují, kde se nachází další klíč se stejným hashem. Jak je vidět tento blok nemusí být spojitý a mezi klíči, které mají shodný hash, můžou být vklíněny klíče s jinými hashi. Záleží jen na pořadí vložení.

Pokud chci klíč smazat, pokusím se nalézt místo, kde se nachází (postup stejný jako v minulém odstavci) a potom, co ho objevím, ve třetím poli označím danou pozici jako smazanou (jde o pole boolů nebo o bitmapu). Kdybych klíče mazal tak, že bych odpovídající pozici v prvních dvou polích jen nastavil na null, rozbil bych hledání, protože bych tak do spojitých bloků vnesl díry a hledání by při nalezení této díry skončilo a nemuselo by se dostat ke správným hodnotám, které se můžou nacházet těsně za touto dírou.3

Konstantní čas pro všechny operace je zaručen opět tím, že hash funkce negeneruje kolize a rovnoměrně rozkládá hashe. V důsledku toho budou sekvence, které je potřeba projít při hledání klíčů, rozumně krátké. Tři zmíněná pole mají délku která je opět úměrná počtu položek v tabulce, i když bývají větší než u chainingu, protože efektivita probingu dramaticky degraduje, když se zaplnění tabulky blíží 100 procentům. Dlouhé běhy obsazených pozic, kam spadne mnoho hashů, která náhodou byly blízko sebe, se ještě více prodlužují a s nimi i počet míst, které je třeba projít. Nestačí jen jako v případě chainingu projít jeden spojový seznam odpovídající jednomu hashi, ale blok odpovídající několika hashům.

(Technicky vzato: Hashovací funkce použitá pro linear probing musí být aspoň 5-independent nebo musí jít o tabulation hashing, jinak není garantován konstantní čas. Více v tomto videu.)


Třetí variantou je takzvaný cuckoo hashing, který je založen na myšlence, že nepoužívám jednu hash funkci, ale dvě (a někdy i více). Výsledek každé funkce odpovídá jedné pozici v tabulce. Při vkládání se pokusím vložit na první pozici a když je obsazená, použiji druhou. Když je na této pozici jiný klíč, ten vezmu a přesunu ho na jeho druhou pozici a tak dále dokud nevystrčím klíč do neobsazeného místa. Pokud jsou hash funkce dobré, kolize řídké a nedojde k vytvoření smyček, počet přestěhovaných elementů při každém vložení je malý a to teoreticky garantuje konstantní časy jednotlivých operací (amortizované, když přihlédnu ke zvětšování tabulky, které probíhá, ne když je tabulka plná, ale když nemůžu vložit nový klíč). Cuckoo hashing se navíc chová rozumně i když je hash tabulka téměř plná a jeho výkon nedegraduje tak drasticky.

Podle toho co jsem četl, Cuckoo hashing se v praxi příliš nepoužívá, protože jeho výkon není z hlediska CPU cache nijak zázračný. Při vkládání musím načíst několik položek, které se nacházejí na prakticky náhodných místech v paměti (to je důsledek dobré hash funkce, která se chová jako PRNG), které s největší pravděpodobností nebudou v cache a z toho důvodu utrpím několik cache-miss. To samé platí při hledání klíče, jen počet načtených míst je v nejhorším případě stejný jako počet hash funkcí.

Chaining hash tabulka na tom také není nijak slavně: Nejprve musím načíst jeden pointer na prakticky náhodné pozici v paměti a pak postupně všechny objekty ze spojového seznamu. Protože články jsou dynamicky alokovány a lezí někde na haldě, každý pointer vede na nepředvídatelné místo v paměti. Pokud je tabulka větší než cache, s největší pravděpodobností toto místo nebude v cache. Navíc je načítání objektů spojového seznamu datově závislé: Musím načíst jeden objekt, abych zjistil pointer na ten další a teprve potom ho můžu načíst. Hardware nemá možnost paralelního přístupu k paměti. Některým z těchto nedostatků se dá předejít, když je spojový seznam částečně odrolovaný (jak je například použito v implementaci hashmap v Go).

Probing verze se vzhledem k hardwaru chová nejlépe ze všech. Sice také potřebuje načíst data ze tří prakticky náhodných míst paměti, ale tyto dotazy jsou nezávislé a hardware je může provést paralelně (paralelní fetch je demonstrován např zde). Změnami algoritmu a vnitřních datových struktur, je možné snížit počet interních polí na dvě nebo dokonce jen na jedno, což má za následek ještě lepší chování z hlediska CPU cache, jak ukazuje například tento článek. Navíc: Kolidující klíče nejsou ve spojovém seznamu, který skáče z místa na místo, ale jsou naskládány v paměti jeden vedle druhého (aspoň v případě lineárního probingu). S velkou pravděpodobností se budou nacházet na jedné cache-line, kterou procesor načte z paměti najednou. A když budu muset projít dlouhý seznam klidujících klíčů, prefetching se postará o to, abych měl data k dispozici včas a bez čekání na pomalou RAM. Počet kolizí může být velký, ale když hardware začne kouzlit, cena jedné kolize bude menší než v případě prolézání kratšího spojového seznamu.

Existují pochopitelně i další způsoby jak implementovat strukturu, která se chová jako asociativní pole. Ty jsou většinou založené na stromech a mají svoje specifické použití, jako například Hash Array Mapped Trie (také tady a tady).


Z předešlých odstavců je jasné, proč je útok na kolizi hashů (v literatuře se používá termín hash flooding) vůbec možný a jak efektivní může být. Všechny varianty hash tabulek se spoléhají na hashovací funkci, která musí být dobrá, rovnoměrná a produkovat málo kolizí. Když tyto předpoklady z nějakých důvodů přestane platit na horizontu se objeví čtyři jezdci apokalypsy a zátěž procesoru vyskočí na 100%. Jak se říká: Tvoje hash tabulka je jenom tak dobrá, jak dobrá je tvoje hash funkce.

Pokud do tabulky vkládám klíče, které mají stejný hash nastane právě takový apokalyptický scénář. A navíc stačí když je hash stejný modulo délka tabulky. V případě řetězení se nové klíče postupně vkládají na konec neustále rostoucího spojového seznamu a v případě probingu se musí projít čím dál tím větší blok obsazených pozic než je nalezeno volné místo, kam se může nováček vložit.

Například PHP pro hashování stringových klíčů používá velice jednoduchou a rychlou funkci DJBX33A, která vypadá takhle:

long DJBX33A(const unsigned char *key) {
  long hash = 5381;
  while(*key) hash = 33*hash + *key++;
  return hash;
}

Stačí chvíli střídavě zírat do kódu a do ASCII tabulky, než člověk přijde na několik krátkých stringů, které mají stejný hash. Jsou to například tyto tři: "~!", "}B", "|c"

Obecně platí, že stačí vybrat takové znaky, pro jejichž pozice v ASCII tabulce platí následující:

 a    * 33 +  b
(a+1) * 33 + (b-33)
(a+2) * 33 + (b-66)
     ...
(a+i) * 33 + (b-i*33)

Z toho, jak je DJBX33A napsána však vyplývá ještě jedna skutečnost: všechny stringy stejné délky poskládané z těchto střípků mají také stejný hash. Takže následující tři stringy mají také stejný hash:

"~!~!~!", "~!}B|c", "|c|cB|"

Najednou si můžu vygenerovat tolik stringů kolik si moje srdce žádá a ve velkém je začít posílat jako GET nebo POST argumenty a provést tak DOS útok nějaké PHP aplikace, jak bylo demonstrováno na Chaos Communication Congress v roce 2011 a 2012.

PHP tuto zranitelnost napravilo po svém a omezilo počet argumentů předaných jako POST nebo GET na něco kolem tisíce. Na(ne)štěstí v současnosti každá druhá webová služba poskytuje JSON API, které může útočník podstrčit speciálně připravený JSON dokument, který obsahuje jen a pouze zákeřné klíče. Když pak webová aplikace zavolá json_decode1, naparsuje JSON do PHP pole a sama se tak DOSne.

Pro představu jak velkou škodu může útok napáchat. Vytvoření pole se 100000 neškodnými klíči zabere (v PHP 5.6) 40 milisekund, se stejným počtem zákeřných klíčů 84000 milisekund.

V Javě/Scale jsem po rychlém experimentu v REPLu naměřil hodnoty: 27 ms a 20000 ms. Na konkrétních číslech však moc nezáleží, jde o to, že jedno roste lineárně, druhé kvadraticky.

V PHP verze 7 mají hash tabulky drasticky pozměněnou implementaci. Přešli z řetězu dynamicky alokovaných spojových seznamů na dvě před-alokovaná fixní pole, která má výkon velice podobný probingu. Stále však používají hashovací funkci DJBX33A. Útok je tedy stále možný, jen díky zvýšenému výkonu bude o něco málo méně účinný. Podle prvních nástřelů to vypadá, že je zhruba 10× méně nebezpečný.

A to nemluvím o tom, že celočíselné klíče se berou jako hash, což vede k možnosti útoků, které jsou triviální a zcela devastující.1

Jedinou skutečnou obranou proti tomuto útoku je použít randomizovanou hashovací funkci. Protože i když používám funkci, která rovnoměrně rozkládá hashe a není možné jednoduše před-počítat kolidující klíče a pak z nich poskládat mnohem víc delších klíčů, dostatečně odhodlaný útočník může brute-force metodou před-počítat dostatečné množství kolizí a pak je opakovaně posílat na cílový server. Přesně v duchu rčení pain is temporary, glory is ethernal.


Jak jsou na tom jiné jazyky/runtime?

HHVM standardně používá Murmur3 hash, ale když běží ne stroji, který podporuje SSE4.2, tak začne hashovat pomocí CRC32, která je implementována v hardware a tudíž velice rychlá (bohužel je rychlá i pro útočníky).

Java k hashování stringů používá funkci, která je velice podobná té z PHP:

public int hashCode() {
  int h = 0;
  for (char v : value) {
    h = 31 * h + v;
  }
  return h;
}

A stejně jako v případě PHP je zcela triviální vygenerovat kolidující stringy:

"1}" // 49 * 31 + 125 = 1644
"2^" // 50 * 31 + 94  = 1644
"3?" // 51 * 31 + 63  = 1644
"4 " // 52 * 31 + 20  = 1644

Integer se (zase podobně jako v PHP) hashuje na sebe sama. Na rozdíl od PHP se však tato slabina nedá tak jednoduše zneužít, protože int má malý rozsah. V javě má 4 bajty a něco přes 4 miliardy možných hodnot, takže v nejhorším případě do hash tabulky můžu vložit 65536 klíčů, jejichž hash je dělitelný 65536 a hash modulo velikost tabulky bude stejný.

Long se hashuje takto:

public static int hashCode(long value) {
  return (int)(value ^ (value >>> 32));
}

To vypadá rozumně, ale

(long) x << 32 | (long) x

je nula pro libovolné x. Když jsou horní a dolní čtyři bajty stejné, jejich xor je 0.

Zkuste si to sami a uvidíte jak se roztočí větrák na procesoru:

val map = new collection.mutable.HashMap[Long, Int]()

for (i <- 0 until 100000) {
  val h = (i.toLong << 32 | i.toLong)
  map.put(h, 1)
}

Problém může také nastat, když obě poloviny longu používají jen pár dolních bitů. Když používají jen 10 bitů (tedy čísla v rozmezí 0 – 1024), můžou reprezentovat milion kombinací, ale jejich hash je jen v onom desetibitovém rozmezí 0 – 1024. V takovém případě se bude hashovat jen do malé částí hash-tabulky a všechno bude (opět) velice pomalé. Tenhle případ jsem poznal na vlastní kůži, když jsem se snažil udělat něco chytrého.

Java v některé nedávné verzi přinesla jednu vychytávku: Pokud jeden slot obsahuje víc jak 8 kolizí, je spojový seznam přestavěn na binární red-black strom. Ten v patologických případech zredukuje lineární časovou složitost na logaritmickou za cenu nepříjemně komplikovaného kó­du.

C# podle všeho používá celkem jednoduchou multiplikativní funkci.

val hash1, hash2 = 5381
while (i < str.length) {
  hash1 = (hash1 * 33) + str(i)
  hash2 = (hash2 * 33) + str(i+1)
  i += 2
}

hash1 + (hash2 * 1566083941)

Když se nechci uchylovat k bruteforce řešení, je spočítání kolidujících stringů poměrně jednoduché. Například tyto dva kolidují: "*_O_", "+_n_".

Pokud se však zahledíte do odkazovaného kódu, uvidíte tam za několika IFDEFy nějaké ty XORy a možná i nějakou tu randomizaci. Těžko říct o co přesně tam jde.


Ať už je ale hash funkce libovolná, útočník může vždycky bruteforce vypočítat kolidující řetězce. Když jich zkusí moře, dostane pár kapek, a to může být víc než dost.

Jediným systematickým řešením je použít randomizovanou hash funkci. Taková funkce má několik skrytých argumentů (jsou např. vygenerovány na začátku běhu programu), které změní výstup hashovací funkce tak, že když s jedním argumentem mají dvě hodnoty stejný hash, s jiným mají vzájemně odlišné hashe.

Jako špatný příklad příklad uvedu PHP funkci DJBX33A. V té se na začátku vyskytuje magická konstanta 5381. Když bych randomizoval tuto konstantu, hashe budou jiné, ale jejich korelace zůstanou.

Pokud platí:

ϝ(a, 5381) = ϝ(b, 5381)

Tak bude také platit

ϝ(a, s) = ϝ(b, s)

Hashe se změní, kolize zůstanou.

Je nutné použít randomizovanou funkci, která tak byla navržena, jako je třeba SipHash (podrobně popsaná v tomto paperu). Ta je jednak rychlá (pomalejší než jednoduchá DJBX33A, ale mnohem rychlejší než třeba MD5) a zároveň poskytuje všechno dobré kolem randomizace. Útočník je schopný provést hash flood útok jen pokud zná náhodných 128 bitů, které SipHash používá jako seed, ale ty se můžou měnit každou minutu nebo při každém HTTP požadavku. Útočník jednak nemá jak zjistit tento seed a i kdyby ho zjistil 2, má jen limitované možnosti, jak ho uplatnit.

V případě PHP je změna hashovací funkce triviální, protože na rozdíl od světa Javy není součástí veřejného API a její změna nemůže rozbít kód, který na ní závisí.


Jako další ukázka proč jsou hashovací funkce důležité může posloužit příklad Bloom filtrů. Bloom filtr je pravděpodobnostní datová struktura pro test členství v množině, která interně používá několik hashovacích funkcí. Nedávno jsem psal vlastní maličkou implementaci bloom filtru do knihovny sketches a použil jsem stejné hashovací schéma jako v knihovně Breeze založené na MurMur3 hashi. Při testování se ukázalo, že výsledný filtr generuje 1000000000× více falešných pozitiv, než by statisticky měl. To byla velice špatná zpráva, protože to dělalo tuto strukturu mnohem méně přesnou a mnohem méně užitečnou. Později se ukázalo, že za tuto vadu může právě použití MurMur3 hashe, který očividně neprodukuje uniformní rozložení a nepatří tedy do rodiny univerzálních hashovacích funkcí. Když jsem použil dobrou univerzální hashovací funkci, problém zmizel.


Na závěr přidám dvě perličky:

Cache procesoru je také (překvapivě) organizována jako hash-tabulka – obsahuje několik slotů (řekněme 64) a každý slot může obsahovat 8 cache line (každá má 64 bajtů). Hashovací funkce vypadá takto: vezme fyzickou adresu, zahodí 6 nejnižších bitů a vezme 6 dalších. Výsledkem je index slotu v němž bude hledat nebo do něhož bude ukládat data.

Když program přistupuje na adresy s velice specifickým rozestupem, cache alising v podstatě udělá hash flood CPU cache, jak je ukázáno ve What Every Programmer Should Know About Memory nebo v Binary Search Is a Pathological Case for Caches. Následující graf pochází právě z WEPSKAM.

Za pozornost stojí také Robin Hood Hashing. Jde o variantu probing hashování, kde v nejhorším případě mají operace logaritmickou složitost. RH hashování toho docílí tím, že průměruje počet slotů, které je potřeba prohledat ve stylu Robina Hooda – bohatým bere (krátká probing cesta) a chudým dává (dlouhá cesta). Elementy v jenom běhu jsou tedy seřazeny podle délky probe cesty (nebo dokonce v některých variantách podle hashe samotnho), při vkládáná se tedy provádí inserion sort a hledání je možné předčasně ukončit i v případě, když element není v tabulce přítomen a dále to vede k velice jednoduchému a rychlému zvětšování tabulky a hromadnému vkládání, kterému předchází radix sort.

Tohle všechno znamená, že jsou všechny operace rychlé a tabulka se chová rozumně, i když je z velké části (klidně i 90%) zaplněná. Jak o tom teď píšu, napadá mě, že bych měl napsat vlastní implementaci ve Scale a přidat ji do sketches.


A to je všechno. Pokud jste se dočetli až sem, nejspíš toho víte o hash funkcích a jejich kolizích víc, než jste kdy chtěli vědět.


Dále k tématu:


  1. Tady patří díky Janu-Sebastianovi Fabíkovi, který krátce potom, co jsem dopřednášel, přišel na to, že json_decode je částečně odolný nejtriviálnějšímu útoku, kdy se serveru pošle pole s speciálně připravenými numerickými klíči. Když json_decode bez extra parametrů narazí na JSON objekt, naparsuje ho jako PHP objekt. PHP objekty konvertují všechny atributy na string a tedy útok numerickými klíči nemá žádný dopad. Když se však json_decode předá druhý argument $assoc, tak JSON objekty naparsuje na pole a jsme zpátky odkud jsme vyšli. V obou případech je však zranitelný, když mu podstrčíme json s kolidujícími stringovými klíči.
  2. např. Python používal podobnou funkci, ale později se ukázalo, že je možné zjistit tajný seed.
  3. Přímo je možné smazat elementy na konci řady, tedy ty, které mají napravo neobsazenou pozici. Stejně tak je možné vložit element do smazané pozice. Tyto dvě zkratky nerozbijí hledání, protože nezpřetrhají řady obsazených pozic.

Flattr this!

Posted in DS, Hashování, Java, PHP | 1 Comment