Tajemnice liczb zmiennoprzecinkowych

English | Polish

Ten artykuł został oryginalnie opublikowany w numerze 5/2024 (115) (listopad/grudzień 2024) magazynu Programista.

W tym artykule omówię liczby zmiennoprzecinkowe zgodne ze standardem IEEE 754, które są dostępne w większości języków programowania, a także opiszę ich budowę, możliwości i ograniczenia. Odniosę się również do powszechnie panujących przekonań, że liczby te są niedokładne lub że są niedeterministyczne. Pokażę ponadto wiele nieoczywistych pułapek, które czyhają na używających ich programistów. Ta wiedza może się przydać każdemu, obojętnie, w jakim języku i na jakie platformy programujesz!

Liczby zmiennoprzecinkowe to świetny wynalazek. Dostępne w językach programowania typy danych tego rodzaju pozwalają na niewielkiej liczbie bitów zapisywać liczby dodatnie, jak i ujemne, całkowite, jak i ułamkowe, bardzo małe (bliskie zeru), jak i bardzo duże. Ponieważ komputery na początku służyły głównie naukowcom do wykonywania obliczeń, historia różnych sposobów kodowania liczb jest tak stara, jak historia samej informatyki. Standard IEEE 754, na którym opierają się używane współcześnie liczby zmiennoprzecinkowe, pojawił się w roku 1985. To jego używają wszystkie typy danych, o których będziemy mówić w tym artykule.

Podczas nauki programowania i w miarę nabierania doświadczenia w używaniu wybranego języka programowania zrozumienie dla liczb zmiennoprzecinkowych można podzielić na trzy etapy: Najpierw programista używa takich liczb bez zastanowienia, przyjmując, że reprezentują one dowolne liczby rzeczywiste. Szybko jednak może się natknąć na trudności z tym związane i zaobserwować różne wynikające z tego błędy. Wówczas przychodzi moment na refleksję nad ograniczeniami tego typu liczb. Stosowanie prostych „zasad ograniczonego zaufania” do liczb zmiennoprzecinkowych pozwala uniknąć wielu błędów. Są to zasady takie jak:

  • „Liczby zmiennoprzecinkowe są niedokładne.”
  • „Liczby zmiennoprzecinkowe są niedeterministyczne.”
  • „Liczb zmiennoprzecinkowych nie wolno porównywać operatorem ==.”

Warto jednak zagłębić się bardziej w szczegóły budowy i działania tych typów danych. Wówczas, dzięki ich lepszemu zrozumieniu, można ich używać świadomie, unikać błędów bądź minimalizować szanse ich wystąpienia, a równocześnie wykorzystywać pełnię ich możliwości. Ten artykuł ma za zadanie wprowadzić czytelnika w tajniki liczb zmiennoprzecinkowych, pokazać różne pułapki i nieoczywiste zjawiska z nimi związane.

Podstawy

Zasada działania liczb zmiennoprzecinkowych przypomina notację naukową używaną na lekcjach fizyki w szkole. O ile liczby, którymi posługujemy się na co dzień, np. licząc pieniądze, można zapisywać zwyczajnie, pisząc np. 12300, o tyle w fizyce często mamy do czynienia z liczbami bardzo dużymi lub bardzo małymi. Na przykład odległość Ziemi od Słońca wynosi około 150 milionów kilometrów, co zamiast 150 000 000 możemy wygodniej zapisać jako 1.5 * 108. Tak zapisana liczba jest niejako „znormalizowana” do postaci, w której przed przecinkiem została tylko pierwsza cyfra, natomiast wykładnik potęgi 10x służy do tego, aby przesunąć przecinek o wybraną liczbę miejsc dziesiętnych w prawo (zwiększając wartość liczbową), kiedy jest wykładnik dodatni, lub w lewo (zmniejszając wartość), kiedy jest ujemny.

Liczby zmiennoprzecinkowe działają na podobnej zasadzie, jak taka notacja naukowa w fizyce, ale na liczbach binarnych, a nie na dziesiętnych. W pamięci komputera mają budowę nieco bardziej skomplikowaną od liczb całkowitych. Liczby całkowite (ang. integer) kolejne kombinacje bitów przeznaczają do reprezentowania kolejnych wartości: 0, 1, 2, … Liczby zmiennoprzecinkowe natomiast dzielą cały ciąg bitów na trzy części. Przykład pokazano na Rysunku 1. Format danych, którego tutaj używamy dla przykładu, to liczba zmiennoprzecinkowa 32-bitowa nazywana liczbą „pojedynczej” precyzji (ang. single precision), która przeznacza 1 bit na znak, 8 bitów na wykładnik i 23 bity na mantysę.

Rysunek 1. Dekodowanie liczby zmiennoprzecinkowej

Zaczynając od najmłodszych bitów, mantysa (ang. mantissa, nazywana też significand, pokazana tu na zielono) koduje cyfry po przecinku. Liczba jest znormalizowana tak, aby przed przecinkiem znajdowała się pojedyncza jedynka. Toteż możemy tę jedynkę pozostawić w domyśle i nie zapisywać jej w pamięci, zyskując tym samym jeden bit precyzji. (Wyjątkiem są liczby zdenormalizowane, w których ta domyślna jedynka nie występuje. Wspomnimy o nich jeszcze, omawiając wartości specjalne.)

