Tekst źródłowy strony MN03LAB

<!-- 
Konwertowane  z pliku LaTeX przez latex2mediawiki, zob. http://www.ii.uj.edu.pl/&nbsp;pawlik1/latex2mediawiki.php.
Niezb�dne rozszerzenia i modyfikacje oryginalnego latex2mediawiki 
wprowadzi� przykry@mimuw.edu.pl
-->
 
=Własności arytmetyki zmiennoprzecinkowej=

{{powrot |Metody numeryczne | do strony głównej 
przedmiotu <strong>Metody numeryczne</strong>}}

<div class="mw-collapsible mw-made=collapsible mw-collapsed">
Oglądaj wskazówki i rozwiązania __SHOWALL__<br>
Ukryj wskazówki i rozwiązania __HIDEALL__
</div>

<div style="margin-top:1em; padding-top,padding-bottom:1em;">
<span  style="display: block; background-color:#fefeee; border-bottom: 1px solid #E5E5E5; line-height: 1.1em; padding-bottom: 0.2em; font-variant:small-caps; color:#1A6ABF;">Ćwiczenie: Równe i równiejsze</span> 
<div class="exercise">

Wyjaśnij, dlaczego w arytmetyce podwójnej precyzji IEEE 754 mamy
<div style="font-family: monospace; white-space: pre; border-style: dashed; border-width: thin; border-color: black; margin: 1em; padding:1em; color: #444; background-color:#fdfdfd;"><nowiki>octave:19> 2006/1e309
ans = 0
octave:20> 2.006/1e306
ans =  2.0060e-306
octave:21> (2006/1000)/(1e309/1000)
ans = 0
</nowiki></div>
 
Oczywiście, "teoretycznie" wszystkie trzy liczby powinny być sobie
równe (i niezerowe). 

<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Wskazówka </span><div class="mw-collapsible-content" style="display:none">
<div style="font-size:smaller; background-color:#f9fff9; padding: 1em"> Kluczem jest wyjaśnienie, jaka jest reprezentacja liczby <math>10^{309}</math> w podwójnej precyzji... A reprezentacja <math>10^{306}</math>? </div>
</div></div>

</div></div>

<div style="margin-top:1em; padding-top,padding-bottom:1em;">
<span  style="display: block; background-color:#fefeee; border-bottom: 1px solid #E5E5E5; line-height: 1.1em; padding-bottom: 0.2em; font-variant:small-caps; color:#1A6ABF;">Ćwiczenie: Szeregi zbieżne(?)</span> 
<div class="exercise">

Podaj przykłady <strong>zbieżnych</strong> szeregów postaci <math>\sum_{n=1}^{\infty} a_n</math>, których <math>N</math>-te sumy częściowe
obliczone w arytmetyce pojedynczej precyzji algorytmem
 <div style="margin: 1em; padding:1em; color: #000; background-color:#fcfcfc;"><pre>suma = 0.0;
for n = 1..N
	suma += <math>a_n</math>;
</pre></div>
 
będą 
* ograniczone niezależnie od <math>N</math>, albo
* numerycznie rozbieżne, to znaczy takie, że dla bardzo dużych <math>N</math> zachodzi <code style="color: #006">suma == Inf</code>.
 
Wykonaj to samo zadanie, ale podając przykłady szeregów  <strong>rozbieżnych</strong> (w
arytmetyce dokładnej).
</div></div>

<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwiązanie </span><div class="mw-collapsible-content" style="display:none"><div style="margin-left:1em">  

