99.00000000000000009 problémů s floating point čísly

24. 1. 2017

Floating point čísla a jejich IEEE 754 varianta jsou jednou z těch věcí, které mě nikdy nepřestanou fascinovat. Jde o užitečnou věc, která člověka skoro přesvědčí, že všechno bude ok, že svět má svůj řád, že se stačí odevzdat floating vlnám a už se není třeba o nic starat. Hned potom ale narazí na zaokrouhlovací chyby, nekonečna, NaNy, negativní nulu, denormální čísla a hned z toho narkotického blaha vystřízliví.

V tomto duchu začal 13. sraz přátel PHP - historkou o klasickém případě, kdy floaty nepoužívat - k reprezentaci peněz.

Floaty pochopitelně nejsou přesné a zaokrouhlovací chyby vznikají i v triviálních případech. Například 0.1 + 0.2 je 0.30000000000000004. Pokud to šlo peníze, najednou jsem z ničeho vytvořil čtyři femtohalíře.

Tyto chyby jsou malé a nenápadné a člověk si je uvědomí až ve chvíli, kdy mu kancelář začnou obléhat auditoři s beranidly a pokřikují něco o tom, že nesedí účetnictví. Kdyby se peníze ukládaly do intů, tak v nejhorším případě částka přeteče do velkého množství záporných peněz a to je chyba, které si nově vzniklý záporný milionář snadno všimne.

Samozřejmě, na peníze a jiné kvantity, které je třeba vést naprosto přesně, se musí použít něco jako BigInteger nebo Big Decimal.

A to je jenom začátek.

Některé jazyky jako třeba PHP nebo JavaScript přidávají svoje speciality, protože v nich celá čísla automaticky přetékají do floatů.

Například tento PHP kód

for ($i = 0; $i != ($i+1); $i++) { echo "lol"; }

vypadá jak nekonečná smyčka, protože se přece nikdy nemůže stát, že $i je stejné jako $i+1. PHP na to má ale jiný názor. Int eventuálně přeteče do double floatu s nižší přesností, které mají krok mezi dvěma sousedními čísly větší než 1, součet nic nezmění a smyčka skončí. V PHP je int jen double, který čeká na svůj okamžik slávy.

Ve slušných jazycích, by int přetekl do záporných hodnot, ale pořád by platilo, že i a i+1 jsou rozdílné.

Pro zajímavost: celá čísla v intervalu [-2²⁴, 2²⁴] můžou být v single precision floatu uloženy přesně bez zaokrouhlovacích chyb.

Asociativita je užitečná vlastnost, která je stěžejní pro paralelizaci. Tomu si, jak doufám, žádný slušný člověk nedovolí odporovat. Floaty pochopitelně asociativní nejsou.

val ds = Array.fill(10000)(util.Random.nextDouble)
assert(ds.sum - util.Random.shuffle(ds).sum == 0.0)

val ds = Array.fill(10000) {
  val a, b, c = util.Random.nextDouble
  ((a+b)+c) - (a+(b+c))
}
assert(ds.forall(d => d == 0.0))

Obě aserce selžou, protože součty prováděné v různých pořadích, uvedou různé zaokrouhlovací chyby a to vede k různým finálním výsledkům.

Tohle je poměrně důležité, protože to kompilátorům svazuje ruce v optimalizacích, které můžou bezpečně udělat. Změna pořadí operací může vést ke zrychlení, ale taky k odlišnému výsledku.¹

IEEE 754 specifikuje bitovou reprezentaci floatů a toho se dá využít k různým trikům. Jedním z nich je možnost transformovat float na int s tím, že vztah mezi floaty je stejný jako mezi nově vytvořenými inty. Tedy, když platí

a < b

tak platí i

floatToInt(a) < floatToInt(b)

a zároveň existuje transformace zpátky

intToFloat(floatToInt(f)) = f

Jde o pár operací, které ve Scale vypadají takto:

def flip(f: Int) = f ^ (-(f >>> 31) & 0x7FFFFFFF)

def floatToInt(f: Float) = flip(Float.floatToRawIntBits(f))
def intToFloat(i: Int)   = Float.intBitsToFloat(flip(i))

Kód je založený na tomto článku, jen jsem ho upravil pro znaménkové inty Javy. V podstatě jde o to, že pokud je float záporný, převrátím všechny bity kromě znaménka. Tím se bitová reprezentace bude chovat stejně jako int.

Když vynesu průběh intů a floatů na intervalu všech bitových hodnot od samých nul nalevo až po samé jedničky napravo, inty mají takovýhle průběh:

Floaty zase takovýto (NaNy by měly být někde za nekonečny, ale ty si teď dovolím ignorovat):

Kód jen překlopí druhou polovinu intervalu.

Tohohle triku se dá využít k jednoduchému řazení floatů radix sortem. Na začátku floaty převedu na inty zachovávající pořadí, rozjedu radix sort pro inty a na konci výsledek přepíšu zpátky na foaty. Existují i přímé způsoby, jak radixsortovat floaty, ale ty jsou komplikovanější.

Dále k tématu

Pozn.

Některé architektury můžou provádět výpočty s vyšší přesností, třeba 80 bitů namísto 64 a výsledek nakonec zaokrouhlit na 64, viz strictpf.