Tak zakodowana liczba leży w zakresie od 1 do 2. Aby ją zwiększyć lub zmniejszyć, przesuwając przecinek, kolejna część koduje wykładnik (ang. exponent, pokazany tu na czerwono) dla potęgi dwójki, przez którą liczba zostaje pomnożona. Mając do dyspozycji 8 bitów, interpretowanych jako liczba całkowita, wykładnik byłby w zakresie 0…255. Jednakże stosuje się do niego przesunięcie (ang. bias), które w przypadku tego typu 32-bitowego wynosi 127. Przesunięcie to jest odejmowane od wartości wykładnika, aby móc zakodować zarówno wartości dodatnie, jak i ujemne. Wartość wykładnika zawierająca same zera lub same jedynki jest zarezerwowana na wartości specjalne, więc pozostaje nam zakres wykładnika [1…254] - 127 = [-126…127].

Wreszcie, najstarszy bit to bit znaku (ang. sign bit, oznaczony na niebiesko), który pozwala na zapisywanie liczb ujemnych. Pokazana na Rysunku 1 liczba -1 podniesione do potęgi s to nic innego, jak mądry sposób na zapisanie prostej zasady, że bit ten odpowiada na pytanie: „Czy jest minus?” Kiedy wartością bitu jest 0, minusa nie ma, więc liczba jest dodatnia. Kiedy bit jest ustawiony na 1, wówczas dopisujemy minus i liczba jest ujemna.

Zauważmy, że mając osobny bit reprezentujący znak, liczby zmiennoprzecinkowe są symetryczne względem zera. Aby zanegować liczbę, wystarczy zanegować jej najstarszy bit. Aby obliczyć wartość bezwzględną, wystarczy wyzerować ten bit. Działa to inaczej, niż w przypadku liczb całkowitych ze znakiem, które stosują notację „uzupełnienia do dwóch” (U2). W liczbach całkowitych ze znakiem najstarszy bit również jest ustawiony wtedy, kiedy liczba jest ujemna, jednak ich zanegowanie nie jest już tak trywialne.

Z takiej budowy liczb zmiennoprzecinkowych wynika jeszcze jedna ich ciekawa właściwość: Zmieniając tylko bit znaku, możemy uzyskać dwie różne reprezentacje zera, nazywane +0 i -0. Takie wartości faktycznie istnieją i mogą nieść pewną informację, np. czy liczba, która po obliczeniach stała się zbyt mała i musiała zostać zapisana jako 0, była dodatnia, czy ujemna. Ponieważ jednak te dwie wartości przy porównaniu są sobie równe, to praktycznie nigdy nie musimy się tym przejmować.

Wsparcie sprzętowe

Liczby zmiennoprzecinkowe, przy całej swojej skomplikowanej budowie, nie miałyby większego sensu w programowaniu, gdyby obliczenia na nich nie były wspierane w sposób sprzętowy, aby mogły być wykonywane szybko. Takiego wsparcia w dawnych czasach dostarczał osobny, instalowany opcjonalnie koprocesor matematyczny: Floating Point Unit (FPU), jak układ 8087 w dawnych pecetach. Obecnie jednak procesory (CPU), jak i układy graficzne (GPU) w naszych komputerach mają wbudowaną obsługę takich liczb.

Jeśli chodzi o wspierane sprzętowo typy danych, to najpopularniejsze są liczby zmiennoprzecinkowe 32-bitowe pojedynczej precyzji (ang. single) oraz 64-bitowe podwójnej precyzji (ang. double). Ich możliwości zaprezentowano w Tabeli 1, której rozszerzoną wersję można też znaleźć online w [1]. Ich dostępność i nazwa bywa różna w różnych językach programowania. Na przykład języki C i C++ definiują typy float i double. Ich szerokość nie jest co prawda ściśle określona, ale na większości platform odpowiadają one wspomnianym formatom 32b i 64b. Niektóre języki skryptowe natomiast (w tym JavaScript, Lua) oferują tylko jeden typ liczbowy, który implementowany jest za pomocą typu double i który może być używany nawet do zapisywania liczb całkowitych. Nazwy typów zmiennoprzecinkowych w różnych językach programowania podsumowuje Tabela 2. Niektóre języki oferują także typy o jeszcze większej precyzji, jak long double w C/C++ czy decimal w C#, nie pokazane w tabeli.

Tabela 1. Możliwości liczb 32- i 64-bitowych

Nazwa Pojedynczej precyzji, single Podwójnej precyzji, double
Bitów: suma = znak + wykładnik + mantysa 32 = 1 + 8 + 23 64 = 1 + 11 + 52
Wykładnik: przesunięcie, zakres 127, -126…127 1023, -1022…1023
Precyzja - dziesiętne cyfry znaczące 7.22 (6…9) 15.95 (15…17)
Najmniejsza liczba zdenormalizowana 2−149 ≈ 1.40 × 10−45 2−1074 ≈ 4.94 × 10−324
Najmniejsza liczba znormalizowana 2−126 ≈ 1.18 × 10−38 2−1022 ≈ 2.23 × 10−308
Następna wartość po 1 1 + 2−23 ≈ 1.00000012 1 + 2−52 ≈ 1.00000000000000022
Liczby całkowite reprezentowane dokładnie 0…224 = 16 777 216 0…253 = 9 007 199 254 740 992
Maksimum (2 − 2−23) × 2127 ≈ 3.40 × 1038 (2 − 2−52) × 21023 ≈ 1.80 × 10308

Tabela 2. Typy zmiennoprzecinkowe w językach programowania

Język32-bit64-bit
C, C++ *floatdouble
C#floatdouble
Javafloatdouble
Delphi / Object PascalSingleDouble
Python * float
JavaScript Number
Lua ** number

* Zazwyczaj, zależnie od implementacji.
** Jedyny obsługiwany typ liczbowy.

