Do arytmetyki zmiennoprzecinkowej

Przykład. Załóżmy, że wartości funkcji
zostały obliczone na komputerze. Może się zdarzyć, że wykres funkcji
ma mnóstwo różnych miejsc zerowych w okolicy
, a sama funkcja nie jest gładka.

Tymczasem, ponieważ także
, to funkcja
ma dokładnie jedno miejsce zerowe (krotność pierwiastka - 4).

Wykonywanie realnych obliczeń na liczbach rzeczywistych w komputerze może być źródłem wielu innych zaskoczeń.

Przykład. W komputerze
co można łatwo sprawdzić:

octave:7> 10 * (1.1 -1) - 1

ans = 8.8818e-16

Jaki wynik dostaniemy kalkulatorze?

Dlatego w praktyce numerycznej należy wystrzegać się testów w rodzaju

if (x == 1.0)

{

....

}

Źródło problemów leży w zbyt uproszczonym modelu obliczeniowym. Matematycznym modelem arytmetyki maszyny cyfrowej jest arytmetyka zmiennoprzecinkowa
.

Niech będzie zadana liczba naturalna
(jej znaczenie wyjaśni się dalej). Dowolną liczbę rzeczywistą
można jednoznacznie przedstawić w postaci

,

gdzie
jest znakiem, liczba całkowita
cechą, a liczba rzeczywista
mantysą liczby
.

Zauważmy, że taki rozkład jest jednoznaczny i odpowiada przesuwaniu przecinka w rozwinięciu binarnym liczby do pierwszej cyfry znaczącej, tj. różnej od zera (stąd nazwa: reprezentacja zmiennoprzecinkowa (floating point)).

Mantysa ma w ogólności nieskończenie wiele cyfr binarnych
w swoim rozwinięciu dwójkowym

.

W komputerach osobistych mamy do czynienia z reprezentacją liczb rzeczywistych, w której do zapisania liczby używa się ściśle określonej liczby bitów
do zapisania mantysy i także określonej liczby bitów
do zapisania cechy danej liczby niezerowej

.

Liczby zapisane przy użyciu powyższej sekwencji bitów nazywa się liczbami maszynowymi. Są to jedyne dokładnie zapisywalne w komputerze liczby rzeczywiste, pozostałe będą musiały zostać wyrażone z wykorzystaniem liczb maszynowych.

Reprezentacją zmiennoprzecinkową niezerowej liczby będziemy nazywać liczbę
taką, że

,

gdzie
jest liczbą dwójkową postaci
, natomiast
jest liczbą naturalną postaci
.

Na znak liczby,
, przeznaczony jest jeden bit. Wartości
i
dobiera się tak, żeby
była tak bliska
jak to możliwe. Stałą całkowitą
dobiera się tak, by uzyskać zbalansowany zakres cechy
(mniej więcej tyle samo wartości ujemnych i dodatnich), a zysk z korzystania z niej jest taki, że nie marnujemy dodatkowego bitu na przechowywanie znaku wykładnika potęgi dwójki
.

Przy takim sposobie reprezentacji, jej błąd względny szacuje się przez

,

gdzie liczbę
nazywa się precyzją arytmetyki. Precyzja arytmetyki zależy wyłącznie od liczby bitów przeznaczonych na reprezentację mantysy. Ostatnią nierówność wygodnie jest zapisać w równoważny sposób jako

.

Przykład. Rozważmy system, w którym zarówno na cechę, jak i mantysę, przeznaczone są jedynie po dwa bity, zatem jedna liczba maszynowa zajmuje 5 bitów. Ponieważ możliwy zakres wartości
jest
, to przyjmiemy korektę
, dzięki czemu
. Z kolei możliwe wartości mantysy to

Wobec tego, jedyne (dodatnie) liczby maszynowe naszej pięciobitowej arytmetyki zmiennopozycyjnej to

Liczby maszynowe: reprezentowane dokładnie w pięciobitowej arytmetyce o precyzji
. (Przedstawiliśmy tylko liczby nieujemne)

