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órego $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 = 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 dlicznik przyjmie wartości ujemne.

Aby doliczyć się "do końca", musielibyśmy licznik potraktować jako liczbę typu long long double, którego 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, zatem 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 [http://www.informatik.uni-trier.de/Reports/TR-08-2004/rnc6_12_markstein.pdf Software Division and Square Root Using Goldschmidt s Algorithms] oraz [http://arith.stanford.edu/techrep.html 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.

Ćwiczenie

Aby obliczyć $S(a,b)=a^2-b^2$ można zastosować dwa algorytmy: ${\mathbf{𝐀}\mathbf{𝐋}\mathbf{𝐆}}_1(a,b)=a*a-b*b$ oraz ${\mathbf{𝐀}\mathbf{𝐋}\mathbf{𝐆}}_2(a,b)=(a+b)*(a-b)$. Pokazać, że oba algorytmy są numerycznie poprawne, ale drugi z nich wywołuje mniejszy błąd względny wyniku w przypadku, gdy $r{d}_{\nu }(a)=a$ i $r{d}_{\nu }(b)=b$.

Ćwiczenie

Pokazać, że naturalny algorytm obliczania cosinusa kąta między dwoma wektorami Parser nie mógł rozpoznać (nieznana funkcja „\inR”): {\displaystyle \displaystyle a, b\inR^n} ,

$\cos (a,b)=\frac{\sum _{j=1}^n{a}_j{b}_j}{\sqrt{\left(\sum _{j=1}^n{a}_j^2\right)\left(\sum _{j=1}^n{b}_j^2\right)}},$

jest numerycznie poprawny. Oszacować błąd względny wyniku w $f{l}_{\nu }$.

Ćwiczenie

Pokazać, że naturalny algorytm obliczania $\Vert Ax{\Vert }_2$ dla danej macierzy Parser nie mógł rozpoznać (nieznana funkcja „\inR”): {\displaystyle \displaystyle A\inR^{n\times n}} i wektora Parser nie mógł rozpoznać (nieznana funkcja „\inR”): {\displaystyle \displaystyle x\inR^n} jest numerycznie poprawny. Dokładniej,

$f{l}_{\nu }(\Vert Ax{\Vert }_2)=(A+E)x,$

gdzie $\Vert E{\Vert }_2\le 2(n+2)\sqrt{n}\nu \Vert A{\Vert }_2$. Ponadto, jeśli $A$ jest nieosobliwa to

$|f{l}_{\nu }(\Vert Ax{\Vert }_2)-\Vert Ax{\Vert }_2|\le 2(n+2)\sqrt{n}\nu (\Vert A{\Vert }_2\Vert A^{-1}{\Vert }_2)\Vert Ax{\Vert }_2.$

Ćwiczenie

Niech $\mathbf{𝐀}\mathbf{𝐋}\mathbf{𝐆}$ będzie algorytmem numerycznie poprawnym w zbiorze danych $f\in F_0$, przy czym dla małych $\nu$, $f{l}_{\nu }(\mathbf{𝐀}\mathbf{𝐋}\mathbf{𝐆}(f))=\phi (y_{\nu })$, gdzie $\Vert y_{\nu }-y\Vert \le K\nu \Vert y\Vert$ i $K$ nie zależy od $\nu$ i $f$ ($y=N(f)$). Pokazać, że w ogólności $\mathbf{𝐀}\mathbf{𝐋}\mathbf{𝐆}$ nie musi być "numerycznie poprawny po współrzędnych", tzn. w ogólności nie istnieje bezwzględna stała $K_1$ taka, że dla małych $\nu$ i dla dowolnej $f\in F_0$

$|y_{\nu ,j}-y_j|\le K_1\nu |y_j|,1\le j\le n,$

gdzie $y=(y_1,\dots ,y_n)$.

Ćwiczenie

Podaj przykład funkcji $f$, której miejsce zerowe $x^{*}$ ma wspólczynnik uwarunkowania

• mały
• duży

Rozwiązanie

Ponieważ nasze zadanie to wyznaczenie $x^{*}=f^{-1}(0)$, to

${\text{cond}}_{abs}(f^{-1},0)=\frac{1}{f^{\prime }(x^{*})}.$

Znaczy to, że im bardziej płaska jest $f$ w otoczeniu pierwiastka $x^{*}$, tym bardziej nawet małe zaburzenia $f$ mogą spowodować duże przemieszczenie jej miejsca zerowego.

Zauważ, że dla wielokrotnych miejsc zerowych, ${\text{cond}}_{abs}(f^{-1},0)=\mathrm{\infty }$, co zgadza się z intuicją, bo może być, że nawet minimalne zaburzenie $f$ spowoduje, iż miejsc zerowych po prostu nie będzie...