Czy te liczby są niedokładne?

Popularny pogląd mówi, że liczby zmiennoprzecinkowe są niedokładne i w związku z tym nie należy ich nigdy porównywać operatorem ==. Stosowanie takiej zasady pomaga uniknąć wielu błędów, ale jest to pewne uproszczenie. Obliczenie 2.0 + 2.0 zawsze zwróci 4.0, a nigdy 3.99999. Przyjrzyjmy się zatem bliżej tematowi precyzji liczb zmiennoprzecinkowych, aby lepiej go zrozumieć.

To oczywiste, że liczby zmiennoprzecinkowe są tylko pewnym przybliżeniem liczb rzeczywistych z matematyki, a zakodowane na określonej liczbie bitów z natury mają pewną skończoną precyzję. Kiedy mówimy o precyzji, mamy na myśli najmniejsze różnice w wartości, które możemy zapisać i rozróżnić w określonym typie danych. W przypadku liczb całkowitych sprawa jest prosta: Precyzja zawsze wynosi 1 w całym zakresie, więc po wartości 5 następna możliwa do zapisania wartość to 6, dalej 7 itd.

Można też sobie wyobrazić format obsługujący ułamki, ale stałoprzecinkowy, na przykład używający starszych 8 bitów do kodowania cyfr binarnych przed przecinkiem, a młodszych 8 bitów do kodowania cyfr po przecinku. Takie formaty stałoprzecinkowe nie są obsługiwane sprzętowo przez typowe procesory w komputerach, więc obliczenia na nich trzeba by emulować programowo i działały by one wolno. Dlatego też nie są zbyt często wykorzystywane. Dla lepszego zrozumienia warto jednak wyobrazić sobie taki format i zauważyć, że miałby on precyzję co do 2-8 = 1/256 ≈ 0.0039, również jednakową w całym swoim obsługiwanym zakresie.

W przeciwieństwie do liczb całkowitych czy stałoprzecinkowych liczby zmiennoprzecinkowe mają precyzję zależną od tego, jak duża jest wartość liczbowa. Symbolicznie ilustruje to Rysunek 2. Kolejne wartości możliwe do zakodowania, zaznaczone na osi liczbowej, są rozmieszczone gęściej w pobliżu zera i tym rzadziej, im większa (co do wartości bezwzględnej) jest dana liczba. Z każdą wyższą wartością wykładnika precyzja spada 2 razy.

Rysunek 2. Liczby zmiennoprzecinkowe na osi liczbowej

Możemy więc tutaj raczej mówić o precyzji w sensie liczby cyfr znaczących. Ile cyfr binarnych precyzji mamy do dyspozycji, to wyznacza oczywiście liczba bitów mantysy. W systemie dziesiętnym można ją oszacować na minimum 6 cyfr znaczących dla typu float i 15 cyfr dla typu double (wg Tabeli 1). Na przykład jeśli wartość liczby 32-bitowej wynosi 10.5, to następna możliwa do zapisania po niej wynosi 10.500001, a więc mamy tu precyzję wynoszącą około 1/1000000. Jeżeli jednak wartość jest większa, np. 32000, to precyzja wynosi już 0.004. Przykładowy kod w języku C ilustrujący ten problem pokazuje Listing 1.

Listing 1. Demonstracja skończonej precyzji

float x = 1111.111111111111f;
printf("%.16g\n", x);
// Wynik: 1111.111083984375

Poprzednie dwa akapity warto przeczytać jeszcze raz i dobrze zrozumieć, bo wynika z nich wiele daleko idących konsekwencji w postaci ograniczeń liczb zmiennoprzecinkowych i pułapek, jakie czyhają na programistę. Zajmijmy się nimi teraz.

Reprezentowanie liczb całkowitych

Po pierwsze, warto zadać sobie pytanie, jak liczby zmiennoprzecinkowe potrafią reprezentować liczby całkowite. Okazuje się, że robią to całkiem dobrze, ale tylko w ograniczonym zakresie. Jeżeli pominiemy fakt, że obsługują one także ułamki, a będziemy posługiwać się wyłącznie liczbami całkowitymi, wykonywać wyłącznie dodawanie, odejmowanie, mnożenie, a wynik dzielenia zawsze zaokrąglać w dół lub w górę do liczby całkowitej, to nasze wyniki będą dokładne. Możemy na nich polegać, porównywać je operatorem == itd. Nic dziwnego więc, że w niektórych językach skryptowych typ zmiennoprzecinkowy 64-bitowy może być używany jako jedyny dostępny typ liczbowy, także do zapisywania liczb całkowitych.

Jednak ze względu na przeznaczenie określonej liczby bitów na znak, wykładnik i mantysę, możliwość dokładnego zapisywania liczb całkowitych przez liczby zmiennoprzecinkowe jest ograniczona w porównaniu z liczbami typu int. Na przykład 32-bitowy int ze znakiem obsługuje wartości aż do maksymalnej 2 147 483 648 (ponad 2 miliardy), podczas gdy 32-bitowy float jako największą wartość całkowitą reprezentowaną dokładnie ma 16 777 216 (ponad 16 milionów). Powyżej tej wartości utrata precyzji jest tak duża, że kolejne możliwe do zakodowania wartości zaczynają „przeskakiwać” co 2, potem co 4, 8 itd.