Standard IEEE 754. Z nielicznymi egzotycznymi wyjątkami (np. Cray C90), współczesne procesory implementują IEEE 754 Floating Point Standard, który definiuje dwa zasadnicze formaty reprezentacji zmiennoprzecinkowej liczb rzeczywistych:

Typ IEEE 754	Pojedynczej precyzji	Podwójnej precyzji
Nazwa typu w C	float	double
Liczba bitów cechy	8	11
Liczba bitów mantysy	23	52
Liczba bajtów dla typu w C	4	8
Bias (liczba powyżej)	127	1023
Orientacyjny zakres
Orientacyjna precyzja

Procesory Intel'a mają jedną z najlepszych implementacji standardu IEEE 754).

W Octave można łatwo podejrzeć reprezentację binarną liczby zmiennoprzecinkowej podwójnej precyzji (jest to domyślny typ numeryczny stosowany w MATLABie i Octave),

octave:9> format bit

octave:10> x = -2

x = 1100000000000000000000000000000000000000000000000000000000000000

octave:11> x = 1/4

x = 0011111111010000000000000000000000000000000000000000000000000000

octave:13> x = 0

x = 0000000000000000000000000000000000000000000000000000000000000000

octave:15> x = 0.1

x = 0011111110111001100110011001100110011001100110011001100110011010

(w MATLABie możemy zobaczyć tę samą liczbę w zapisie szestnastkowym).

Rozwinięcie dwójkowe liczby 0.1 jest nieskończone:

Ten banalny fakt jest bardzo często przeoczony przez programistów, a w 1991 roku doprowadził do awarii systemu antyrakietowego Patriot. Okazało się, że rakiety Patriot traciły skuteczność, gdy przez wiele godzin pozostawały w stanie gotowości. Jak zbadano, w celu pomiaru czasu, zliczano kolejne tyknięcia zegara rakiety, które następowały dokładnie co 0.1 sekundy. Następnie, w celu wyznaczenia prawdziwego czasu, mnożono liczbę tknięć zegara przez 0.1 (które właśnie było niedokładnie reprezentowane). Gdy cykli zegara było bardzo dużo, błąd bezwzględny wyznaczenia czasu stawał się na tyle poważny, że uniemożliwiał precyzyjne wyznaczenie parametrów toru lotu nieprzyjacielskiego obiektu.

Nie wszystkie maszyny liczące wykorzystują reprezentację dwójkową.

Nadmiar i niedomiar. W pierwszym przypadku liczba jest tak duża (co do modułu), że nie zawiera się w przedziale liczb reprezentowalnych, a w drugim jest tak mała, że musi być reprezentowana przez zero.

Np. próżnia wokół zera (na przykładzie 5-bitowej arytmetyki)

Arytmetyka IEEE 754 przyjmuje, że liczby dla których następuje overflow są reprezentowane przez specjalną wartość Inf (nieskończoność, ze znakiem), która propaguje się w obliczeniach zgodnie z powszechnie przyjętymi regułami, np. 1+Inf daje Inf, 1/Inf daje 0, Inf-Inf daje NaN, itd.

Np. wszystkie liczby większe od największej zapisywalnej liczby są reprezentowane przez Inf (na przykładzie 5-bitowej arytmetyki)

W dalszych rozważaniach zjawiska nadmiaru i niedomiaru będziemy dla uproszczenia zaniedbywać, jednak nie zawsze jest to uzasadnione, o czym niech świadczy poniższy przykład.

Przykład. Jedną z najczęściej wykonywanych operacji na wektorze
jest obliczenie jego normy euklidesowej

Jak widać, możemy tu łatwo zetknąć się ze zjawiskiem zarówno niedomiaru, jak i nadmiaru, gdyż może się na przykład tak złożyć, że mimo iż
jest reprezentowana, to
już nie (np. w arytmetyce podwójnej precyzji
i
). Łatwym wyjściem z tej sytuacji jest wstępna normalizacja danych tak, by wszystkie nie były większe od 1: niech
i wtedy

