MN03LAB

Ćwiczenia. Własności arytmetyki zmiennopozycyjnej.

Ćwiczenie: Równe i równiejsze

Wyjaśnij, dlaczego w arytmetyce podwójnej precyzji IEEE 754 mamy

 
octave:19> 2006/1e309
ans = 0
octave:20> 2.006/1e306
ans =  2.0060e-306
octave:21> (2006/1000)/(1e309/1000)
ans = 0

Oczywiście, teoretycznie wszystkie trzy liczby powinny być sobie równe (i niezerowe).

Wskazówka

Kluczem jest stwierdzenie, jaka jest reprezentacja liczby $10^{309}$ w podwójnej precyzji? A $10^{306}$?

Ćwiczenie: Szeregi zbieżne(?)

Podaj przykłady zbieżnych szeregów postaci $\sum _{n=1}^{\mathrm{\infty }}{a}_n$, których $N$-te sumy częściowe obliczone w arytmetyce pojedynczej precyzji algorytmem

 
suma = 0.0;
for n = 1..N
	suma += <math>\displaystyle a_n</math>;

będą

• ograniczone niezależnie od $N$, albo
• numerycznie rozbieżne, to znaczy takie, że dla bardzo dużych $N$ zachodzi

suma == Inf.

Wykonaj to samo zadanie, ale podając przykłady szeregów rozbieżnych (w arytmetyce dokładnej).

Rozwiązanie

Rozpatrz na przykład takie szeregi:

• Szereg zerowy (numerycznie dostajesz oczywiście też 0)
• $1+10^{2006}-10^{2006}+0+\dots =1$ oraz $1+10^{2006}-1+0+\dots =10^{2006}$

(numerycznie dostaniesz, odpowiednio, NaN (dlaczego?) i Inf.

• Szereg harmoniczny (numerycznie sumy częściowe w pewnym momencie przestaną

rosnąć).

• Szereg jedynkowy.

Ćwiczenie

Dla kolejnych $N$ wyznacz $N$-tą sumę częściową szeregu Taylora dla funkcji wykładniczej, gdy $x=-90$:

$e^x\approx \sum _{n=0}^N\frac{x^n}{n!},$

korzystając z algorytmu podanego w poprzednim zadaniu. Porównaj błąd względny i bezwzględny w porównaniu do wartości wyznaczonej z wykorzystaniem funkcji bibliotecznej exp(). Powtórz następnie dla $x=10$.

Czy --- zgodnie z teorią matematyczną --- sumy te dążą do wartości $e^x$. Objaśnij dokładnie, co się stało.

Rozwiązanie

Zastosowanie naszego algorytmu dla $x=-90$ daje w wyniku (dla arytmetyki podwójnej precyzji) sumę równą około $10^{30}$, gdy oczywiście naprawdę $\exp (-90)\approx 0$). Problem jest skutkiem zjawiska redukcji cyfr przy odejmowaniu: aby nasza suma była równa zeru, w pewnym momencie musimy odjąć duży składnik od innego podobnego, by w rezultacie dostać coś małego. Spróbuj wychwycić ten moment i wartość dodawanego składnika.

Dla $x=10$ mamy, dla rosnącego $N$, całkiem niezły błąd względny ostatecznego wyniku.

Ćwiczenie

Już wcześniej stwierdziliśmy, że wyznaczanie $e\approx (1+1/n{)}^n$ dla dużego $n$ nie jest dobrym pomysłem. Przeprowadź eksperyment numeryczny potwierdzający to stwierdzenie i objaśnij jego wyniki.

Ćwiczenie

Jak wiadomo, szereg harmoniczny $\sum _{n=1}^{\mathrm{\infty }}1/n$ jest rozbieżny. Spróbuj przewidzieć, jaki będzie efekt numerycznego wyznaczenia tej sumy w arytmetyce podwójnej precyzji przy użyciu poniższego kodu.

 
	int dlicznik;
	double dsuma, dstarasuma;
	double dskladnik;
	
	dstarasuma = 0.0; dsuma = 1.0; dlicznik = 1;
	while(dstarasuma != dsuma) 
	{
		dskladnik = 1.0/dlicznik;
		dstarasuma = dsuma;
		dsuma += dskladnik;
		dlicznik++;
	}
	printf("Suma = %e. Liczba składników = %d, składnik = %e\n", 
		dsuma, dlicznik-1, dskladnik);
	 

Wskazówka

Policz, ile mniej więcej obrotów pętli potrzeba, by suma przestała rosnąć.

Rozwiązanie

Nasza suma przekroczy wartość 20, więc aby 20 + dskladnik == 20, musi być dskładnik $\approx 10^{-15}$ (lub więcej). Tymczasem zakres liczb typu integer wynosi niewiele ponad $10^9$, więc po jakimś czasie licznik przyjmie wartości ujemne.