To zachowanie może być akceptowalne w wielu zastosowaniach, ale wyklucza stosowanie takich typów danych wszędzie tam, gdzie nawet dla dużych wartości musimy zachować precyzję co do jednostki, np. zapisując rozmiar pliku w bajtach czy ilość pieniędzy na koncie w banku. Okazuje się jednak, że dla 64-bitowego typu double maksymalna wartość całkowita reprezentowana dokładnie to już 253, a więc wystarczająco dużo do większości zastosowań - chyba, że pracujemy z plikami większymi, niż… 8 PB (petabajtów).

Skończona dokładność

Po drugie, skończona dokładność oznacza, że niektórych liczb nie jesteśmy w stanie zapisać precyzyjnie. Dla liczb takich, jak stałe matematyczne π i e czy ułamki typu 1/3 = 0.333333… jest to dla nas normalne, że posługujemy się ich przybliżeniami, ponieważ mają nieskończone rozwinięcie po przecinku. Okazuje się jednak, że pewne ułamki, które mają skończone rozwinięcie w systemie dziesiętnym, nie są możliwe do dokładnego zapisania w systemie binarnym. Taką liczbą jest np. 1/10 = 0.1. Listing 2 prezentuje prosty kod w C, który oblicza wartość 0.1 na dwa sposoby, a następnie je porównuje.

Listing 2. Problem dokładności liczb zmiennoprzecinkowych

Zła wersja ❌

float a = 1.0f / 10.0f;
float b = 1.0f - 0.9f;
printf("a=%g, b=%g\n", a, b);
// Wynik: a=0.1, b=0.1
if (a == b)
    printf("Są równe.\n");
else
    printf("Nie są równe!\n"); // <--

Dobra wersja ✔

float a = 1.0f / 10.0f;
float b = 1.0f - 0.9f;
printf("a=%g, b=%g\n", a, b);
// Wynik: a=0.1, b=0.1
if (fabsf(b - a) < 0.000001)
    printf("Są równe.\n"); // <--
else
    printf("Nie są równe!\n");

Okazuje się, że mimo wydrukowania obu zmiennych jako wartość 0.1, program ten równocześnie drukuje „Nie są równe!” Dzieje się tak dlatego, że te dwie wartości wyliczone na różne sposoby różnią się nieco wartością na dalszych miejscach po przecinku. Jeśli sprawdzimy ich wartość binarną oraz wydrukujemy je w systemie dziesiętnym, ale z większą dokładnością, to okaże się, że:

  • a = 0x3DCCCCCD ≈ 0.10000000149
  • b = 0x3DCCCCD0 ≈ 0.10000002384

Wniosek z tego przykładu płynie taki, że liczby zmiennoprzecinkowe faktycznie mają ograniczoną precyzję i nie powinniśmy polegać na ich dokładności co do ostatniego bitu. Dotyczy to szczególnie wyników obliczeń. Co do wartości stałych wpisywanych wprost do zmiennej, możemy jednak mieć pewność, że pozostaną identyczne, jak nasza stała, dopóki w jakiś sposób ich nie zmienimy.

Prawdą jest więc również, że wyników obliczeń nie należy porównywać operatorem ==, który zwróci prawdę wyłącznie, kiedy dwie liczby są identyczne co do wartości. Zamiast tego warto sprawdzać, czy liczby są dostatecznie blisko siebie, korzystając z wartości bezwzględnej (funkcja abs) ich różnicy i porównując ją z pewną małą wartością (nazywaną często ε – epsilon), jak pokazuje druga wersja kodu w Listingu 2. Ten program już prawidłowo wydrukuje: „Są równe.”

Jeszcze bardziej zaawansowanym sposobem porównywania dwóch liczb jest ten oferowany w Pythonie przez funkcję math.isclose, gdzie podajemy dwie wartości tolerancji: bezwzględną (abs_tol, od absolute tolerance) i względną (rel_tol, od relative tolerance), a wzór na obliczenie, czy liczby a i b są sobie mniej więcej równe, to:

abs(a - b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol)

Catastrophic cancellation

Po trzecie wreszcie, porównywanie lub odejmowanie od siebie dwóch liczb zmiennoprzecinkowych o dużej wartości może sprawić, że drobna różnica między nimi będzie odzwierciedlona niedokładnie lub całkowicie zginie. Zjawisko to nazywamy z ang. catastrophic cancellation.

Jako przykład, w Listingu 3 rozpatrzmy kod mierzący czas trwania jakiejś operacji. Na różnych platformach mamy zazwyczaj dostępną funkcję, która zwraca liczbę milisekund, nanosekund lub jakiś cykli od określonego momentu w czasie, na przykład od uruchomienia komputera albo od określonej daty. Takie funkcje często zwracają wynik jako liczbę całkowitą 64-bitową. Zmierzony czas trwania chcemy wyświetlić w sekundach, więc musimy także odpytać nasz system o częstotliwość tego zegara (liczbę cykli na sekundę) i podzielić przez tę częstotliwość. Ważne, żeby prędkość upływu tego czasu była stała, a nie zależny od aktualnej częstotliwości taktowania procesora. Wartości wyrażone w sekundach mogą być ułamkowe. Naturalne wydaje się więc użycie liczb zmiennoprzecinkowych.

Listing 3. Kod mierzący czas trwania operacji

Zła wersja ❌

LARGE_INTEGER freq, t;
QueryPerformanceFrequency(&freq);
QueryPerformanceCounter(&t);
double begin = (double)t.QuadPart / (double)freq.QuadPart;
LongOperation();
QueryPerformanceCounter(&t);
double end = (double)t.QuadPart / (double)freq.QuadPart;
double duration = end - begin;
printf("Duration in seconds: %g\n", duration);

Dobra wersja ✔