Rozpatrz na przykład takie szeregi:
* Szereg zerowy (numerycznie dostajesz oczywiście też 0)
* <math>1 + 10^{2006} - 10^{2006} + 0 + \ldots = 1</math> oraz <math>1 + 10^{2006} - 1 + 0 + \ldots = 10^{2006}</math> (numerycznie dostaniesz, odpowiednio, <code style="color: #006">NaN</code> (dlaczego?) i <code style="color: #006">Inf</code>.
* Szereg harmoniczny (numerycznie sumy częściowe w pewnym momencie przestaną rosnąć).
* Szereg jedynkowy.
 
</div></div></div>

<div style="margin-top:1em; padding-top,padding-bottom:1em;">
<span  style="display: block; background-color:#fefeee; border-bottom: 1px solid #E5E5E5; line-height: 1.1em; padding-bottom: 0.2em; font-variant:small-caps; color:#1A6ABF;">Ćwiczenie</span> 
<div class="exercise">

Dla kolejnych <math>N</math>, wyznacz <math>N</math>-tą sumę częściową szeregu Taylora dla funkcji wykładniczej, gdy <math>x=-90</math>:
<center><math>
e^x \approx \sum_{n=0}^N \frac{x^n}{n!}</math>,</center>

korzystając z algorytmu podanego w poprzednim zadaniu. Oblicz błąd względny
i bezwzględny wyznaczonego przybliżenia, w porównaniu 
do wartości wyznaczonej z wykorzystaniem funkcji bibliotecznej <code>exp()</code>. Powtórz następnie dla <math>x=10</math>.

Czy --- zgodnie z teorią matematyczną --- sumy te dążą do wartości
<math>e^{x}</math>. Objaśnij dokładnie, co się stało.
</div></div>

<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwiązanie </span><div class="mw-collapsible-content" style="display:none"><div style="margin-left:1em">  
Zastosowanie naszego algorytmu dla <math>x=-90</math> daje w wyniku (dla arytmetyki
podwójnej precyzji) sumę równą około <math>10^{30}</math>, 
gdy oczywiście naprawdę <math>\exp(-90)
\approx 0</math>). 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 <math>x=10</math> mamy, dla rosnącego <math>N</math>, całkiem niezły błąd względny ostatecznego
wyniku.

</div></div></div>

<div style="margin-top:1em; padding-top,padding-bottom:1em;">
<span  style="display: block; background-color:#fefeee; border-bottom: 1px solid #E5E5E5; line-height: 1.1em; padding-bottom: 0.2em; font-variant:small-caps; color:#1A6ABF;">Ćwiczenie</span> 
<div class="exercise">

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

<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Wskazówka </span><div class="mw-collapsible-content" style="display:none">
<div style="font-size:smaller; background-color:#f9fff9; padding: 1em"> Kiedy da znać o sobie epsilon maszynowy? </div>
</div></div>

</div></div>

<div style="margin-top:1em; padding-top,padding-bottom:1em;">
<span  style="display: block; background-color:#fefeee; border-bottom: 1px solid #E5E5E5; line-height: 1.1em; padding-bottom: 0.2em; font-variant:small-caps; color:#1A6ABF;">Ćwiczenie</span> 
<div class="exercise">

Jak wiadomo, szereg harmoniczny <math>\sum_{n=1}^\infty 1/n</math> 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.

 <div style="margin: 1em; padding:1em; color: #000; background-color:#fcfcfc;"><pre>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 = &#37;e. Liczba składników = &#37;d, składnik = &#37;e\n", 
		dsuma, dlicznik-1, dskladnik);
	 
</pre></div>
 
<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Wskazówka </span><div class="mw-collapsible-content" style="display:none">
<div style="font-size:smaller; background-color:#f9fff9; padding: 1em"> Policz, ile mniej więcej obrotów pętli potrzeba, by suma przestała rosnąć. </div>
</div></div>

</div></div>

<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwiązanie </span><div class="mw-collapsible-content" style="display:none"><div style="margin-left:1em">  
Nasza suma przekroczy wartość 20, więc aby <code>20 + dskladnik == 20</code>, musi
być dskładnik <math>\approx 10^{-15}</math> (lub więcej). Tymczasem zakres liczb typu integer wynosi niewiele ponad <math>10^9</math>, więc po jakimś czasie <code>dlicznik</code> przyjmie wartości ujemne.

</div></div></div>

<!-- 

<div style="margin-top:1em; padding-top,padding-bottom:1em;">
<span  style="display: block; background-color:#fefeee; border-bottom: 1px solid #E5E5E5; line-height: 1.1em; padding-bottom: 0.2em; font-variant:small-caps; color:#1A6ABF;">Ćwiczenie: Zadanie o wyznaczaniu pierwiastka kwadratowego metodą Newtona</span> 
<div class="exercise">

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

<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwiązanie </span><div class="mw-collapsible-content" style="display:none"><div style="margin-left:1em">  

Jak pamiętamy, <math>a = (1+f)2^p = \widetilde{f} 2^{p+1}</math>,
gdzie <math>\frac{1}{2}\leq \widetilde{f} < 1</math> oraz <math>p\in N</math>. Stąd

<center><math>
\sqrt{a} = \sqrt{\widetilde{f}}2^{(p+1)/2} = \sqrt{g} 2^k</math>,</center>

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

Jak można jeszcze bardziej poprawić <math>x_0</math>?

</div></div></div>

-->
 
<div style="margin-top:1em; padding-top,padding-bottom:1em;">
<span  style="display: block; background-color:#fefeee; border-bottom: 1px solid #E5E5E5; line-height: 1.1em; padding-bottom: 0.2em; font-variant:small-caps; color:#1A6ABF;">Ćwiczenie: Zadanie o wyznaczaniu odwrotności bez dzielenia metodą Newtona</span> 
<div class="exercise">

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

