0xDEADBEEF

Floating point triky

IEEE 754 nařizuje, že 32 bitové číslo v plovoucí řádové čárce, také známé jako float, v paměti vypadá takhle:

+eeeeeeee.......................
^ znaménko (1 bit)
 ^ exponent (8 bitů)
         ^ zbytek (23 bitů)

Hodnota se z binární reprezentace vypočítá podle vzorce -1^znaménko * 2^{(exponent-127)} + (1.zbytek).

Zbytek neboli mantisa představuje binární cifry za čárkou. To je důvod, proč float nemůže přesně reprezentovat desetinná čísla jako 0.1. Žádná kombinace polovin, čtvrtin a 2ⁿ-tin nedá dohromady přesně 0.1.

Za pozornost stojí drobnost, že sčítání dvou floatů se stejným exponentem, sečte 23 bitů zbytku. Sčítání floatů s různými exponenty se udělá posunem jednoho na společný exponent, pak se sečtou bity za čárkou a výsledek se normalizuje na společný exponent. Podobně tak to platí pro násobení, jen mantisa se násobí a exponent sčítá. Ve výsledku to znamená, že v každém floatu se nachází malý 23-bit integer.

union IF { uint i; float f; }

float intToFloat(uint i) { IF x = { i: i }; return x.f; }
uint floatToInt(float f) { IF x = { f: f }; return x.i; }

uint a = args[1].to!uint;
uint b = args[2].to!uint;

float f = intToFloat(a);
float g = intToFloat(b);

auto intSum   = a + b;
auto floatSum = floatToInt(f + g);

assert(isum == fsum);

Jaký to má užitek? Prakticky žádný. Proč to píšu? Je to jenom zajímavé.

Může se to hodit v krajních případech, kdy mám procesor se SIMD rozšířením, které umí pracovat s floaty, ale ne s inty a zrovna bych potřeboval vynásobit/sečíst mnoho intů. Taková situace existovala v době před příchodem AVX2. Teď už je to jen teoretické cvičení.

Nebo možná to může najít uplatnění na určitých grafických kartách, které mají víc ALU pro floaty než pro inty. Gen7.5 grafika v mém Haswellu má 2 SIMD4 ALU, obě zvládají floaty, jen jedno int operace.

Na druhou stranu, že to není úplná zbytečnost, dosvědčuje jedno z mnoha rozšíření AVX-512. To přidává instrukce pro sčítání a násobení vektorů 8 celých čísel o šířce 52 bitů. Proč zrovna 52? Double má právě tolik bitů v mantise. Úvaha byla nejspíš taková, že hardware tak jako tak nese sčítačky a násobičky pro 52 bitů, ale využívá je jen pro floaty. Takhle přihodí pár instrukcí a pár drátů a existující hardware najde využití i v dalších programech.

V podobném duchu RDNA 2 GPU architektura specifikuje instrukci V_MUL_I32_I24, která násobí dva 24-bitové integery za použití efektivních násobiček pro floaty.

Ok, když už nic, přidám další trik, který může někdě najít okrajové uplatnění.

Pokud mám velké množství floatů v rozmezí mezi 0 a 1, ale málo místa a na přesnosti příliš nezáleží, můžu je velice jednoduše smrsknout na 2 bajty:

ushort smush(float f) {
  union IF { float f; uint i; }
  IF x = { f: f };
  return cast(ushort) (x.i >> 14);
}

nebo následujícím způsobem, to ale (myslím) není validní D. Type punning se může dělat jen přes union.

ushort smush(float f) {
  return cast(ushort) (*(cast(uint*) &f) >> 14);
}

Co to dělá? Funkce na to jde z druhé strany, od exponentu. Exponent určuje dynamický rozsah, mantisa přesnost. Když to posunu doprava, ukousnu dolní bity jemné přesnosti. Něco z mantisy zůstane a určitá přesnost bude zachována. Při posunu jen o 14 bitů se navíc ztratí horní dva bity – znaménko a horní bit exponentu. Ty pro čísla mezi 0-1 budou oba vždy nulové a o žádnou informaci nepřijdu.

+eeeeeeee.......................
  ^^^^^^^^^^^^^^^^
  tyto bity zbydou

Dovedu si představit, že když bych chtěl na celkem primitivním hardwaru jako RasPi bez podpory polovičních floatů¹ rozběhnout nějakou podstatnou aplikaci, tohle by se mohlo v určité formě hodit. Možná. Nejde ale o ideální mapování.

1:    512 hodnot mezi 1 a 1/2
2:    512 hodnot mezi 1/2 a 1/4
3:    512 hodnot mezi 1/4 a 1/8
4:    512 hodnot mezi 1/8 a 1/16

...

128:  512 hodnot mezi 2.3/100000000000000000000000000000000000000000 a 0

Větší přesnost se dá docílit použitím menší části exponentu a větší části mantisy. To ale není tak triviální jako trik výše a v podstatě jde o 16 bit half float. (Pro konverzi floatů na half floaty a zpět má x86 instrukce VCVTPS2PH a VCVTPH2PS).

Tipoval bych, že to bude rychlejší než kvantizace floatů na inty² . Na x86 jsou cvtsi2ss, rcpss pomalejší než logický shift, ať se dívám na jakoukoli verzi procesoru.

Dodatek: Ještě jeden trik, který hraničí se zločinem. Dolní bity určují nejjemnější rozlišení přesnosti, když na něm extra nezáleží, můžeme jich pár použít k něčemu jinému. Třeba do nich zakódovat id. Aritmetika na takto poskvrněných floatech pořád funguje (jen s malou nepřesností). Můžeme například seřadit pole a dostat správný výsledek a přehazovat jen 4B floatů, namísto 8B, které by zabral celý pár.

1. Na x86 se o konverzi mezi floaty a polovičními floaty stará rozšíření F16C, které je dostupné ve všem hardware uvedeném od roku 2013.
2. Q (number format)