LARGE_INTEGER freq, begin, end;
QueryPerformanceFrequency(&freq);
QueryPerformanceCounter(&begin);
LongOperation();
QueryPerformanceCounter(&end);
double duration = (double)(end.QuadPart - begin.QuadPart) /
    (double)freq.QuadPart;
printf("Duration in seconds: %g\n", duration);

W pierwszej wersji kodu czai się jednak pewne niebezpieczeństwo. Jeżeli wartość zwrócona przez zegar będzie duża, a czas trwania zmierzonej operacji mały, to może się okazać, że po odjęciu tych dwóch wartości – jako liczby zmiennoprzecinkowe – otrzymana różnica reprezentująca czas trwania operacji będzie bardzo niedokładna lub wręcz za każdym razem wyniesie 0.

Rozwiązaniem jest odjęcie czasu początkowego od końcowego jeszcze jako liczb całkowitych, które są dokładne co do jednego cyklu, a następnie zamiana na liczbę zmiennoprzecinkową samej różnicy, która już będzie mniejsza. Takie podejście pozwala zachować dobrą precyzję obliczeń. Poprawiony kod przedstawia druga wersja w Listingu 3.

Czy te liczby są niedeterministyczne?

Inny popularny pogląd mówi, że liczby zmiennoprzecinkowe są niedeterministyczne, a więc nie powinniśmy się spodziewać za każdym razem takiego samego wyniku. Stosowanie takiej „zasady ograniczonego zaufania” również pomaga uniknąć wielu błędów, jednak warto poznać to zagadnienie bardziej szczegółowo. Okazuje się bowiem, że żaden generator liczb losowych tutaj nie działa za naszymi plecami. Ogólnie można więc powiedzieć, że ten sam program, uruchomiony na tej samej platformie, dla tych samych danych wejściowych powinien dać ten sam wynik. Istnieje jednak kilka nieoczywistych czynników, które mogą na ten wynik wpłynąć.

Po pierwsze, wyniki mogą być różne na różnych platformach. Standard IEEE 754 definiuje format liczb zmiennoprzecinkowych i operacji na nich, ale dla niektórych operacji nie gwarantuje zawsze identycznego wyniku. Podczas dodawania czy mnożenia wartości, o których wiemy, że mogą być zapisane precyzyjnie, możemy się spodziewać dokładnego wyniku. Jednak bardziej skomplikowane operacje, szczególnie funkcje niealgebraiczne (ang. transcendental), jak sinus i cosinus, na różnych platformach mogą dać wyniki różniące się najmniej znaczącą cyfrą (ULP - ang. unit in the last place), a więc nie identyczne co do reprezentacji bitowej i nie zwracające prawdy przy porównaniu operatorem ==. Wyniki mogą się różnić zależnie od producenta danego procesora czy układu graficznego (np. AMD, Intel, Nvidia) czy generacji danego chipa, co pozostaje nadal poprawne i zgodne ze standardem.

Po drugie, różne optymalizacje i inne przekształcenia wykonywane na naszym kodzie przez kompilator mogą zaowocować wykonywaniem innych instrukcji lub ich wykonywaniem w innej kolejności, co może wpłynąć na wynik. Tradycyjny zestaw instrukcji zmiennoprzecinkowych w architekturze x86, czyli FPU, wartości pośrednie obliczeń może przechowywać w wyższej, 80-bitowej precyzji. To może się wydawać korzystne, ale w praktyce grozi powstaniem różnych problemów, jeśli będziemy zbytnio polegali na określonym wyniku obliczeń.

Nowszych zestawów instrukcji zmiennoprzecinkowych może dotyczyć podobny problem. Na przykład wiele architektur CPU jak i GPU posiada instrukcję Fused Multiply-Add (FMA) lub Fused Multiply-Accumulate (FMAC), która wykonuje naraz 3-argumentową operację d = a * b + c. Logicznie są to dwie operacje, ale często są wykonywane w czasie jednej. Ma to ogromne znaczenie dla wydajności obliczeń, ponieważ z takich mnożeń i dodawań składa się operacja mnożenia macierzy, a ona z kolei leży u podstaw całej grafiki, jak i deep learningu.

Oprócz zwiększonej wydajności taka instrukcja często oblicza wynik dokładniej, niż dwie osobe operacje, ponieważ dokonuje zaokrąglenia wyniku tylko raz. To może spowodować, że dana funkcja matematyczna zwraca pewien wynik, kiedy kompilator widzi całe to wyrażenie i optymalizuje je do instrukcji FMA, a inny, kiedy inna część programu zawiera to samo wyrażenie wykonywane jako oddzielne mnożenie, a oddzielne dodawanie.

Z optymalizacjami wykonywanymi przez kompilator na obliczeniach liczb zmiennoprzecinkowych związana jest dostępna w wielu kompilatorach flaga nazywana Fast Math. Na przykład, spośród kompilatorów C++, w Visual Studio jest to parametr /fp:fast, w GCC i Clang: --ffast-math. Kiedy jest użyta, pozwala kompilatorowi robić dodatkowe przekształcenia, które – zgodnie z nazwą – mogą przyspieszyć kod. Są one tożsame z matematycznego punktu widzenia, ale mogą łamać standard IEEE 754 i nie gwarantują, że wynik takiego zoptymalizowanego wyrażenia będzie identyczny z oryginalnym.

Mogłoby się wydawać, że użycie tej flagi jest pożądane wszędzie tam, gdzie dokładność obliczeń nie jest kluczowa, a liczy się przede wszystkim wydajność – np. w grafice. Okazuje się jednak, że nawet pośród programistów gier i silników graficznych czy fizycznych, których powinno cechować bezkompromisowe dążenie do jak najwyższej wydajności kodu, panuje powszechny pogląd, że flagi Fast Math nie należy używać, bo choć przyspiesza ona w pewnym stopniu kod, to powoduje w zamian mnóstwo problemów. Więcej informacji na ten temat można znaleźć choćby w [2].