i teraz suma pod pierwiastkiem jest zawsze pomiędzy 1 a
.

Liczby denormalizowane. Wymaganie, że mantysa jest postaci
,
, powoduje, że wokół zera pojawia się coś w rodzaju próżni. Formalnie, liczby mniejsze niż
powinny być reprezentowane przez 0, lecz zazwyczaj zamiast tego

octave:16> format bit

octave:17> x = 2^(-1022)

x = 0000000000010000000000000000000000000000000000000000000000000000

octave:18> x = 2^(-1023)

x = 0000000000001000000000000000000000000000000000000000000000000000

octave:19> x = 2^(-1028)

x = 0000000000000000010000000000000000000000000000000000000000000000

W ten sposób można (w podwójnej precyzji) zbliżyć się do zera na odległość około
.

Liczby denormalizowane trochę wypełniają próżnię wokół zera

Działania w arytmetyce zmiennoprzecinkowej. Standard IEEE 754 określa także "prawidłowy" sposób realizacji działań arytmetycznych w arytmetyce zmiennoprzecinkowej. W arytmetyce
implementującej standard IEEE 754, działania arytmetyczne na liczbach rzeczywistych (a raczej na ich reprezentacjach) muszą być implementowane tak, jakby działanie było wykonywane dokładnie i tylko wynik reprezentowano w zbiorze liczb maszynowych. Mamy więc

,

gdzie
, Ogólniej, jeśli
i
są wyrażeniami o wartościach rzeczywistych, to dla dowolnych wartości zmiennych

,

Zwykle dla prostoty będziemy również zakładać podobną zależność dla niektórych funkcji standardowych, o ile należą one do zbioru operacji elementarnych (chociaż w rzeczywistości są one obliczane przez procedury używające czterech podstawowych operacji arytmetycznych). I tak będziemy mieć np.

,

,

gdzie
i
są "niewielkimi" stałymi.

Przypuśćmy, że w naszym prościutkim pięciobitowym systemie spełniającym wymogi standardu IEEE 754 zechcemy wykonać mnożenie

.

Poniżej możemy przekonać się, jak będzie ono przebiegać (na przykładzie 5-bitowej arytmetyki).

Liczba 1.3 nie jest dokładnie reprezentowalna w naszym systemie

Jej reprezentacja to najbliższa jej liczba maszynowa 1.25

Również drugi czynnik, 2.4, nie jest liczbą maszynową

A więc jego reprezentacją będzie znów najbliższa mu liczba maszynowa.

Mnożenie odbywa się już na reprezentacjach obu czynników

Wynik dokładnego mnożenia tych liczb maszynowych to 3.125 i znowu musi być zaokrąglony ... do najbliższej liczby maszynowej.

Ostatecznie, błąd względny wyniku wynosi około
i jest znacznie mniejszy niż pesymistyczne oszacowanie (czasem osiągalne, lecz nie tym razem)
.

Podobnie, jeśli
, to wartością wyrażenia logicznego
w
jest dokładna wartość wyrażenia
.

Dla specjalnej liczby NaN (not-a-number), która pojawia się jako wynik zabronionych operacji matematycznych, np.
,
, Inf - Inf, itp., zawsze zachodzi, że NaN
NaN.

Standard IEEE 754 nie gwarantuje, że działania arytmetyczne będą łączne, co widać na poniższym przykładzie:

octave:9> 7.1 - (7+0.1)

ans = 0

octave:10> (7.1 - 7) - 0.1

ans = -3.6082e-16

Praktyczne wyznaczanie precyzji arytmetyki. Aby wyznaczyć precyzję używanej przez nas arytmetyki możemy wykonać prosty test. Pomyślmy, jaka jest najmniejsza dodatnia liczba
, która dodana do jedności da w wyniku liczbę większą od 1.0 (liczbę
nazywa się epsilonem maszynowym, macheps). Jasne jest, że w przypadku arytmetyki IEEE 754 jest to liczba równa podwojonej precyzji arytmetyki,
, gdzie
jest liczbą cyfr mantysy
. Stąd dostajemy prosty algorytm wyznaczania epsilona maszynowego:

x = 1.0;

while ( 1.0 + x > 1.0 )

{

x = x / 2.0;

}

printf("Macheps = %g", 2.0*x);

}

Jednak, w rzeczywistości musimy być bardziej ostrożni. Implementując ten algorytm w C następująco

/* Wyznaczanie epsilona maszynowego, wersja 1 */

#include <stdio.h>

int main(void)

{

int dt;

double dx;

dt = 0; dx = 1.0;

while(1.0 + dx > 1.0)

{

dx *= 0.5;

dt++;

}

printf("Macheps (double) = %g. Liczba bitów mantysy = %d\n", 2*dx, dt);

return(0);

}

dostajemy wynik niezgodny z oczekiwaniami:

Macheps = 1.0842e-19. Liczba bitów mantysy = 64.

Wynika to stąd, że w C obliczenia wykonują się zawsze z maksymalną możliwą precyzją. W procesorach x86 jest to precyzja arytmetyki extended double precision, wykorzystującej 80 bitów do reprezentacji liczb. Dlatego działanie

1.0 + dx > 1.0

wykona się w arytmetyce nie podwójnej (64-bitowej), ale rozszerzonej podwójnej precyzji. Aby sprawić, by działanie zostało wykonane z wykorzystaniem typu double, musimy nasz program trochę zmodyfikować:

/* Wyznaczanie epsilona maszynowego, wersja 2 */

#include <stdio.h>

int main(void)

{

int dt;

double dx, dxp1;

dt = 0; dx = 1.0; dxp1 = 2.0;

while(dxp1 > 1.0)

{

dx *= 0.5;

dxp1 = 1.0 + dx; /* tym razem wynik działania zostanie zapisany

do zmiennej typu double */

dt++;

}

printf("Macheps (double) = %g. Liczba bitów mantysy = %d\n", 2*dx, dt);

}

Tym razem wynik jest prawidłowy:

Macheps = 2.22045e-16. Liczba bitów mantysy = 53

Sprawdzić, jak zmienią się wyniki, gdy wykorzystuje się opcje kompilacji:

• gcc -O3

• gcc -ffast-math

• gcc -O3 -ffast-math

Objaśnić te wyniki, wspomagając się ewentualnie dokumentacją kompilatora.

Wpływ błędu zaokrągleń na wyniki obliczeń. Redukcja cyfr.

Korzystając z wprowadzonego powyżej modelu arytmetyki zmiennoprzecinkowej możemy spróbować ocenić wpływ błędu reprezentacji i błędu zaokrągleń na wynik konkretnego algorytmu.

Rozważmy zadanie wyznaczenia iloczynu
liczb z tablicy

W tym celu stosujemy banalny algorytm:

s = 1.0;

for (i=0; i < N; i++)

s *= x[i];

Sprawdźmy, jak będzie on realizowany w arytmetyce zmiennoprzecinkowej. Dla uproszczenia założymy, że nie wystąpiło ani zjawisko nadmiaru, ani zjawisko niedomiaru (w przeciwnym razie dostaniemy w wyniku, odpowiednio,
Inf lub 0).

Zamiast dokładnych wartości
, będziemy mieli w komputerze jedynie ich reprezentacje,
, przy czym
.

Oznaczając
wyznaczoną numerycznie wartość iloczynu po
-tym kroku pętli, mamy, że

gdzie znów
. Ostatecznie więc, wyznaczona wartość iloczynu,
spełnia

Ponieważ
, gdzie, z pominięciem małych wyższego rzędu,
, dostajemy ostatecznie

