Problem z operacjami arytmetycznymi realizowanymi przez komputery polega na tym, że „po cichu" zakładamy, że argumenty równe są swoim reprezentacjom zmiennoprzecinkowym, oraz, że podczas obliczeń, nie wystąpi ani nadmiar1 ani niedomiar2.
Niech $a=\operatorname{rd}(a)=2^{c_a}m_a$ i $b=\operatorname{rd}(b)=2^{c_b}m_b$ będą argumentami. Przyjmijmy, że $a\ge b>0$. Ich (dokładna) suma wynosi:
$$a+b=2^{c_a}(ma_a+m_b2^{-(c_a-c_b)})=2^{c_a}m_s$$Mantysa wyniku wyliczana jest przez dodanie do mantysy liczby $a$ odpowiednio przekształconej (to znaczy tak, żeby cechy liczb $a$ i $b$ były jednakowe) mantysy liczby $b$.
Wszystko to wygląda jakoś tak:
Jakość wyniku zależy od możliwości poprawnego zaokrąglenia wyniku. A zatem od liczby bitów, na których zapisywany jest wynik. Gdy jest ich tylko $t$ — nie jest dobrze. Historycznie, liczba dodatkowych bitów wewnętrznego rejestru arytmometru zmieniała się od zera aż do $t$ (wynikowy rejestr podwójnej długości). Dopiero norma IEEE-754 [1] wprowadziła „obowiązek" korzystania z dodatkowego bitu.
Jeżeli symbolem $\operatorname{fl}$ oznaczać będziemy realizację działania w arytmetyce zmiennoprzecinkowej, to:
$$\operatorname{fl}(a+b)=\operatorname{rd}(a+b)$$Ogólnie, dla dowolnej operacji $\diamond$
$$\operatorname{fl}(a \diamond b) = \operatorname{rd}(a\diamond b)= (a \diamond b)(1+\varepsilon),\quad \varepsilon=\varepsilon(a,b,\diamond), \quad |\varepsilon|\le 2^{-t}$$uzyskany (komputerowo) wynik, nieznacznie, różni się od wyniku „prawdziwego".
Wydaje się, że niewielkie ($2^{-t} \in [10^{-15}, 10^{-6}]$) względne błędy nie powinny mieć wielkiego wpływu na wyniki. Okazuje się, że jest inaczej, a sposób realizacji obliczeń może mieć ogromny wpływ na wyniki.
Rozpatrzmy dwa różne algorytmy wyliczania wartości różnicy kwadratów dwu liczb:
$A1(a^2-b^2) = a^2-b^2$
$A2(a^2-b^2) = (a-b)(a+b)$
W przypadku gdy $a^2$ jest bliskie $b^2$, a $\varepsilon_1$ i $\varepsilon_2$ mają przeciwne znaki — błąd względny $\delta$ moze być dowolnie duży.
W przypadku drugiego algorytmu:
$$\begin{split} \operatorname{fl}((a-b)(a+b)) & = ((a-b)(1+\varepsilon_1)(a+b)(1+\varepsilon_2))(1+\varepsilon_3)\\ & = (a^2-b^2)(1+\delta) \end{split}$$gdzie $\delta \le |\varepsilon_1| + |\varepsilon_2| + |\varepsilon_3|$, zatem $\delta$ jest zawsze mniejszy od $3\cdot 2^{-t}$.
W szczególności, ze względu na algorytm dodawania, gdy $a$ i $b$ są takie, że $|b|\le \dfrac{1}{2}2^{-t}|a|$, zachodzi zależność:
$$\operatorname{fl}(a+b) \equiv a.$$
(Wspominałam o tym na zajęciach z Technologii Informacyjnych.)
Konsekwencje tego mogą być bardzo poważne w przypadku numerycznego aproksymowania ilorazem różnicowym, pochodnej funkcji w punkcie. Zamiast otrzymywania coraz to lepszego przybliżenia pochodnej (gdy różnica wartości pomiędzy dwiema wartościami zmiennej niezależnej) dostaniemy wartości nie mające nic wspólnego z pochodną.
Cytowania
IEEE Standard for Floating-Point Arithmetic,
IEEE, 2008. doi:10.1109/IEEESTD.2008.4610935