Drugą stroną medalu jest zablokowanie kompilatorowi możliwości optymalizowania naszych wyrażeń matematycznych. Kiedy nie włączamy flagi Fast Math, wówczas sygnalizujemy kompilatorowi, że operacje zmiennoprzecinkowe muszą się odbywać dokładnie tak, jak zapisaliśmy je w kodzie źródłowym, ponieważ nie chcemy, aby jakiekolwiek optymalizacje miały wpływ na wynik. To sprawia, że kompilator nie może wykonać nawet prostych przekształceń i uproszczeń, jakie wykonałby na kartce uczeń szkoły podstawowej, korzystając ze znanych tożsamości matematycznych. W związku z tym wyrażenia matematyczne powinniśmy upraszczać sami i do kodu wpisywać w postaci ręcznie zoptymalizowanej. Prostym przykładem może być zapisanie wielomianu d = ax2 + bx + c jak w Listingu 4:

Listing 4. Optymalizacja wielomianu

// Wersja najwolniejsza
d = a*pow(x, 2.0) + b*x + c;
// Wersja szybsza
d = a*x*x + b*x + c;
// Wersja najszybsza z prezentowanych
d = ((a*x) + b)*x + c;

Wreszcie, co może chyba najbardziej przerażające, jednostkę zmiennoprzecinkową w procesorze można konfigurować ustawiając jej różne tryby działania, które stanowią globalny stan dla danego wątku. One również mogą wpływać na wyniki obliczeń. Na przykład dla tradycyjnej jednostki FPU ustawia je funkcja _controlfp, a na działanie nowszych instrukcji wektorowych SSE wpływają flagi w rejestrze MXCSR ustawiane za pomocą funkcji _mm_setcsr. Może się więc okazać, że jakaś zewnętrzna biblioteka, której używamy, przestawia te tryby dla naszego programu.

Podsumowując, można chyba powiedzieć, że przy mnogości różnych czynników mogących wpływać na wyniki obliczeń, bezpieczniej jest założyć, że wyniki te są niedeterministyczne i nie powinniśmy się spodziewać dokładnego, konkretnego lub za każdym razem takiego samego wyniku. Te różnice nie biorą się jednak znikąd. Okazuje się, że rozumiejąc dogłębnie te zagadnienia oraz przygotowując starannie kod, można uzyskać pełen determinizm obliczeń na liczbach zmiennoprzecinkowych, nawet pisząc silnik fizyki obsługujący wiele platform i działający wielowątkowo, co udowodnili twórcy biblioteki Box2D [3].

Wartości specjalne

Wspomnieliśmy na początku, że wykładnik złożony z samych zer lub samych jedynek sygnalizuje wartości traktowane w specjalny sposób. Zajmijmy się nimi teraz nieco dokładniej. Rodzaje wartości specjalnych liczb zmiennoprzecinkowych przedstawiono w Tabeli 3. Te wartości specjalne obowiązują obojętnie, z jakim konkretnym typem zmiennoprzecinkowym mamy do czynienia, a więc ile bitów ma wykładnik czy mantysa. Nie zajmujemy się tutaj też bitem znaku, przyjmując, że każda z tych wartości może być zarówno dodatnia, jak i ujemna.

Tabela 3. Rodzaje wartości specjalnych liczb zmiennoprzecinkowych

WykładnikMantysaZnaczenieReprezentacja
000…0000…0Zero0
000…0Inna wartośćLiczba zdenormalizowanaWartość liczbowa
111…1000…0Nieskończonośćinf, Infinity
111…1Inna wartośćNot a NumberNaN, IND

O dwóch pierwszych już mówiliśmy. Kiedy zarówno wykładnik, jak i mantysa składają się z samych zer, jest to sposób na zakodowanie liczby zero. Mamy dwie takie wartości zależnie od bitu znaku, które oznaczają -0 i +0, ale ponieważ przy porównaniu one są sobie równe, nie musimy się tym martwić.

Liczby zdenormalizowane stanowią specjalny tryb, ponieważ w ich przypadku podczas dekodowania wartości nie obowiązuje dopisanie domyślnej jedynki przed przecinkiem. Dzięki temu możemy kodować wartości jeszcze mniejsze, niż pozwalają na to wartości normalne (znormalizowane). To również jest funkcja obsługiwana automatycznie i nie musimy się nią zajmować.

Kiedy wykładnik składa się z samych jedynek, a mantysa z samych zer, jest to sposób na zakodowanie wartości nieskończoności ∞. Zależnie od bitu znaku może to być plus lub minus nieskończoność. Po zamianie na łańcuch znaków przedstawiana jest zazwyczaj jako „inf” lub „Infinity”. Taka wartość wychodzi z obliczeń, kiedy nie mieści się już w zakresie możliwych do zapisania liczb, np. kiedy wartość jest już bardzo duża, a my nadal zwiększamy ją, dodając, mnożąc przez inną dużą liczbę lub (co bywa chyba najczęstszą przyczyną pojawienia się nieskończoności) dzieląc przez liczbę bardzo małą (bliską zeru) lub przez zero.

Ostatni rodzaj wartości specjalnych to każda taka wartość, która ma w wykładniku same jedynki, a w mantysie wartość różną od zera. Oznacza ona Not a Number, a więc wartość nieprawidłową, błąd obliczeń. Po zamianie na string przedstawiana bywa jako „NaN” lub radziej „IND” (od ang. indeterminate – nieokreślony).