gdzie
. Ponieważ w arytmetyce podwójnej precyzji
, to nawet biorąc iloczyn tysiąca liczb, dostajemy wynik, którego błąd względny będzie (o ile tylko nie wystąpi nadmiar/niedomiar) bardzo mały, rzędu
.

Trzeba wyznaczyć różnicę dwóch liczb:

Prowadząc analizę jak poprzednio, widzimy, że obliczona numerycznie wartość to

Stąd po prostych oszacowaniach

A więc, gdy
, to
i w efekcie możemy utracić nawet wszystkie znaczące cyfry wyniku! To zjawisko nosi żargonową nazwę utraty cyfr przy odejmowaniu, choć precyzyjnie powinno się mówić o "zmniejszeniu liczby dokładnych cyfr znaczących wyniku przy odejmowaniu dwóch bardzo bliskich sobie liczb".

Zauważmy, że prowadząc identyczną analizę dla sumy dwóch liczb
, gdzie
i
są tego samego znaku, dostajemy oszacowanie błędu względnego równe
, niezależnie od wartości liczbowych
i
.

Niech
. Korzystając ze szkolnego wzoru na pierwiastki równania kwadratowego

gdzie
, możemy natknąć się na trudności, gdy jeden z pierwiastków jest bardzo bliski zera (tzn. gdy
).

Niestety, skoro
, to wyznaczając mniejszy pierwiastek
ryzykujemy utratę cyfr przy odejmowaniu. Ale na szczęście, można to ominąć zgodnie z często stosowaną w numeryce regułą:

Jeśli problem jest trudny --- najlepiej zmienić problem!

W naszym wypadku ratunkiem jest matematyczna transformacja problemu tak, by już nie było w nim odejmowania bliskich sobie liczb. Rzeczywiście, przecież wciąż mamy dobry wzór na większy z pierwiastków,
! Dokładając do tego wzór Viete'a,

dostajemy inny wzór na
, nie zawierający odejmowania.

O pakietach obliczeń symbolicznych

Czasem możemy spotkać się z twierdzeniem, że nie warto zawracać sobie głowy metodami numerycznymi, gdyż są dostępne pakiety obliczeń symbolicznych (Maple, Mathematica, MuPAD, Maxima), które potrafią "wszystko" "policzyć z dowolną precyzją".

Precyzja, o której mowa, jest jedynie precyzją używanej arytmetyki (rzeczywiście, softwarowo można emulować dowolną precyzję), ale dokładność wyniku nie może być w nich a priori zadana. Wiele bowiem może zależeć od właściwości samego zadania obliczeniowego.

Przykład. Precyzja w pakietach symbolicznych. Oto fragment sesji pakietu obliczeń symbolicznych MUPAD:

>> ((4/3)*3 - 3) - 1

0

>> DIGITS := 10

10

>> ((4/3.0)*3 - 3) - 1

-2.168404345e-19

>> subs(((4/a)*3 - 3) - 1, a = 3.0)

-4.33680869e-19

>> subs(((4/a)*3 - 3) - 1, a = 3)

0

Najpierw wyznaczona została wartość wyrażenia algebraicznego, przy użyciu rachunków symbolicznych -wyrażenie upraszcza się do zera.

Następnie zażądaliśmy, by DIGITS (parametr sterujący: "liczbą cyfr znaczących dla liczb zmiennoprzecinkowych") przyjął wartość równą 10.

Dalej, wymuszając (przez wpisanie 3.0, zamiast 3) stosowanie w obliczeniach arytmetyki zmiennoprzecinkowej w miejsce rachunków symbolicznych dostajemy wynik, który nie ma ani jednej cyfry znaczącej dokładnej. Z drugiej strony, widzimy także, iż faktycznie stosowana precyzja obliczeń jest znacznie większa i wynosi około 19 cyfr znaczących, chociaż żądaliśmy jedynie 10...

Jak wynika z powyższego, w praktyce pakiety symboliczne stosują znacznie większą niż żądana precyzję obliczeń, by ustrzec się najbardziej typowych patologii.

---