Wersja z 18:11, 1 wrz 2006

Arytmetyka zmiennoprzecinkowa

Metody iteracyjne mają czasem kłopoty, które nie są związane z samą naturą problemu matematycznego. Przyrzyjmy się bowiem, jak w dużym zbliżeniu wygląda wykres funkcji $w (x) = x^{4} - 4 x^{3} + 6 x^{2} - 4 x + 1$ , której wartości zostały obliczone na komputerze PC. Nietrudno sprawdzić, że $w$ ma dokładnie jedno miejsce zerowe, gdyż $w (x) = (x - 1)^{4}$ . Tymczasem, wykres $w$ (wyznaczony oryginalnym wzorem) zdaje się mieć mnóstwo różnych miejsc zerowych w okolicy $x = 1$ . Co gorsza, wygląda na to, że $w$ wcale nie jest gładka!

Wartości funkcji $w (x) = x^{4} - 4 x^{3} + 6 x^{2} - 4 x + 1$ obliczone według wzoru. Na marginesie: $w (x) = (x - 1)^{4}$ . Prawdziwe wartości zaznaczone na czerwono.

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

10 \cdot (1.1 - 1) \neq 1

co możesz łatwo sprawdzić:

 
octave:7> 10 * (1.1 -1) - 1
ans =  8.8818e-16

Dlatego

W praktyce numerycznej należy wystrzegać się testów w rodzaju
 
if (x == 1.0) 
{
	....
}

Przedstawiony wcześniej model obliczeniowy jest modelem idealistycznym, tzn. zakłada on, że wszystkie operacje są wykonywane bezbłędnie. Dlatego w tym przypadku będziemy mówić o arytmetyce idealnej. W praktyce jednak, np. wykonując obliczenia na maszynie cyfrowej, operacje arytmetyczne na liczbach rzeczywistych wykonywane są z pewnym błędem. Matematycznym modelem arytmetyki maszyny cyfrowej jest arytmetyka $f l_{ν}$ (albo arytmetyka zmiennoprzecinkowa), którą teraz przypomnimy.

Niech będzie zadana liczba naturalna $b$ (jej znaczenie wyjaśni się w następnym rozdziale). Dowolną liczbę rzeczywistą $x \neq 0$ można jednoznacznie przedstawić w postaci

x = s \cdot 2^{c - b} \cdot m,

gdzie $s \in {- 1, 1}$ jest znakiem, liczba całkowita $(c - b)$ cechą, a liczba rzeczywista $m \in [1, 2)$ mantysą liczby $x$ . 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, ang. floating point). Mantysa ma w ogólności nieskończenie wiele cyfr binarnych $f_{j}$ w swoim rozwinięciu dwójkowym,

m = 1 + f \equiv 1 + \sum_{j = 1}^{\infty} f_{j} 2^{- j} = (1 . f_{1} f_{2} f_{3} \dots)_{2},

gdzie $f_{j} \in {0, 1}$ . Wobec tego najczęściej nie będzie mogła być zapamiętana dokładnie w pamięci komputera, gdyż możemy przechować jedynie ograniczoną liczbę cyfr cechy i mantysy.

Reprezentacja zmiennoprzecinkowa

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 $t$ do zapisania mantysy i także określonej liczby bitów $p$ do zapisania cechy danej liczby niezerowej $x$ :

s c_{1} c_{2} \dots c_{p} f_{1} f_{2} \dots f_{t}

(łącznie $1 + p + t$ bitów). 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 $x$ będziemy nazywać liczbę $r d_{ν} (x)$ taką, że

r d_{ν} (x) = (- 1)^{s} \cdot (1 + f) \cdot 2^{c - b},