<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwiązanie </span><div class="mw-collapsible-content" style="display:none"><div style="margin-left:1em">  
Jak pamiętamy, <math>a = (1+f)2^p</math>, skąd <math>\frac{1}{a} = \frac{1}{1+f}2^{-p}</math>.
Wystarczy więc przybliżyć <math>\frac{1}{1+f}</math>. Ponieważ dla <math>0 \leq f < 1</math>,

<center><math>\frac{1}{2} \leq \frac{1}{1+f} \leq 1</math>,</center>

to możemy jako pierwsze przybliżenie wybrać <math>x_0 = \frac{3}{4}2^{-p}</math> (jak
wybrać jeszcze lepsze <math>x_0</math>?)

Wtedy mamy następujące oszacowania:

<center><math>|x_0a - 1| = |\frac{3}{4}2^{-p}\cdot (1+f)2^p - 1| = |\frac{3}{4}(1+f) - 1|
\leq \frac{1}{2}</math></center>

Ponieważ iteracja Newtona spełnia

<center><math>|x_{k+1}a - 1| = |x_ka-1|^2 = |x_0a-1|^{2^k}</math>,</center>

to, jeśli chcemy, by błąd względny przybliżenia <math>x_{k+1}</math> spełniał <math>\frac{|x_{k+1} -
\frac{1}{a}|}{\frac{1}{a}} = |x_{k+1}a - 1| \leq \epsilon = 2^{-53}</math>, 
musimy wykonać co najwyżej <math>6</math> iteracji. Ponieważ koszt jednej iteracji
to 3 działania arytmetyczne,  wyznaczenie odwrotności <math>a</math> 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ć, [http://www.cs.utexas.edu/users/moore/publications/divide_paper.pdf  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]
</div></div></div>

<div style="margin-top:1em; padding-top,padding-bottom:1em;">
<span  style="display: block; background-color:#fefeee; border-bottom: 1px solid #E5E5E5; line-height: 1.1em; padding-bottom: 0.2em; font-variant:small-caps; color:#1A6ABF;">Ćwiczenie</span> 
<div class="exercise">

Niech <math>0 < a_1 < a_2 < \cdots < a_n</math>. Czy z punktu widzenia 
błędów w <math>fl_\nu</math> lepiej jest policzyć sumę tych liczb w kolejności 
od najmniejszej do największej, czy odwrotnie? 
</div></div>

<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwiązanie </span><div class="mw-collapsible-content" style="display:none"><div style="margin-left:1em">  
Od najmniejszej do największej.
</div></div></div>

<div style="margin-top:1em; padding-top,padding-bottom:1em;">
<span  style="display: block; background-color:#fefeee; border-bottom: 1px solid #E5E5E5; line-height: 1.1em; padding-bottom: 0.2em; font-variant:small-caps; color:#1A6ABF;">Ćwiczenie</span> 
<div class="exercise">

Dlaczego w algorytmie bisekcji rozwiązywania równania <math>f(x)=0</math>, sprawdzając warunek różnych znaków <math>f</math> w krańcach przedziału <math>[a,b]</math>, używamy kryterium 
 <div style="margin: 1em; padding:1em; color: #000; background-color:#fcfcfc;"><pre>if ( sign(f(x)) != sign(f(xl)) )
...	
</pre></div>
 
a nie, matematycznie równoważnego, wyrażenia
 <div style="margin: 1em; padding:1em; color: #000; background-color:#fcfcfc;"><pre>if ( f(x)*f(lewy) < 0 )	
...
</pre></div>
 
<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Wskazówka </span><div class="mw-collapsible-content" style="display:none">
<div style="font-size:smaller; background-color:#f9fff9; padding: 1em"> Nie daj się zwieść, nie chodzi tu o zaoszczędzenie jednego mnożenia... </div>
</div></div>

</div></div>

<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwiązanie </span><div class="mw-collapsible-content" style="display:none"><div style="margin-left:1em">  
Gdy zbliżamy się do miejsca zerowego <math>f</math>, wartości <code>f(x)</code> i <code>f(lewy)</code> mogą być na tyle małe, że ich iloczyn --- poprzez niedomiar --- będzie miał numeryczną wartość zero i w efekcie program może przeskoczyć do niewłaściwej gałęzi pętli <code>if</code>, a tym samym --- wybrać niewłaściwą lokalizację przedziału zawierającego poszukiwany pierwiastek.
</div></div></div>