Cechą charakterystyczną wartości NaN jest jej propagacja przez kolejne wyrażenia. Prawie każda operacja otrzymująca na wejściu NaN zwraca w wyniku NaN. Powoduje to, że kiedy taka nieprawidłowa wartość dostanie się w „tryby” obliczeń, wyniki już do końca są nieprawidłowe. Taką sytuację zilustrowano na Rysunku 3. Zadaniem programisty jest wtedy prześledzić, gdzie w danych wejściowych lub wyrażeniach pośrednich po raz pierwszy pojawiają się błędne wartości. Często przyczyną bywa wartość ujemna, zerowa lub bardzo mała (bliska zeru) tam, gdzie się tego nie spodziewamy lub odczytywanie przypadkowych danych z pamięci albo z pliku i traktowanie ich jako liczby zmiennoprzecinkowe.

Rysunek 3. Propagacja wartości NaN w obliczeniach

Osobliwe jest też zachowanie wartości NaN podczas porównania. Porównanie NaN z dowolną inną liczbą zawsze zwraca false – nawet z samą sobą! Dzięki temu warunek if(a == a) to najprostszy sposób na sprawdzenie, czy zmienna przenosi wartość liczbową, a nie NaN. Biblioteki standardowe w językach programowania oferują jednak dedykowane funkcje do takiego sprawdzania. Na przykład w C nagłówek <math.h> definiuje funkcje: isinf, isnan, a także isfinite, która zwraca prawdę, kiedy podany argument nie jest ani nieskończonością, ani NaN.

Często przyjmuje się, że pojawienie się INF lub NaN w wynikach obliczeń sygnalizuje błąd, który należy znaleźć i naprawić. Najczęściej tak właśnie jest, jednakże te wartości specjalne powstają i zachowują się w obliczeniach w pewien logiczny sposób, zgodny z zasadami matematyki. Niektóre z tych zasad ilustrują przykładowe wyrażenia i ich wyniki przedstawione w Listingu 5.

Listing 5. Zachowanie wartości INF i NaN

 2 / 0    ==  inf
-2 / 0    == -inf
 0 / 0    ==  NaN
log ( 0)  == -inf
log (-2)  ==  -NaN
sqrt(-2)  ==  -NaN
 2 + inf  ==  inf
 2 * inf  ==  inf
-2 * inf  == -inf
 0 * inf  ==  NaN
 2 / inf  ==  0
 0 / inf  ==  0
inf + inf ==  inf
inf - inf ==  -NaN
inf * inf ==  inf
inf / inf ==  -NaN
 2 + NaN  ==  NaN
 2 * NaN  ==  NaN
 2 / NaN  ==  NaN
 inf > 2  ==  true
-inf < 2  ==  true
(  2 == NaN) == false
(NaN == NaN) == false

Na przykład warto zauważyć, że:

  • Dzielenie przez zero zwraca ∞, a dzielenie liczby ujemnej przez zero: -∞.
  • Logarytm z zera zwraca -∞, ale już logarytm z liczby ujemnej, jak również pierwiastek z liczby ujemnej daje w wyniku NaN.
  • Próba dalszego zwiększenia ∞ przez dodanie czy pomnożenie jej przez jakąś liczbę zwraca nadal ∞, ale pomnożenie jej przez liczbę ujemną zwraca -∞.
  • Podczas porównania, -∞ jest mniejsza, a +∞ jest większa od każdej skończonej liczby.

Znając te zasady, możemy spróbować świadomie wykorzystać wartości specjalne w programie. W Listingu 6 pokazano napisaną w C++ funkcję, która znajduje wartość największego elementu z podanej tablicy liczb typu float. Zmienna tymczasowa przechowuje największą znalezioną dotychczas wartość. Jest ona na początku inicjalizowana wartością -∞, co gwarantuje, że każda skończona liczba będzie od niej większa. Taka implementacja algorytmu ma tę ciekawą właściwość, że gwarantuje określone zachowanie po podaniu tablicy pustej (count == 0) - zwraca wtedy -∞. Zachowanie to można by opisać w dokumentacji tej funkcji.

Listing 6. Znajdowanie maksimum w tablicy

float FindMaxValue(const float* arr, size_t count)
{
    float maxVal = -std::numeric_limits<float>::infinity();
    for(size_t i = 0; i < count; ++i)
        if(arr[i] > maxVal)
            maxVal = arr[i];
    return maxVal;
}

Na koniec warto jeszcze raz wrócić do problemu dzielenia przez zero, aby podkreślić jedną rzecz. Dzięki temu, że dzielenie liczb zmiennoprzecinkowych przez zero zwraca nieskończoność, mamy szansę wychwycić ten błąd w programie i jakoś go obsłużyć. Inaczej jest z liczbami całkowitymi, gdzie próba dzielenia przez zero stanowi błąd krytyczny i przerywa cały program, więc nie można dopuścić, aby takie dzielenie się wykonało.

Mniejsze formaty

W przykładach w tym artykule posługiwaliśmy się liczbami 32-bitowymi pojedynczej precyzji, omawiając ich ograniczenia. Wspomnieliśmy też o liczbach 64-bitowych podwójnej precyzji. Do pewnych zastosowań nadają się jednak liczby jeszcze mniejsze. Istnieje na przykład format 16-bitowy „połówkowej precyzji” (ang. half-float, nazywany też fp16), który przeznacza 5 bitów na wykładnik i 10 bitów na mantysę.