gdzie $f$ jest liczbą dwójkową postaci $(0 . f_{1} \dots f_{t})_{2}$ , natomiast $c$ jest liczbą naturalną postaci $(c_{1} \dots c_{p})_{2}$ . Na znak liczby, $s$ , przeznaczony jest jeden bit. Wartości $c$ i $f$ dobiera się tak, żeby $r d_{ν} (x)$ była tak bliska $x$ jak to możliwe. Stałą całkowitą $b$ dobiera się tak, by uzyskać zbalansowany zakres cechy $c - b$ (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 $c - b$ .

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

| \frac{r d_{ν} (x) - x}{x} | \leq \frac{1}{2^{t + 1}} .

Liczbę $ν = \frac{1}{2^{t + 1}}$ nazywa się precyzją arytmetyki. Jak widać, ma ma ona wpływ na to, jak wiele cyfr znaczących liczby jest reprezentowanych dokładnie. 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

r d_{ν} (x) = x (1 + ϵ), gdzie | ϵ | \leq ν .

Przykład

Rozważmy bardzo prościutki 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ż w konsekwencji możliwy zakres wartości $c$ to $0, \dots, 3$ , rozsądne jest więc przyjęcie korekty $b = 1$ , dzięki czemu $- 1 \leq c - b \leq 2$ . Z kolei możliwe wartości mantysy to

(1.00)_{2} = 1, (1.01)_{2} = 1.25, (1.10)_{2} = 1.5, (1.11)_{2} = 1.75 .

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

0.500, 0.625, 0.750, 0.8751.000, 1.250, 1.500, 1.7502.000, 2.500, 3.000, 3.5004.000, 5.000, 6.000, 7.000

Liczby maszynowe: reprezentowane dokładnie w pięciobitowej arytmetyce o precyzji

2^{- 2}

. (Przedstawiliśmy tylko liczby nieujemne)

Standard IEEE 754

Plik:Kahan.jpg

William Kahan
Guru arytmetyki zmiennoprzecinkowej i współtwórca [http://www.cs.berkeley.edu/ wkahan/ieee754status/754story.html standardu IEEE754] Zobacz biografię

Z nielicznymi egzotycznymi wyjątkami (np. Cray C90), współczesne procesory używane w komputerach osobistych lub większych, implementują IEEE 754 Floating Point Standard, który definiuje dwa zasadnicze formaty reprezentacji zmiennoprzecinkowej liczb rzeczywistych:


Typ IEEE 754	Pojedycznej 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 $b$ powyżej)	127	1023
Orientacyjny zakres	$1 0^{- 38} \dots 1 0^{+ 38}$	$1 0^{- 308} \dots 1 0^{+ 308}$
Orientacyjna precyzja	$6 \cdot 1 0^{- 8}$	$1 0^{- 16}$

(maksymalna i minimalna wartość cechy $c$ ma specjalne znaczenie). Dodatkowo, w procesorach x86 mamy typ podwójnej rozszerzonej precyzji (także zdefiniowany w IEEE 754 i odpowiadający dokładnie ówczesnym możliwościom procesora Intel 8087; procesory Intela mają zresztą jedną z najlepszych implementacji IEEE 754). Wszystkie operacje arytmetyczne na procesorach x86 są faktycznie wykonywane w takiej precyzji (korzystając z 64 bitów dla reprezentacji mantysy i 15 bitów dla cechy). Należy pamiętać, że odpowiadający mu typ w C long double zajmuje w pamięci 12 bajtów (a nie 80 bitów).

Uwaga

Producenci niektórych procesorów świadomie rezygnują z pełnej implementacji IEEE 754 dla zwiększenia szybkości działania, niestety czasem kosztem dokładności wyniku. Tak dawno temu było w procesorach Cray; tak też działają niektóre instrukcje wektorowe (z tzw. zestawu 3DNow!) w procesorach AMD, które np. wynik dzielenia wektorowego zwracają z precyzją tylko 14 bitów mantysy. Procesor IBM Cell (stosowany w Sony Playstation 3) również także [http://domino.watson.ibm.com/comm/research.nsf/pages/r.arch.innovation.html nie

w pełni implementuje 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:12> x = NaN
x = 1111111111111000000000000000000000000000000000000000000000000000
octave:13> x = 0
x = 0000000000000000000000000000000000000000000000000000000000000000
octave:14> x = Inf
x = 0111111111110000000000000000000000000000000000000000000000000000
octave:15> x = 0.1
x = 0011111110111001100110011001100110011001100110011001100110011010

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

Przykład: Nawet liczba 0.1 nie jest reprezentowana dokładnie!

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

0.1 = (0.0001100110011001 \dots)_{2} .

Ten banalny fakt jest bardzo często przeoczany przez programistów, a w 1991 roku doprowadził nawet do spektakularnej awarii systemu antyrakietowego Patriot. Okazało się, że --- w tajemniczy sposób --- zazwyczaj bezbłędnie trafiające w cel rakiety Patriot traciły skuteczność, gdy przez wiele godzin pozostawały w stanie gotowości.

Plik:Patriot missile launch.jpg

System rakietowy Patriot

Wyjaśnienie zagadki leżało na styku pomiędzy hardware a software rakiety. 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ę tyknięć 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!

Na marginesie zauważmy, że np. liczba $0.125$ jest reprezentowana dokładnie w arytmetyce zmiennoprzecinkowej (dlaczego?) i nie powodowałaby już tego problemu. Podobnie dawałoby się uniknąć problemu, działając wyłącznie na różnicach liczby tyknięć...

Więcej informacji o najróżniejszych katastrofach spowodowanych błędami w programowaniu można przeczytać na stronach Thomasa Huckle.

Uwaga

Nie wszystkie maszyny liczące wykorzystują reprezentację dwójkową. Kiedyś zdarzały się komputery reprezentujące liczby w postaci

x = s \cdot β^{c} \cdot m,

gdzie $β = 8$ lub 16, a nawet 3 (sic!). Do dzisiaj, w podręcznych kalkulatorach najczęściej spotykaną podstawą reprezentacji liczb jest $β = 10$ .

Są także takie realizacje arytmetyki zmiennoprzecinkowej, które nie realizują w pełni standardu IEEE (np. stare komputery Cray) i np. zamiast zaokrąglenia, stosują obcięcie wyniku.

Standard IEEE 754 jest obecnie uaktualniany, jego nowa wersja powinna ukazać się pod koniec 2006 roku.

Nadmiar i niedomiar

W maszynie cyfrowej cecha $c$ liczby rzeczywistej nie może oczywiście mieć dowolnie dużej wartości bezwzględnej, $| c | \leq c_{\max}$ , dlatego nie wszystkie liczby rzeczywiste są w ogóle reprezentowalne. Powoduje to powstanie zjawiska nadmiaru gdy dla liczby $x c > c_{\max}$ , oraz zjawiska niedomiaru gdy $c < - c_{\min}$ . 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ć reprentowana przez zero, przy czym błąd względny reprezentacji wynosi wtedy $1$ a nie $ν$ .

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.

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: Wyznaczanie normy euklidesowej wektora

Jedną z najczęściej wykonywanych operacji na wektorze $x = (x_{1}, \dots, x_{n})^{T} \in R^{n}$ jest obliczenie jego normy euklidesowej,

| | x | |_{2} = \sqrt{x_{1}^{2} + \dots x_{n}^{2}} .

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ż $| | x | |_{2}$ jest reprezentowana, to $x_{1}^{2}$ już nie (np. w arytmetyce podwójnej precyzji $x_{1} = 1 0^{200}$ i $x_{2} = 1$ ).

Łatwym wyjściem z tej sytuacji jest wstępna normalizacja danych tak, by wszystkie nie były większe od 1: niech $M = \max {| x_{i} | : i = 1, \dots, n}$ i wtedy

| | x | |_{2} = \sqrt{x_{1}^{2} + \dots x_{n}^{2}} = M \cdot \sqrt{{(\frac{x_{1}}{M})}^{2} + \dots + {(\frac{x_{1}}{M})}^{2}} .

i teraz suma pod pierwiastkiem jest zawsze pomiędzy 1 a $N$ . Wadą omówionego rozwiązania jest to, że wymaga ono dwukrotnego przejrzenia całego wektora (raz, by znaleźć $M$ , drugi raz --- by policzyć sumę. Na szczęście można go zmodyfikować tak, by działał w jednym przebiegu. Zupełnie inny algorytm podał Moler.

Liczby denormalizowane

Wymaganie, że mantysa jest postaci $1 + f$ , $f \geq 0$ , powoduje, że wokół zera pojawia się coś w rodzaju próżni. Formalnie, liczby mniejsze niż $2^{1 - 1023}$ 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 $1 0^{- 323}$ .

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

Działania arytmetyczne w arytmetyce $f l_{ν}$

W arytmetyce $f l_{ν}$ implementującej standard IEEE 754, działania arytmetyczne na liczbach rzeczywistych (a raczej na ich reprezentacjach) są wykonywane dokładnie i tylko wynik jest zaokrąglany. Mamy więc

f l_{ν} (x ◻ y) = r d_{ν} (r d_{ν} (x) ◻ r d_{ν} (y)),

gdzie $◻ \in {+, -, \times, \div}$ , Ogólniej, jeśli $𝒲_{1}$ i $𝒲_{2}$ są wyrażeniami o wartościach rzeczywistych, to dla dowolnych wartości zmiennych

f l_{ν} (𝒲_{1} ◻ 𝒲_{2}) = r d_{ν} (f l_{ν} (𝒲_{1}) ◻ f l_{ν} (𝒲_{2})) .

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.

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned fl_\nu\Big(\sqrt{{\cal W}}\Big) &= \left(\sqrt{fl_\nu({\cal W})}\right) (1+\beta_1),\\ fl_\nu(\cos({\cal W})) &= \left(\cos(fl_\nu({\cal W}))\right)(1+\beta_2), \endaligned}

gdzie $| ϵ_{j} | \leq ν$ , oraz $β_{j} \leq K_{j} ν$ i $K_{j}$ są "niewielkimi" stąłymi.

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

1.3 \cdot 2.4

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

<flash>file=binarysystem4.swf</flash><div.thumbcaption>Mnożenie dwóch liczb rzeczywistych (na przykładzie 5-bitowej arytmetyki)

Podobnie, jeśli $△$ jest operatorem porównania, $△ \in {<, \leq, =, \neq}$ , to wartością wyrażenia logicznego $𝒲_{1} △ 𝒲_{2}$ w $f l_{ν}$ jest dokładna wartość wyrażenia Parser nie mógł rozpoznać (nieznana funkcja „\trianglefl”): {\displaystyle \displaystyle fl_\nu({\cal W}_1)\trianglefl_\nu({\cal W}_2)} .

Dosyć dziwnie w porównaniach zachowuje się specjalna liczba NaN (ang. not-a-number), dla której zawsze zachodzi, że NaN $\neq$ NaN. Liczba NaN pojawia się jako wynik zabronionych operacji matematycznych, np. $0 / 0, \sqrt{- 2},$ Inf - Inf, itp., i także propaguje się w dalszych obliczeniach.

Działania arytmetyczne nie są łą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

Wbrew pozorom, fakt, że nie mamy dostępu do arytmetyki nieskończonej precyzji może mieć daleko idące konsekwencje, o czym przekonaliśmy się na początku wykładu.

Praktyczne wyznaczanie precyzji arytmetyki

Aby wyznaczyć precyzję używanej przez nas arytmetyki możemy wykonać prosty test. Pomyślmy, jaka jest najmniejsza dodatnia liczba $ϵ_{mach}$ , która dodana do jedności da w wyniku liczbę większą od 1.0 (liczbę $ϵ_{mach}$ nazywa się czasem epsilonem maszynowym, macheps). Jasne jest, że w przypadku arytmetyki IEEE 754 jest to liczba równa podwojonej precyzji arytmetyki, $2^{- t}$ , gdzie $t$ jest liczbą cyfr mantysy $f$ . 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, 80-bitowej. 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

Ćwiczenie

Sprawdź, jak zmienią się wyniki, gdy wykorzystasz w swoim programie (zarówno w wersji 1, jak w wersji 2) opcje kompilacji:

gcc -O3
gcc -ffast-math
gcc -O3 -ffast-math

Spróbuj objaśnić te wyniki, wspomagając się ewentualnie dokumentacją kompilatora.

Dodaj link: LAPACK daje gotową funkcję, DLAMCH (dla liczb podwójnej precyzji) i SLAMCH (dla pojedynczej precyzji), pozwalającą stwierdzić eksperymentalnie, jakie są parametry używanej arytmetyki, m.in. zakres liczb reprezentowalnych, w jakim systemie reprezentowana jest mantysa, oraz oczywiście precyzję arytmetyki i liczbę cyfr mantysy. Polecamy analizę kodu źródłowego LAPACK/dlamch1.f oraz lekturę prac

Malcolm M. A. (1972) Algorithms to reveal properties of floating-point arithmetic. Comms. of the ACM, 15, 949-951.
Gentleman W. M. and Marovich S. B. (1974) More on algorithms that reveal properties of floating point arithmetic units.

Comms. of the ACM, 17, 276-277.

na których oparto tę funkcję LAPACKa. Poniżej przykład zastosowania tej funkcji i wyniki uzyskane na procesorze x86.

 
/* Wyznaczanie epsilona maszynowego i innych charakterystyk arytmetyki podwójnej
precyzji z wykorzystaniem funkcji DLAMCH z LAPACKa */
#include <stdio.h>
#include <math.h>

double dlamch_(char *CMACH); /* funkcja DLAMCH z LAPACKa */

int main(void)
{
char CMACH;

	CMACH = 'e';
	printf("Epsilon maszynowy: %g\n", dlamch_(&CMACH));
	CMACH = 'b';
	printf("Podstawa arytmetyki: %g\n", dlamch_(&CMACH));
	CMACH = 'n';
	printf("Liczba bitów mantysy: %g\n", dlamch_(&CMACH));
	CMACH = 'u';
	printf("Zakres: %g ", dlamch_(&CMACH));
	CMACH = 'o';
	printf("... %g\n", dlamch_(&CMACH));
	CMACH = 'r';
	return(0);
}

 
Epsilon maszynowy: 2.22045e-16
Podstawa arytmetyki: 2
Liczba bitów mantysy: 53
Zakres: 2.22507e-308 ... 1.79769e+308

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

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

Przykład

Rozważmy banalne zadanie wyznaczenia iloczynu $N$ liczb z tablicy $x$ ,

s = x_{0} \cdot \dots \cdot x_{N - 1} .

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, $\pm$ Inf lub 0).

Naturalnie, zamiast dokładnych wartości $x_{0}, \dots x_{N - 1}$ , będziemy mieli w komputerze jedynie ich reprezentacje, ${\tilde{x}}_{i} = r d_{ν} (x_{i}) = x_{i} (1 + δ_{i})$ , przy czym $| δ_{i} | \leq ν$ .

Oznaczając ${\tilde{s}}_{i}$ wyznaczoną numerycznie wartość iloczynu po $i$ -tym kroku pętli, mamy, że

{\tilde{s}}_{i + 1} = f l_{ν} ({\tilde{s}}_{i} \times {\tilde{x}}_{i}) = {\tilde{s}}_{i} \cdot {\tilde{x}}_{i} \cdot (1 + ϵ_{i}),

gdzie znów $| ϵ_{i} | \leq ν$ . Ostatecznie więc, wyznaczona wartość iloczynu, $\tilde{s}$ spełnia

\tilde{s} = x_{0} \dots x_{N - 1} \cdot Π_{i = 0}^{N - 1} (1 + ϵ_{i}) (1 + δ_{i}) .

Ponieważ $Π_{i = 0}^{N - 1} (1 + ϵ_{i}) = (1 + ℰ)$ , gdzie, z pominięciem małych wyższego rzędu, $| ℰ | \leq N ν$ , dostajemy ostatecznie

\tilde{s} = s \cdot (1 + E),

gdzie $| E | \leq 2 N ν$ . Ponieważ w arytmetyce podwójnej precyzji $ν \approx 1 0^{- 16}$ , 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 $1 0^{- 13}$ !

Powyższe rozumowanie, a także intuicja często wyrażana przez osoby postronne, prowadzi do przypuszczenia, że:

"Duży błąd względny wyniku jest możliwy dopiero po kumulacji błędów zaokrągleń po przeprowadzeniu bardzo wielu działań arytmetycznych."

Jednak to jest to całkowicie fałszywy pogląd, o czym świadczy kolejny, bardzo znamienny przykład.

Przykład: Redukcja cyfr przy odejmowaniu

Tym razem nasze zadanie jest znacznie prostsze od poprzedniego. Trzeba wyznaczyć po prostu różnicę dwóch liczb:

s = a - b .

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

\tilde{s} = f l_{ν} (r d_{ν} (a) - r d_{ν} b) = (a (1 + δ_{a}) - b (1 + δ_{b})) (1 + ϵ),

Stąd po prostych oszacowaniach

| \frac{\tilde{s} - s}{s} | \leq 2 \frac{| a | + | b |}{| a - b |} \cdot ν .

A więc, gdy $a \approx b$ , to $\frac{| a | + | b |}{| a - b |} \approx \infty$ i w efekcie możemy utracić nawet wszystkie znaczące cyfry wyniku! To zjawisko właśnie 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".

Przy okazji zauważmy, że prowadząc identyczną analizę dla sumy dwóch liczb $a + b$ , gdzie $a$ i $b$ są tego samego znaku, dostajemy oszacowanie błędu względnego równe $2 ν$ , niezależnie od wartości liczbowych $a$ i $b$ !

Skutki zjawiska redukcji cyfr przy odejmowaniu mogą być dramatyczne i ujawnić się w nadspodziewanie elementarnych sytuacjach.

Numeryczne kłopoty z wyznaczaniem pierwiastków trójmianu kwadratowego

Niech $a, p, q > 0$ . Korzystając ze szkolenego wzoru na pierwiastki równania kwadratowego $a x^{2} - 2 p x + q = 0$

x_{1, 2} = \frac{1}{a} (p \pm \sqrt{Δ}),

gdzie $Δ = p^{2} - q a > 0$ , możemy natknąć się na trudności, gdy jeden z pierwiastków jest bardzo bliski zera (tzn. gdy $p \approx \sqrt{Δ}$ ).

Taka sytuacja np. powstaje gdy rozważać jednostajnie opóźniony ruch pocisku wystrzelonego z działka przeciwlotniczego do celu lecącego na małej wysokości. Czas trafienia w cel jest --- przy pominięciu oporu powietrza --- rozwiązaniem równania kwadratowego, przy czym czs krótki odpowiada bezpośredniemu trafieniu w cel, a czas długi --- wystrzeleniu pocisku z wyprzedzeniem wysoko w górę i poczekaniu, aż spadając trafi w cel od góry. Oczywiście, żaden artylerzysta nie będzie się interesował tym ostatnim przypadkiem: interesujące jest więc dokładne (bo cel leci szybko) wyznaczenie mniejszego pierwiastka.

Niestety, skoro $p \approx \sqrt{Δ}$ , to wyznaczając mniejszy pierwiastek $x_{1}$ 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, $x_{2} = \frac{1}{a} (p + \sqrt{Δ})$ ! Dokładając do tego wzór Viete'a,

x_{1} x_{2} = \frac{q}{a},

dostajemy inny wzór na $x_{1}$ , nie zawierający feralnego odejmowania. Poniżej demonstrujemy program w C, testujący jakość obu podejść.

# include <stdio.h>
# include <math.h>

/* w(x) = ax^2 - 2px + q = 0  */
/* delta = 4(p^2 - qa) */

double const a = 2.1, q = 1e-6, p=1.1;

double w(double x) /* wartość wielomianu w punkcie x */
{
	return(a*x*x - 2.0*p*x + q);
}

int main(void)
{
	double x1, x2, x1v, X1, X1v, X2;
	double Delta; /* wartość Delty liczymy w podwójnej   precyzji */
	float  delta; /* wartość delty liczymy w pojedynczej precyzji */
		
	delta = Delta = sqrt(p*p - q*a);
	printf("Wielomian w(x) = %e x^2 - %e x + %e.\nDelta = %e\n", a, 2*p, q, delta);

	/* pierwiastki liczone wzorem szkolnym, z niedokładną deltą */
	x1 = (p - delta)/a;
	x2 = (p + delta)/a;
	
	/* mniejszy pierwiatek, liczony z mało dokładną deltą, ale lepszym
	wzorem: Viete'a */
	x1v = (q/a)/x2;
	
	/* pierwiastki liczone wzorem szkolnym, z dokładniejszą Deltą */
	X1 = (p - Delta)/a;
	X2 = (p + Delta)/a;
	
	/* mniejszy pierwiatek, liczony z dokładniejszą Deltą, ale lepszym
	wzorem: Viete'a */
	X1v = (q/a)/X2;

	printf("\nPierwiastki z mało dokładną deltą:\n");
	printf(" Wzór szkolny: x1  = %e x2 = %e\n Wzór Viete'a: x1v = %e x2 = j.w.\n", 
	x1,x2,x1v);
	printf("\nPierwiastki z dokładniejszą Deltą:\n");
	printf(" Wzór szkolny: X1  = %e X2 = %e\n Wzór Viete'a: X1v = %e X2 = j.w.\n", 
	X1,X2,X1v);
	printf("\nWzględna zmiana wartości pierwiastka:\n");
	printf("   (x1 - x1v)/x1v = %e\n", (x1-x1v)/x1v);
	printf("   (x1v -X1v)/X1v = %e\n", (x1v-X1v)/X1v);
	printf("   (x2 -  X2)/X2  = %e\n", (x2-X2)/X2);
	
	printf("\nWartość wielomianu w wyznaczonych punktach:\n");
	printf(" w(x1)  = %e\n w(x1v) = %e w(X1v) = %e\n w(x2)  = %e\n w(X2) = %e\n ", 
		w(x1),w(x1v),w(X1v),w(x2),w(X2));
	
	return(0);
}

 
Wielomian w(x) = 2.100000e+00 x^2 - 2.200000e+00 x + 1.000000e-06.
Delta = 1.099999e+00

Pierwiastki z mało dokładną deltą:
 Wzór szkolny: x1  = 4.427774e-07 x2 = 1.047619e+00
 Wzór Viete'a: x1v = 4.545456e-07 x2 = j.w.

Pierwiastki z dokładniejszą Deltą:
 Wzór szkolny: X1  = 4.545457e-07 X2 = 1.047619e+00
 Wzór Viete'a: X1v = 4.545457e-07 X2 = j.w.

Względna zmiana wartości pierwiastka:
   (x1 - x1v) / x1v = -2.589022e-02
   (x1v -X1v) / X1v = -1.123337e-08
   (x2 -  X2) / X2  = 1.123337e-08

Wartość wielomianu w wyznaczonych punktach:
 w(x1)  = 2.589022e-08
 w(x1v) = 1.123337e-14, w(X1v) = -3.194985e-23
 w(x2)  = 2.589022e-08
 w(X2) = -1.357688e-17

Jak więc widzimy, nawet z niezbyt dokładnie wyznaczoną deltą, mniejszy pierwiastek jesteśmy w stanie wyznaczyć bardzo precyzyjnie --- o ile tylko unikniemy redukcji cyfr przy odejmowaniu.

Przykład: Powiększenie?

Ten obrazek już widzieliśmy na początku wykładu.

Wykres funkcji $w (x) = x^{4} - 4 x^{3} + 6 x^{2} - 4 x + 1 = (x - 1)^{4}$ wyznaczony (w bardzo bliskiej okolicy zera) na dwa sposoby w arytmetyce podwójnej precyzji.

Wyjaśnieniem tej niepokojącej obserwacji jest znowu zjawisko redukcji cyfr przy odejmowaniu: wartości $f$ są bliskie zera, a powstają jako suma dużych liczb z przeciwnymi znakami.

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ą".

To oczywiście też nie jest do końca prawdą. 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, o czym mówimy w następnej części wykładu. Na razie prosty przykład:

Przykład: Co oznacza 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, korzystając z manipulacji symbolicznej --- oczywiście bez trudu system stwierdził, że to wyrażenie upraszcza się do zera.

Następnie zażądaliśmy, by DIGITS --- parametr sterujący "liczbą cyfr znaczących dla liczb zmiennoprzecinkowych", jak to określa manual MUPADa --- przyjął wartość równą 10.

Wymuszając, przez podanie 3.0 zamiast 3 stosowanie w obliczeniach arytmetyki zmiennoprzecinkowej zamiast symbolicznej (pamiętasz, jak to jest w C?) 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. I faktycznie, zazwyczaj taka strategia (choć kosztowna) jest satysfakcjonująca!

MN03: Różnice pomiędzy wersjami

Wersja z 18:11, 1 wrz 2006

Spis treści

Arytmetyka zmiennoprzecinkowa

Reprezentacja zmiennoprzecinkowa

Standard IEEE 754

Nadmiar i niedomiar

Liczby denormalizowane

Działania arytmetyczne w arytmetyce $f l_{ν}$

Praktyczne wyznaczanie precyzji arytmetyki

Numeryczne kłopoty z wyznaczaniem pierwiastków trójmianu kwadratowego

Menu nawigacyjne

Działania na stronie

Opcje strony

Narzędzia osobiste

Nawigacja

Szukaj

Narzędzia

@@ Linia 1: / Linia 1: @@
 =Arytmetyka zmiennoprzecinkowa=
@@ Linia 924: / Linia 923: @@
 kosztowna) jest
 satysfakcjonująca!