Aby doliczyć się "do końca", musielibyśmy licznik potraktować jako liczbę typu long long double, której zakres sięga $10^{18}$. Ale i tak mielibyśmy kłopot z doczekaniem do końca działalności programu.

Zakładając optymistycznie, że na sekundę zliczamy $10^7$ składników, potrzebowalibyśmy razem około $10^7$ sekund, czyli ponad trzy miesiące...

Ćwiczenie

Jak szybko i na jakiej wysokości musiałby lecieć samolot npla, aby pocisk wystrzeliwany z działka z prędkością 7500 km/h nie trafił w cel, gdy potrzebne pierwiastki liczone są wzorem szkolnym?

Ćwiczenie: Zadanie o wyznaczaniu pierwiastka kwadratowego metodą Newtona

Dla zadanej liczby $a>1$ należy wyznaczyć przybliżenie $\sqrt{a}$, stosując metodę Herona. Zaproponuj dobre przybliżenie początkowe $x_0$ wiedząc, że $a$ jest liczbą maszynową typu double. Ile iteracji wystarczy, by osiągnąć $ϵ$-zadowalające przybliżenie?

Rozwiązanie

Jak pamiętamy, $a=(1+f)2^p=\tilde{f}{2}^{p+1}$, gdzie $\frac{1}{2}\le \tilde{f}<1$ oraz $p\in N$. Stąd

$\sqrt{a}=\sqrt{\tilde{f}}{2}^{(p+1)/2}=\sqrt{g}{2}^k,$

dla pewnego (znanego) $k\in N$, przy czym $\frac{1}{2}\le \sqrt{g}\le 1$. Wobec tego obliczenie $\sqrt{a}$ wymaga po prostu obliczenia $\sqrt{g}$. Najprostsze przybliżenie $\sqrt{g}$ to $\sqrt{g}\approx 1$. Błąd takiego przybliżenia to $|\sqrt{g}-1|\le 0.5$.

Jak można jeszcze bardziej poprawić $x_0$?

Ćwiczenie: Zadanie o wyznaczaniu odwrotności bez dzielenia metodą Newtona

Należy wyznaczyć przybliżenie $\frac{1}{a}$ stosując metodę Newtona do równania $\frac{1}{x}-a=0$. Zaproponuj dobre przybliżenie początkowe $x_0$ wiedząc, że $a$ jest liczbą maszynową typu double. Ile iteracji wystarczy, by osiągnąć $ϵ$-zadowalające przybliżenie?

Rozwiązanie

Jak pamiętamy, $a=(1+f)2^p$, skąd $\frac{1}{a}=\frac{1}{1+f}{2}^{-p}$. Wystarczy więc przybliżyć $\frac{1}{1+f}$. Ponieważ dla $0\le f<1$,

$\frac{1}{2}\le \frac{1}{1+f}\le 1,$

to możemy jako pierwsze przybliżenie wybrać $x_0=\frac{3}{4}{2}^{-p}$. Jak wybrać lepsze $x_0$?

Wtedy mamy następujące oszacowania:

$|x_{0}a-1|=|\frac{3}{4}{2}^{-p}\cdot (1+f)2^p-1|=|\frac{3}{4}(1+f)-1|\le \frac{1}{2}.$

Ponieważ iteracja Newtona spełnia

$|x_{k+1}a-1|=|x_{k}a-1{|}^2=|x_{0}a-1{|}^{2^k},$

to, jeśli chcemy, by błąd względny przybliżenia $x_{k+1}$ spełniał $\frac{|x_{k+1}-\frac{1}{a}|}{\frac{1}{a}}=|x_{k+1}a-1|\le ϵ=2^{-53}$, to musimy wykonać co najwyżej $6$ iteracji. Ponieważ koszt jednej iteracji to 3 flopy, wyznaczenie odwrotności $a$ byłoby 18 razy droższe od mnożenia i dodawania, które mają podobny koszt.

Tymczasem we współczesnych realizacjach dzielenie jest tylko około 4 razy droższe od mnożenia! Możesz sprawdzić, jak to się dzieje np. w procesorach AMD.

Więcej na temat implementacji algorytmów dzielenia oraz wyciągania pierwiastka kwadratowego możesz dowiedzieć się np. z artykułu Petera Marksteina Software Division and Square Root Using Goldschmidt s Algorithms oraz raportów technicznych Stanford Computer Architecture and Arithmetic Group

Ćwiczenie

Niech $0<a_1<a_2<\cdots <a_n$. Czy z punktu widzenia błędów w $f{l}_{\nu }$ lepiej jest policzyć sumę tych liczb w kolejności od najmniejszej do największej, czy odwrotnie?

Rozwiązanie

Od najmniejszej do największej.