Taki typ danych ma oczywiście bardzo ograniczoną precyzję (do około 3 znaczących cyfr dziesiętnych) oraz zakres (maksymalna wartość to 65 504, po przekroczeniu której pojawia się już nieskończoność). Typ ten bywa wykorzystywany w grafice i obsługiwany sprzętowo przez niektóre układy GPU. Jeżeli weźmiemy pod uwagę, że komponenty RGB kolorów pikseli tradycyjnie zapisywane są za pomocą 8 bitów o wartości 0…255, to zakres i precyzja liczb fp16 okazują się aż nadto wystarczające do reprezentowania kolorów, nawet na monitorach z obsługą HDR.

Używanie mniejszych formatów danych tam, gdzie ich możliwości są wystarczające, pozwala oszczędzić pamięć, a także przyspieszyć obliczenia. Na przykład jednostka obliczeniowa potrafiąca w jednym cyklu wykonać operację na wektorze 4 liczb zmiennoprzecinkowych 32-bitowych może być w stanie w takim samym czasie i wykorzystując te same tranzystory, wykonać operację na 8 liczbach 16-bitowych upakowanych w rejestrze o tej samej szerokości. Oszczędność pamięci ma też niebagatelne znaczenie dla ilości przesyłanych danych, bo we współczesnych algorytmach to często przepustowość odczytów i zapisów do pamięci (ang. bandwidth) ogranicza wydajność programu, a nie szybkość samych obliczeń.

Mające ostatnio ogromne znaczenie w informatyce algorytmy sztucznej inteligencji, a konkretniej deep learning i związane z nim sieci neuronowe, mają charakter obliczeń jeszcze bardziej regularnych i wykonywanych na jeszcze większej ilości danych, niż to ma miejsce w grafice. Badacze tej dziedziny wykazali w pracach naukowych, że do różnego rodzaju modeli sieci neuronowych wystarczające są obliczenia na liczbach jeszcze niższej precyzji, przy czym do niektórych operacji (konkretnie do treningu – zapisywania gradientów) ważniejszy jest większy zakres, niż precyzja. Zaproponowano więc drugi format 16-bitowy nazywany bf16 (od ang. brain float), który przeznacza 8 bitów na wykładnik i 7 bitów na mantysę.

Czasami w obliczeniach AI liczby bywają kwantyzowane do typów całkowitych (np. int8, int4, a nawet liczb… 1-bitowych). Zaproponowano jednak również formaty zmiennoprzecinkowe o szerokości 8 bitów nazywane fp8 i bf8, mające odpowiednio 4 bity wykładnika i 3 bity mantysy lub 5 bitów wykładnika i tylko 2 bity mantysy, co bywa też oznaczane jako E4M3 i E5M2. Więcej informacji na ich temat można znaleźć w [4]. Wsparcie sprzętowe dla takich typów nie jest jednak póki co zbyt rozpowszechnione za wyjątkiem chipów przeznaczonych typowo do obliczeń sztucznej inteligencji.

Podsumowanie

Liczby zmiennoprzecinkowe są używane w programowaniu wszędzie tam, gdzie potrzebujemy zapisywać nie tylko wartości całkowite, ale i ułamki, liczby bardzo małe lub bardzo duże. Takie typy danych nazywają się różnie w różnych językach programowania, ale najczęściej sprowadzają się do tych samych formatów zgodnych ze standardem IEEE 754, a obliczenia na nich są wspierane sprzętowo przez procesory i układy graficzne.

Zrozumienie budowy i działania liczb zmiennoprzecinkowych, ich możliwości i ograniczeń pozwala wykorzystywać je w sposób świadomy i uniknąć różnych błędów, które czyhają na początkujących czy nieuważnych programistów. W tym artykule omówiliśmy wiele związanych z nimi zagadnień – problem ograniczonej precyzji, niedeterministycznych wyników, obsługiwanych wartości specjalnych (INF, NaN) i ich zachowanie.

Temat liczb zmiennoprzecinkowych można zgłębiać dalej. Wiele jest zagadnień, których nie omówiliśmy lub omówiliśmy tylko pobieżnie. W ogóle nie rozważaliśmy wydajności obliczeń. Warto wiedzieć, że nawet mając w sprzęcie dostępne instrukcje do obliczania różnych funkcji, niektóre z nich działają szybko, a inne wolniej. Najszybsze jest zwykle dodawanie, odejmowanie i mnożenie, a funkcje niealgebraiczne (ang. transcendental) jak sinus, cosinus, potęga, pierwiastek, logarytm, a także… dzielenie mogą potrzebować więcej cykli na wykonanie. Nie powiedzieliśmy też nic o zestawach instrukcji wektorowych w procesorach - MMX, SSE, AVX, które pozwalają wykonać tę samą operację za jednym razem na wielu liczbach, znacznie przyspieszając obliczenia.

Także liczby zdenormalizowane na niektórych platformach działają wolniej, niż normalne. Bywają też dostępne tryby, które całkowicie wyłączają obsługę liczb zdenormalizowanych, sprowadzając je zawsze do zera (co bywa nazywane flush to zero). Podobnie, niektóre platformy obsługują różne tryby zaokrągleń - np. w dół (w kierunku -∞), w górę (w kierunku +∞), w kierunku zera lub w kierunku najbliższej liczby parzystej. Nie wspomnieliśmy też o rozróżnieniu między quiet NaN a signaling NaN.

Bibliografia

Adam Sawicki
Listopad 2024

Comments

[Download] [Dropbox] [pub] [Mirror] [Privacy policy]
Copyright © 2004-2025