-=Własności zadania obliczeniowego i algorytmu numerycznego=
-==Uwarunkowanie zadania obliczeniowego==
-Jak zobaczyliśmy w poprzednich przykładach, dane, jakimi dysponujemy wykonując
-zadanie obliczeniowe są z natury rzeczy wartościami zaburzonymi. Źródłem tych
-zaburzeń są:
-* błąd reprezentacji danych w arytmetyce zmiennoprzecinkowej (0.1 nie jest
-równe dokładnie <math>\displaystyle 1/10</math>)
-* błędy w parametrach będących rezultatem poprzednich obliczeń (chcemy
-rozwiązać równanie <math>\displaystyle f(x) = a</math>, ale <math>\displaystyle a</math> jest rezultatem innej symulacji), a także
-* błędy pomiarowe w zadaniach praktycznych (chcemy policzyć
-numeryczną prognozę pogody, ale temperaturę wyjściową znamy tylko z dokładnością
-do 0.1 stopnia, co gorsza --- z niektórych stacji w ogóle nie mamy danych)
-Okazuje się, że powszechna intuicja, że małe zaburzenia danych powinny dawać
-małe zaburzenia wyniku, nie znajduje potwierdzenia nawet w bardzo prostych
-przypadkach. Z drugiej strony, umiejętność oceny jakościowego <strong>wpływu
-zaburzenia danych na wynik</strong> jest kapitalna w świecie obliczeń numerycznych w
-ogólności, a w szególności --- inżynierskich.
-Wprowadza się pojęcie <strong>uwarunkowania</strong> zadania, to znaczy jego podatności na
-zaburzenia danych. Dla przejrzystości, przypuśćmy, że nasze zadanie obliczeniowe
-polega na wyznaczeniu <math>\displaystyle f(x)</math> dla danego <math>\displaystyle x</math>.
-<div class="thumb tright"><div><flash>file=XXX.swf</flash><div.thumbcaption>Zadanie obliczeniowe i jego odporność na zaburzenia</div></div></div>
-<!--
-[[Image:MNcondition.png|thumb|450px|center|Naszym zadaniem jest wyznaczenie, dla <math>\displaystyle x\in X</math>, wartości
-<math>\displaystyle f(x)\in Y</math>.]]
-[[Image:MNcondition2.png|thumb|450px|center|Jaki będzie rozrzut wyników, gdy <strong>lekko</strong> zaburzymy
-dane?]]
-[[Image:MNcondition3.png|thumb|450px|center|Jeśli równie mały, co zaburzenie, powiemy, że zadanie
-jest dobrze uwarunkowane (jego wynik jest mało podatny na zaburzenia danych).]]
-[[Image:MNcondition4.png|thumb|450px|center|Może jednak zdarzyć się, że zadanie jest źle
-uwarunkowane, i małe zaburzenie danych skutkuje dużym rozrzutem wyników.]]
-[[Image:MNcondition5.png|thumb|450px|center|Wtedy nawet bliskie sobie punkty w X, przekształcenie
-<math>\displaystyle f</math> może odwzorowywać w punkty bardzo od siebie odległe. Jest to sytuacja
-skrajnie niekorzystna w zastosowaniach, a zwłaszcza --- w obliczeniach numerycznych.]]
--->
-Jak bardzo będzie odległe
-<math>\displaystyle f(\widetilde{x})</math>, gdy <math>\displaystyle \widetilde{x}\approx x</math>? Rozważa się dwa przypadki:
-* uwarunkowanie względne: jak względne zaburzenie danych wpływa na błąd
-względny wyniku:
-<center><math>\displaystyle
-\frac{||f(x) - f(\widetilde{x})||}{||f(x)||} \leq  \mbox{cond} _{rel}(f,x) \cdot \frac{||x - \widetilde{x}||}{||x||}
-</math></center>
-Najmniejszy mnożnik <math>\displaystyle  \mbox{cond} _{rel}(f,x)</math> spełniający powyższą nierówność
-nazywamy współczynnikiem uwarunkowania (względnego) zadania obliczenia <math>\displaystyle f(x)</math>
-dla danego <math>\displaystyle x</math>.
-* uwarunkowanie bezwzględne: jak bezwzględne zaburzenie danych wpływa na błąd
-bezwzględny wyniku:
-<center><math>\displaystyle
-||f(x) - f(\widetilde{x})|| \leq  \mbox{cond} _{abs}(f,x) \cdot ||x - \widetilde{x}||
-</math></center>
-Najmniejszy mnożnik <math>\displaystyle  \mbox{cond} _{abs}(f,x)</math>  spełniający powyższą nierówność
-nazywamy współczynnikiem uwarunkowania (bezwzględnego) zadania obliczenia <math>\displaystyle f(x)</math>
-dla danego <math>\displaystyle x</math>.
-Powiemy, że zadanie jest
-* dobrze uwarunkowane w punkcie <math>\displaystyle x</math>, gdy <math>\displaystyle  \mbox{cond} (f,x) \approx 1</math>,
-* źle uwarunkowane w punkcie <math>\displaystyle x</math>, gdy <math>\displaystyle  \mbox{cond} (f,x) \gg 1</math>,
-* źle postawione w punkcie <math>\displaystyle x</math>, gdy <math>\displaystyle  \mbox{cond} (f,x) = +\infty</math>.
-Teraz już rozumiemy, że redukcja cyfr przy odejmowaniu jest tylko
-odzwierciedleniem faktu, że zadanie odejmowania dwóch bliskich liczb jest po
-prostu zadaniem źle uwarunkowanym!
-<div style="margin-top:1em; padding-top,padding-bottom:1em;">
-<span  style="font-variant:small-caps;">Przykład: Uwarunkowanie zadania obliczenia sumy</span>
-<div class="solution">
-Właściwie już wcześniej sprawdziliśmy, że zadanie obliczenia <math>\displaystyle s(x,y) = x + y</math> ma
-<center><math>\displaystyle
- \mbox{cond} _{abs}(s, (a,b)) = 1, \qquad  \mbox{cond} _{rel}(s, (a,b)) = \frac{|a|+|b|}{|a+b|}
-</math></center>
-Tak więc, gdy <math>\displaystyle a\approx -b</math>, to <math>\displaystyle  \mbox{cond} _{rel}(s, (a,b)) \approx +\infty</math> i zadanie
-jest bardzo źle uwarunkowane. Nawet małe zaburzenie względne danych może
-skutkować bardzo dużym zaburzeniem wyniku. Zgodnie z prawem Murphy'ego,
-najczęściej rzeczywiście tak będzie...
-</div></div>
-<div style="margin-top:1em; padding-top,padding-bottom:1em;">
-<span  style="font-variant:small-caps;">Przykład</span>
-<div class="solution" style="margin-left:1em;">
-Dla różniczkowalnej funkcji skalarnej <math>\displaystyle f : R \rightarrow R</math> mamy
-<center><math>\displaystyle
-|f(x) - f(\widetilde{x}| \approx |f'(x) | | x - \widetilde{x} |
-</math></center>
-i w konsekwencji dla zadania obliczenia <math>\displaystyle f(x)</math> dla danego <math>\displaystyle x</math> mamy, przy
-założeniu małych zaburzeń,
-<center><math>\displaystyle
- \mbox{cond} _{abs}( f, x) = |f'(x)|, \qquad  \mbox{cond} _{rel}( f, x) =
-\frac{|f'(x)|\cdot|x|}{|f(x)|}.
-</math></center>
-</div></div>
-Możnaby myśleć, że złe uwarunkowanie zawsze jest szkodliwe w praktyce
-numerycznej. Najczęściej właśnie tak jest istotnie. Jednak w praktyce
-numerycznej sporadycznie zdarza się, że [[sec:invit|Dodaj link: złe uwarunkowanie pewnego podzadania nie tylko
-nie pogarsza sytuacji, ale wręcz pomaga]] szybciej rozwiązać zadanie główne!
-==Rozkład algorytmu względem informacji==
-<strong>Algorytm</strong> to dokładnie określona i dozwolona w danym modelu
-obliczeniowym sekwencja akcji, pozwalająca na rozwiązanie danego
-zadania (w sposób dokładny lub przybliżony).
-Z każdym algorytmem związany jest operator
-<center><math>\displaystyle {\bf ALG}:\,F\longrightarrowG,
-</math></center>
-taki że <math>\displaystyle {\bf ALG}(f)</math> jest wynikiem działania algorytmu
-w arytmetyce idealnej dla danej <math>\displaystyle f</math>.
-Zauważmy, że wynik <math>\displaystyle {\bf ALG}(f)</math> działania algorytmu nie
-zależy bezpośrednio od <math>\displaystyle f</math>, ale raczej od <strong>informacji</strong>
-o <math>\displaystyle f</math> (uzyskanej dzięki poleceniu <math>\displaystyle {\cal IN}</math>). Informacja
-ta może być <strong>pełna</strong> albo tylko <strong>częściowa</strong>.
-Informacja jest pełna gdy, np.
-<math>\displaystyle f=(f_1,\ldots,f_n)\inR^n</math> i wczytamy wszystkie
-współrzędne <math>\displaystyle f_i</math>. Informacja może być częściowa, gdy
-<math>\displaystyle f</math> jest funkcją. Wtedy wiele danych może posiadać tę
-samą informację, co łatwo zaobserwować na przykładzie
-zadania całkowania.
-Niech <math>\displaystyle N:F\to\cup_{n=0}^\inftyR^n</math> będzie
-<strong>operatorem informacji</strong>, tzn.
-<center><math>\displaystyle N(f)\,=\,(y_1,y_2,\ldots,y_n)
-</math></center>
-jest informacją o <math>\displaystyle f</math> zebraną przy idealnej realizacji
-algorytmu. Zauważmy, że nformacja jest pełna gdy <math>\displaystyle N</math> jest
-przekształceniem różnowartościowym, tzn. jeśli
-<math>\displaystyle f_1\nef_2</math> implikuje <math>\displaystyle N(f_1)\neN(f_2)</math>.
-W przeciwnym przypadku mamy do czynienia z informacją
-częściową.
-Każdy algorytm <math>\displaystyle {\bf ALG}</math> może być przedstawiony jako złożenie
-operatora informacji i pewnego operatora
-<math>\displaystyle \varphi:N(F)\toG</math> zdefiniowanego równością
-<center><math>\displaystyle \varphi\left(N(f)\right)\,=\,{\bf ALG}(f).
-</math></center>
-Zauważmy, że w przypadku informacji częściowej zwykle nie
-istnieje algorytm dający dokładne rozwiązanie zadania dla
-każdej danej <math>\displaystyle f\inF</math>, ponieważ dla danych o tej samej
-informacji mogą istnieć różne rozwiązania.
-==Problem wyboru algorytmu==
-Wybór algorytmu jest najistotniejszą częścią całego procesu
-numerycznego rozwiązywania zadania. Kierujemy się przy tym przede
-wszystkim następującymi kryteriami:
-* dokładnością algorytmu,
-* złożonością algorytmu,
-* własnościami numerycznymi algorytmu.
-Przez dokładność algorytmu rozumiemy różnicę między
-rozwiązaniem dokładnym <math>\displaystyle S(f)</math>, a rozwiązaniem
-<math>\displaystyle {\bf ALG}(f)</math> dawanym przez algorytm w arytmetyce idealnej.
-Jeśli <math>\displaystyle {\bf ALG}(f) = S(f)</math>,
-<math>\displaystyle  \forall f \in F</math>,
-to algorytm nazywamy <strong>dokładnym</strong>.
-Mówiąc o złożoności, mamy na myśli złożoność pamięciową
-(zwykle jest to liczba stałych i zmiennych używanych przez
-algorytm), jak również złożoność obliczeniową.
-Na złożoność obliczeniową algorytmu dla danej <math>\displaystyle f</math> składa
-się koszt uzyskania infomacji <math>\displaystyle y=N(f)</math> (zwykle jest on
-proporcjonalny do liczby wywołań polecenia <math>\displaystyle {\cal IN}</math>), oraz
-koszt <strong>kombinatoryczny</strong> przetworzenia tej informacji, aż do
-uzyskania wyniku <math>\displaystyle \varphi(y)</math>. Koszt kombinatoryczny zwykle
-mierzymy liczbą operacji arytmetycznych wykonywanych przez
-algorytm.
-Przez własności numeryczne algorytmu rozumiemy jego
-własności przy realizacji w arytmetyce <math>\displaystyle fl_\nu</math>. Temu
-ważnemu tematowi poświęcimy teraz osobny paragraf.
-==Numeryczna poprawność algorytmu==
-Pożądane jest, aby algorytm dawał "dobry" wynik zarówno
-w arytmetyce idealnej, jak i w arytmetyce <math>\displaystyle fl_\nu</math>. Niestety,
-jak zobaczymy, nie zawsze jest to możliwe. Nawet jeśli algorytm
-jest dokładny to w wyniku jego realizacji w <math>\displaystyle fl_\nu</math> możemy
-otrzymać wynik <math>\displaystyle fl_\nu({\bf ALG}(f))</math> daleko odbiegający od
-<math>\displaystyle S(f)</math>. W szczególności, prawie zawsze mamy
-<center><math>\displaystyle S(f)\,\ne\,fl_\nu\left({\bf ALG}(f)\right).
-</math></center>
-Zauważmy również, że o ile z reguły znamy dokładne zachowanie
-się algorytmu w arytmetyce idealnej dla danej informacji, to nie
-można tego samego powiedzieć o jego zachowaniu się w arytmetyce
-<math>\displaystyle fl_\nu</math>. W związku z tym powstaje pytanie, jak kontrolować błąd
-algorytmów wynikający z błędów zaokrągleń i jakie algorytmy
-uznamy za te o najwyższej jakości numerycznej.
-Istnienie błędów reprezentacji liczb rzeczywistych powoduje,
-że informacja <math>\displaystyle y=N(f)</math> o danej <math>\displaystyle f</math> nie jest w
-ogólności reprezentowana dokładnie. Znaczy to, że zamiast na
-informacji dokładnej, dowolny algorytm będzie operować na
-informacji <strong>nieco zaburzonej</strong> <math>\displaystyle y_\nu</math>, tzn. zaburzonej na
-poziomie błędu reprezentacji. Tak samo wynik dawany przez algorytm
-będzie w ogólności zaburzony na poziomie błędu reprezentacji.
-W najlepszym więc wypadku wynikiem działania algorytmu w <math>\displaystyle fl_\nu</math>
-będzie <math>\displaystyle (\varphi(y_\nu))_\nu</math> zamiast <math>\displaystyle \varphi(y)</math>. Algorytmy
-dające tego rodzaju wyniki uznamy za posiadające najlepsze
-własności numeryczne w arytmetyce <math>\displaystyle fl_\nu</math> i nazwiemy numerycznie
-poprawnymi.
-Dokładniej, powiemy, że ciąg rzeczywisty
-<math>\displaystyle a_\nu=(a_{\nu,1},\ldots,a_{\nu,n})</math>
-(a właściwie rodzina ciągów <math>\displaystyle \{a_\nu\}_\nu</math>) jest
-<strong>nieco zaburzonym</strong> ciągiem <math>\displaystyle a=(a_1,\ldots,a_n)</math>, jeśli
-istnieje stała <math>\displaystyle K</math> taka, że dla wszystkich dostatecznie
-małych <math>\displaystyle \nu</math> zachodzi
-<center><math>\displaystyle
-  |a_{\nu,j} - a_j|\,\le\,K\,\nu\,|a_j|,\qquad 1\le j\le n,
-</math></center>
-albo ogólniej
-<center><math>\displaystyle
-  \|a_\nu - a\|\,\le\,K\,\nu\,\|a\|,
-</math></center>
-gdzie <math>\displaystyle \|\cdot\|</math> jest pewną normą w <math>\displaystyle R^n</math>. W pierwszym
-przypadku mówimy o zaburzeniu "po współrzędnych", a w drugim
-o zaburzeniu w normie <math>\displaystyle \|\cdot\|</math>.
-Zauważmy, że niewielkie zaburzenia po współrzędnych pociągają
-za sobą niewielkie zaburzenia w normie. Rzeczywiście, zachodzi wówczas
-<center><math>\displaystyle \|a_\nu - a\|_\infty \,=\, \max_{1\le j\le n} |a_{\nu,j} - a_j|
-  \,\le\,K\,\nu\,\max_{1\le j\le n} |a_j|\,=\,K\,\nu\,\|a\|_\infty,
-</math></center>
-i korzystając z faktu, że w przestrzeni skończenie wymiarowej
-wszystkie normy są równoważne otrzymujemy dla pewnych stałych
-<math>\displaystyle K_1</math> i <math>\displaystyle K_2</math>
-<center><math>\displaystyle \|a_\nu - a\|\,\le\,K_1\|a_\nu-a\|_\infty\,\le\,
-    K_1 K\,\nu\,\|a\|_\infty\,\le\,K_2 K_1 K\,\nu\,\|a\|,
-</math></center>
-czyli nierówność dla zaburzenia w normie, ze stałą <math>\displaystyle K = K_2 K_1 K</math>.
-{{definicja|Algorytm numerycznie poprawny||
-Algorytm <math>\displaystyle {\bf ALG}</math> rozwiązywania zadania
-nazywamy <strong>numerycznie poprawnym</strong> w zbiorze danych
-<math>\displaystyle F_0\subsetF</math>, jeśli dla każdej danej <math>\displaystyle f\inF_0</math>
-wynik <math>\displaystyle fl_\nu({\bf ALG}(f))</math> działania algorytmu w arytmetyce
-<math>\displaystyle fl_\nu</math> można zinterpretować jako nieco zaburzony wynik
-algorytmu w arytmetyce idealnej dla nieco zaburzonej informacji
-<math>\displaystyle y_\nu=(N(f))_\nu\inN(F)</math> o <math>\displaystyle f</math>, przy czym
-poziom zaburzeń nie zależy od <math>\displaystyle f</math>.
-Formalnie znaczy to, że istnieją stałe <math>\displaystyle K_1</math>, <math>\displaystyle K_2</math>, oraz
-<math>\displaystyle \nu_0>0</math> takie, że spełniony jest następujący warunek.
-Dla dowolnej <math>\displaystyle \nu\le\nu_0</math> oraz informacji <math>\displaystyle y\inN(F_0)</math>
-można dobrać <math>\displaystyle y_\nu\inN(F)</math> oraz
-<math>\displaystyle \left(\varphi(y_\nu)\right)_\nu</math> takie, że
-<center><math>\displaystyle \|y_\nu-y\|\,\le\,K_1\,\nu\,\|y\|,
-</math></center>
-<center><math>\displaystyle \|\left(\varphi(y_\nu)\right)_\nu - \varphi(y_\nu)\|\,\le\,
-     K_2\,\nu\,\|\varphi(y_\nu)\|,
-</math></center>
-oraz
-<center><math>\displaystyle fl_\nu\left({\bf ALG}(f)\right)\,=\,
-      fl_\nu\left(\varphi(N(f))\right)\,=\,
-       \left(\varphi(y_\nu)\right)_\nu.
-</math></center>
-}}
-[[Image:MNcondition7.png|thumb|450px|center|Numerycznie poprawny algorytm daje w arytmetyce <math>\displaystyle fl_\nu</math> wynik <math>\displaystyle ALG(N(x))</math>, który daje
-się zinterpretować jako mało zaburzony wynik <math>\displaystyle f(y)</math> zadania na mało zaburzonych
-danych <math>\displaystyle x</math>.]]
-Zauważmy,że jeśli <math>\displaystyle f\inR^n</math>,
-<math>\displaystyle N(f)=(f_1,\ldots,f_n)</math>, oraz algorytm jest
-dokładny, <math>\displaystyle {\bf ALG}\equiv\varphi\equivS</math>, to numeryczną
-poprawność algorytmu można równoważnie zapisać jako
-<center><math>\displaystyle fl_\nu\left({\bf ALG}(f)\right)\,=\,
-  \left(S(f_\nu)\right)_\nu.
-</math></center>
-==Rola uwarunkowania zadania==
-Niech <math>\displaystyle {\bf ALG}(\cdot)=\varphi(N(\cdot))</math> będzie algorytmem numerycznie
-poprawnym dla danych <math>\displaystyle F_0\subsetF</math>. Wtedy jego błąd w <math>\displaystyle fl_\nu</math>
-można oszacować następująco:
-<center><math>\displaystyle \aligned \lefteqn{ \|S(f)-fl_\nu({\bf ALG}(f))\| \;=\;
-     \|S(f)-\left(\varphi(y_\nu)\right)_\nu\| } \\
-   &\le & \|S(f)-\varphi(y)\|\,+\,
-          \|\varphi(y)-\varphi(y_\nu)\|\,+\,
-          \|\varphi(y_\nu)-\left(\varphi(y_\nu)\right)_\nu\| \\
-   &\le & \|S(f)-{\bf ALG}(f)\|\,+\,
-          \|\varphi(y)-\varphi(y_\nu)\|\,+\,
-          K_2\,\nu\,\|\varphi(y_\nu)\| \\
-   &\le & \|S(f)-{\bf ALG}(f)\|\,+\,
-          (1 + K_2 \nu) \|\varphi(y)-\varphi(y_\nu)\|\,+\,
-          K_2\,\nu\,\|\varphi(y)\|,
-\endaligned</math></center>
-przy czym <math>\displaystyle \|y_\nu-y\|\,\le\,K_1\,\nu\,\|y\|</math>. Stąd
-w szczególności wynika, że jeśli algorytm jest numerycznie
-poprawny i ciągły ze względu na informację <math>\displaystyle y</math>, to
-<center><math>\displaystyle \lim_{\nu\to 0}\,\|S(f)-fl_\nu({\bf ALG}(f))\|\,=\,
-      \|S(f)-{\bf ALG}(f)\|.
-</math></center>
-To znaczy, że dla dostatecznie silnej arytmetyki algorytm będzie
-się zachowywał w <math>\displaystyle fl_\nu</math> prawie tak jak w arytmetyce idealnej.
-Z powyższych wzorów wynika, że błąd w <math>\displaystyle fl_\nu</math> algorytmu
-numerycznie poprawnego zależy w dużym stopniu od:
-* dokładności algorytmu w arytmetyce idealnej,
-* dokładności <math>\displaystyle \nu</math> arytmetyki <math>\displaystyle fl_\nu</math>,
-* wrażliwości algorytmu na małe względne zaburzenia informacji <math>\displaystyle y</math>.
-Ponieważ dwa pierwsze punkty są raczej oczywiste, poświęcimy
-trochę więcej uwagi jedynie trzeciemu.
-Jeśli <math>\displaystyle \varphi</math> spełnia warunek Lipschitza ze stałą <math>\displaystyle L</math>,
-a dokładniej
-<center><math>\displaystyle \|\varphi(y_\nu)-\varphi(y)\|\,\le\,L\,\|y_\nu-y\|,
-</math></center>
-to
-<center><math>\displaystyle \aligned \lefteqn{ \|S(f)-fl_\nu({\bf ALG}(f))\|} \\
-     &\le & \|S(f)-{\bf ALG}(f)\|\,+\,
-       (1+K_2\nu)L\|y_\nu-y\|\,+\,K_2\nu\|\varphi(y)\|  \\
-     &\le & \|S(f)-{\bf ALG}(f)\|\,+\,
-        (1+K_2\nu)LK_1\nu\|y\|\,+\,K_2\nu\|\varphi(y)\|.
-\endaligned</math></center>
-W tym przypadku błędy zaokrągleń zwiększają błąd bezwzględny
-algorytmu proporcjonalnie do <math>\displaystyle \nu</math>.
-Bardziej jednak interesuje nas błąd <strong>względny</strong>. Wybierzmy
-"małe" <math>\displaystyle \eta\ge 0</math> i przypuśćmy, że
-<center><math>\displaystyle \|\varphi(y_\nu)-\varphi(y)\|\;\le\;
-    M\,K_1\,\nu\,\max(\eta,\|\varphi(y)\|),
-</math></center>
-dla pewnej <math>\displaystyle M</math> niezależnej od <math>\displaystyle y</math>, tzn. błąd względny informacji,
-<math>\displaystyle \|y_\nu-y\|\le K_1\nu\|y\|</math>, przenosi się na błąd względny
-wyniku (w arytmetyce idealnej) ze "współczynnikiem wzmocnienia"
-<math>\displaystyle M</math>, albo na błąd bezwzględny ze współczynnikiem <math>\displaystyle M\eta</math>.
-(Zauważmy, że gdybyśmy wzięli <math>\displaystyle \eta=0</math> to dla <math>\displaystyle y</math> takiej, że
-<math>\displaystyle \varphi(y)=0</math> musiałoby być <math>\displaystyle \varphi(y_\nu)=0</math>, co zwykle, choć
-nie zawsze, nie jest prawdą.) Wtedy
-<center><math>\displaystyle \aligned \|S(f)-fl_\nu({\bf ALG}(f))\|
-   & \le & \|S(f)-{\bf ALG}(f)\|+
-     (1 + K_2 \nu) M K_1 \nu \max (\eta, \|\varphi(y)\|)+
-       K_2 \nu \|\varphi(y)\| \\
-    &= \|S(f)-{\bf ALG}(f)\|\,+\,\nu\,
-        \Big(\,MK_1(1+K_2\nu)+K_2\Big)\max(\eta,\|\varphi(y)\|).
-\endaligned</math></center>
-W szczególności, gdy algorytm jest dokładny i korzysta z pełnej
-informacji o <math>\displaystyle f</math>, tzn. <math>\displaystyle S\equiv{\bf ALG}\equiv\varphi</math>, to
-błąd
-<center><math>\displaystyle \frac{\|S(f)-fl_\nu({\bf ALG}(f))\|}
-       {\max (\eta, \|S(f)\|)} \;\le\;
-        \Big( M K_1 (1+K_2\nu) + K_2\Big)\,\nu
-        \,\approx\,(M\,K_1\,+\,K_2)\,\nu.
-</math></center>
-Stąd wynika, że jeśli <math>\displaystyle (MK_1+K_2)\nu\ll 1</math> to błąd względny
-algorytmu w <math>\displaystyle fl_\nu</math> jest mały, o ile <math>\displaystyle \|S(f)\|\ge\eta</math>.
-Błąd względny jest proporcjonalny do dokładności <math>\displaystyle \nu</math>,
-arytmetyki <math>\displaystyle fl_\nu</math>, współczynników proporcjonalności <math>\displaystyle K_i</math>
-algorytmu numerycznie poprawnego, oraz do wrażliwości <math>\displaystyle M</math>
-zadania <math>\displaystyle S</math> na małe względne zaburzenia danych.
-Zwróćmy uwagę na istotny fakt, że interesują nas właściwie
-tylko te zaburzenia danych (informacji), które powstają przy
-analizie algorytmu numerycznie poprawnego. I tak, jeśli algorytm
-jest numerycznie poprawny z pozornymi zaburzeniami danych w normie,
-to trzeba zbadać wrażliwość zadania ze względu na zaburzenia
-danych w normie. Jeśli zaś mamy pozorne zaburzenia
-"po współrzędnych" (co oczywiście implikuje pozorne zaburzenia
-w normie) to wystarczy zbadać uwarunkowanie zadania ze względu na
-zaburzenia "po współrzędnych", itd.
-<div style="margin-top:1em; padding-top,padding-bottom:1em;">
-<span  style="font-variant:small-caps;">Przykład: Iloczyn skalarny</span>
-<div class="solution">
-Załóżmy. że dla danych ciągów rzeczywistych o ustalonej
-długości <math>\displaystyle n</math>, <math>\displaystyle a_j</math>, <math>\displaystyle b_j</math>, <math>\displaystyle 1\le j\le n</math>, chcemy obliczyć
-<center><math>\displaystyle S(a,b)\,=\,\sum_{j=1}^n a_j b_j.
-</math></center>
-Rozpatrzymy algorytm dokładny zdefiniowany powyższym wzorem
-i korzystający z pełnej informacji o kolejnych współrzednych.
-Oznaczmy przez <math>\displaystyle \tilde a_j</math> i <math>\displaystyle \tilde b_j</math> reprezentacje liczb
-<math>\displaystyle a_j</math> i <math>\displaystyle b_j</math> w <math>\displaystyle fl_\nu</math>, <math>\displaystyle \tilde a_j=a_j(1+\alpha_j)</math>,
-<math>\displaystyle \tilde b_j=b_j(1+\beta_j)</math>, oraz przez <math>\displaystyle \gamma_j</math> i <math>\displaystyle \delta_j</math>
-błędy względne powstałe przy kolejnych mnożeniach i dodawaniach.
-Oczywiście <math>\displaystyle |\alpha_j|,|\beta_j|, |\gamma_j|, |\delta_j|\le\nu</math>.
-Otrzymujemy
-<center><math>\displaystyle \aligned fl_\nu\left(\sum_{j=1}^n a_jb_j\right) &=
-    \Big(\,(fl_\nu(\sum_{j=1}^{n-1}a_jb_j)\,+\,\tilde a_n\tilde b_n
-       (1+\gamma_n)\,\Big)(1+\delta_n)\,=\,\ldots \\
-   &= \bigg(\cdots\Big(
-       \tilde a_1\tilde b_1(1+\gamma_1)+\tilde a_2\tilde b_2
-        (1+\gamma_2)\Big)(1+\delta_2) \\
-   & & \qquad\qquad\qquad\qquad +\cdots+
-       \tilde a_n\tilde b_n(1+\gamma_n)\bigg)(1+\delta_n) \\
-   &= \tilde a_1\tilde b_1(1+\gamma_1)(1+\delta_2)
-                 \cdots(1+\delta_n)\\
-   & & \qquad\qquad\qquad\qquad +\cdots+\tilde a_j
-       \tilde b_j(1+\gamma_j)(1+\delta_j)\cdots(1+\delta_n) \\
-   &= \sum_{j=1}^n a_jb_j(1+e_j),
-\endaligned</math></center>
-gdzie w przybliżeniu (tzn. gdy <math>\displaystyle \nu\to 0</math>) mamy <math>\displaystyle |e_1|\leq (n+2)\nu</math>
-i <math>\displaystyle |e_j|\leq (n-j+4)\nu</math>, <math>\displaystyle 2\le j\le n</math>. Algorytm naturalny jest więc
-numerycznie poprawny w całym zbiorze danych, gdyż wynik otrzymany
-w <math>\displaystyle fl_\nu</math> można zinterpretować jako dokładny wynik dla danych
-<math>\displaystyle a_{\nu,j}=a_j</math> i <math>\displaystyle b_{\nu,j}=b_j(1+e_j)</math>, przy czym
-<math>\displaystyle \|b_\nu-b\|_p\leq (n+2)\nu\|b\|_p</math>.
-Zobaczmy teraz, jak błąd we współrzędnych <math>\displaystyle b_j</math> wpływa
-na błąd wyniku. Mamy
-<center><math>\displaystyle \aligned \Big|\sum_{j=1}^n a_jb_j-fl_\nu\Big(\sum_{j=1}^n a_jb_j\Big)\Big|
-    &= \Big|\sum_{j=1}^n a_jb_j-\sum_{j=1}^n a_jb_j(1+e_j)\Big| \\
-    &= \Big|\sum_{j=1}^n e_ja_jb_j\Big|
-      \,\le\, \sum_{j=1}^n |e_j||a_jb_j| \\
-    &\leq  (n+2)\nu\sum_{j=1}^n |a_jb_j|.
-\endaligned</math></center>
-Stąd dla <math>\displaystyle \eta\ge 0</math>
-<center><math>\displaystyle \frac{|\sum_{j=1}^n a_jb_j-fl_\nu(\sum_{j=1}^n a_jb_j)|}
-       {\max(\eta,|\sum_{j=1}^n a_jb_j|)} \,\leq\,
-         K_{\eta}\,(n+2)\,\nu,
-</math></center>
-gdzie
-<center><math>\displaystyle K_\eta\,=\,K_\eta(a,b)\,=\,\frac{\sum_{j=1}^n |a_jb_j|}
-             {\max(\eta,|\sum_{j=1}^n a_jb_j|) }.
-</math></center>
-Zauważmy, że jeśli iloczyny <math>\displaystyle a_jb_j</math> są wszystkie dodatnie
-albo wszystkie ujemne, to <math>\displaystyle K_\eta=1</math>, tzn. zadanie jest dobrze
-uwarunkowane, a błąd względny jest zawsze na poziomie co najwyżej
-<math>\displaystyle n\nu</math>. W tym przypadku algorytm zachowuje się bardzo dobrze, o ile
-liczba <math>\displaystyle n</math> składników nie jest horendalnie duża. W ogólności
-jednak <math>\displaystyle K_\eta</math> może być liczbą dowolnie dużą i wtedy nie możemy
-być pewni uzyskania dobrego wyniku w <math>\displaystyle fl_\nu</math>.
-</div></div>
-<!--
-<div style="margin-top:1em; padding-top,padding-bottom:1em;">
-<span  style="font-variant:small-caps;">Przykład: Pierwiastki trójmianu</span>
-<div class="solution">
-Rozpatrzymy teraz zadanie obliczenia wszystkich pierwiastków
-rzeczywistych równania kwadratowego.
-Będziemy zakładać, że model obliczeniowy dopuszcza obliczanie
-pierwiastków kwadratowych z liczb nieujemnych oraz
-<math>\displaystyle fl_\nu(\sqrt{x})=rd_\nu(\sqrt{rd_\nu(x)})</math>.
-Okazuje się, że nie umiemy pokazać numerycznej poprawności
-"szkolnego" algorytmu obliczającego pierwiastki równania
-bezpośrednio ze wzorów omawianych powyżej. Można jednak pokazać
-numeryczną poprawność drobnej jego modyfikacji wykorzystującej
-wzory Viete'a.
-{{algorytm|||
-<pre>
-Delta = p*p - q;
-if  (Delta == 0)
-      OUT(p);
-else
-	if  (Delta > 0)
-	{
-		Delta1 = sqrt(d);
-		if  (p >= 0)
-		{
-			x1 = p + Delta1;
-			x2 = q/z1;
-		}
-		else
-		{
-			x2 = p - Delta1;
-			x1 = q/ź2;
-		}
-		OUT(x1);  OUT(x2);
-	}
-</pre>}}
-Mamy bowiem
-<center><math>\displaystyle \aligned fl_\nu(\Delta(p,q)) &= \Big(p^2(1+\alpha)^2(1+\epsilon_1)-q(1+\beta)\Big)
-                          (1+\epsilon_2) \\
-    &= \left( p^2-q\frac{(1+\beta)}{(1+\alpha)^2(1+\epsilon_1)}\right)
-           (1+\epsilon_2)(1+\alpha)^2(1+\epsilon_1) \\
-    &= \Big(p^2-q(1+\delta)\,\Big)(1+\gamma) \,=\,
-          \Delta(p,q(1+\delta))(1+\gamma),
-\endaligned</math></center>
-gdzie <math>\displaystyle |\delta|,|\gamma|\leq 4\nu</math>. Wyróżnik obliczony w <math>\displaystyle fl_\nu</math>
-jest więc nieco zaburzonym wyróżnikiem dokładnym dla danych
-<math>\displaystyle p</math> i <math>\displaystyle q_\nu=q(1+\delta)</math>. W szczególności
-<center><math>\displaystyle  \mbox{sgn} (fl_\nu(\Delta(p,q)))= \mbox{sgn} (\Delta(p,q_\nu)).
-</math></center>
-Jeśli <math>\displaystyle p\ge 0</math> to
-<center><math>\displaystyle \aligned fl_\nu(x1(p,q)) &= \Big(p(1+\alpha)+
-         \sqrt{fl_\nu(\Delta(p,q))}(1+\epsilon_3)\Big)(1+\epsilon_4) \\
-   &= \Big(p(1+\alpha)+\sqrt{\Delta(p,q_\nu)
-       (1+\gamma)}(1+\epsilon_3)\Big)(1+\epsilon_4) \\
-   &= \left(p+\sqrt{\Delta(p,q_\nu)}\frac{\sqrt{1+\gamma}(1+\epsilon_3)}
-         {1+\alpha}\right)(1+\epsilon_4)(1+\alpha) \\
-   &= \left(p+\sqrt{\Delta(p,q_\nu)}\right)(1+e_1),
-\endaligned</math></center>
-gdzie <math>\displaystyle |e_1|\leq 6\nu</math>. Zauważmy, że ostatnia równość
-zachodzi dlatego, że dodajemy liczby tego samego znaku. (Inaczej
-<math>\displaystyle |e_1|</math> mogłaby być dowolnie duża i tak byłoby w algorytmie
-szkolnym.) Dla drugiego pierwiastka mamy
-<center><math>\displaystyle fl_\nu(x2(p,q))\,=\,\frac {q(1+\beta)}{fl_\nu(x1(p,q))}(1+\epsilon_5)
-  \,=\,\frac{q_\nu}{fl_\nu(x1(p,q))}(1+e_2),
-</math></center>
-gdzie <math>\displaystyle |e_2|\le 8\nu</math>.
-Podobny wynik otrzymalibyśmy dla <math>\displaystyle p<0</math>. Algorytm zmodyfikowany
-jest więc numerycznie poprawny, gdyż otrzymane w <math>\displaystyle fl_\nu</math> pierwiastki
-są nieco zaburzonymi dokładnymi pierwiatkami dla danych
-<math>\displaystyle p_\nu=p</math> i <math>\displaystyle q_\nu=q(1+\delta)</math>.
-Aby oszacować błąd algorytmu, wystarczy zbadać uwarunkowanie
-zadania ze względu na zaburzenie danej <math>\displaystyle q</math>, ponieważ pokazaliśmy,
-że zaburzenia <math>\displaystyle p</math> można przenieść na zaburzenia <math>\displaystyle q</math> i wyniku.
-Niestety, choć algorytm jest numerycznie poprawny, zaburzenia
-<math>\displaystyle q</math> mogą sprawić, że nawet znak wyróżnika <math>\displaystyle \Delta</math> może być
-obliczony nieprawidłowo. Na przykład dla <math>\displaystyle p=1</math> i <math>\displaystyle q=1\pm 10^{t+1}</math>
-mamy <math>\displaystyle \Delta(p,q)=\mp 10^{t+1}</math>, ale
-<math>\displaystyle \Delta(rd_\nu(p),rd_\nu(q))=\Delta(1,1)=0</math>. Ogólnie
-<center><math>\displaystyle |fl_\nu(\Delta(p,q))-\Delta(p,q)|\,\leq\,4\nu(p^2+2|q|),
-</math></center>
-a więc tylko dla <math>\displaystyle |\Delta(p,q)|=|p^2-q|>4\nu (p^2+2|q|)</math>
-możemy być pewni obliczenia właściwego znaku <math>\displaystyle \Delta</math>. Przy
-tym warunku oraz <math>\displaystyle \Delta>0</math> błąd danych przenosi się w
-normie euklidesowej na błąd wyniku następująco:
-<center><math>\displaystyle \aligned \lefteqn{ \Big( (x1(p,q) - x1(p,q_\nu))^2
-                 +(x2(p,q) - x2(p,q_\nu))^2 \Big)^{1/2} } \\
-  &= \frac{\sqrt 2 |\delta q|} {\sqrt{p^2-q}+\sqrt{p^2-q_\nu}}
-   \,\leq\, 2\sqrt 2 \nu\frac{|q|}{\sqrt{p^2-q}} \\
-  &= 2\sqrt 2 \nu\frac{|q|/p^2}{\sqrt{1-q/p^2}
-        \max(\eta/|p|,\sqrt{2(1+(1-q/p^2))}) } \\
-  & & \qquad\qquad\qquad\cdot\max(\eta,(x1(p,q)^2+x2(p,q)^2)^{1/2}).
-\endaligned</math></center>
-Stąd widać, że zadanie jest dobrze uwarunkowane dla <math>\displaystyle q/p^2\ll 1</math>
-i może być źle uwarunkowane dla <math>\displaystyle q/p^2\approx 1</math>. W ostatnim
-przypadku nie możemy być pewni otrzymania dobrego wyniku w <math>\displaystyle fl_\nu</math>.
-</div></div>
--->

MN03: Różnice pomiędzy wersjami

Wersja z 18:11, 1 wrz 2006

Arytmetyka zmiennoprzecinkowa

Reprezentacja zmiennoprzecinkowa

Standard IEEE 754

Nadmiar i niedomiar

Liczby denormalizowane

Działania arytmetyczne w arytmetyce flν

Praktyczne wyznaczanie precyzji arytmetyki

Numeryczne kłopoty z wyznaczaniem pierwiastków trójmianu kwadratowego

Menu nawigacyjne

Szukaj

Działania arytmetyczne w arytmetyce $f l_{ν}$