Laboratorium 4: Arytmetyka zmiennoprzecinkowa komputerów

Wstęp

Te zajęcia nawiązują do sposobu zapisu liczb opisanego w „części teoretycznej”. Ich celem jest odświeżenie tej wiedzy, zaprojektowanie prostego eksperymentu pozwalającego sprawdzić jak obliczenia są wykonywane (zwłaszcza w arkuszu kalkulacyjnym).

Druga część ma pokazać w jaki sposób można te ograniczenia omijać korzystając ze specjalnych bibliotek.

Proste obliczenia

Jest taki program [1]¹

#include <stdio.h>
int main()
{
    float s;
    double d;
    long double e;
    int i;
    s = d = e = 0.5;
    for(i=1;i<=100;i++)
    {
        s = 3.8F * s * (1.F - s);
        d = 3.8 * d * (1. - d);
        e = 3.8L * e * (1.L - e);
        if (i%10==0) printf("%10i %16.5f %16.5lf %16.5Lf\n", i, s, d, e);
    }
//    printf("%ld\n", sizeof(long double));
    return 0;
}

#include <stdio.h>

int main()

{

float s;

double d;

long double e;

int i;

s = d = e = 0.5;

for(i=1;i<=100;i++)

{

s = 3.8F * s * (1.F - s);

d = 3.8 * d * (1. - d);

e = 3.8L * e * (1.L - e);

if (i%10==0) printf("%10i %16.5f %16.5lf %16.5Lf\n", i, s, d, e);

}

// printf("%ld\n", sizeof(long double));

return 0;

}

nie prowadzi on żadnych skomplikowanych obliczeń przeprowadzając jedynie proste działanie: \[x_{i+1}=\alpha x_i(1-x_i)\] gdzie $\alpha=3.8$, a $x_0=0.5$.

Wykonuje on obliczenia na liczbach typu float (32 bity), double (64 bity) i long double (128 bitów). Program wykonuje 100 iteracji drukując co dziesiątą z nich. Wyniki programu (można go łatwo zapisać na dysku jako caos.c i skompilować poleceniem gcc -o caos caos.c.

Wyniki wyglądają jakoś tak:

        10          0.18510          0.18510          0.18510
        20          0.23951          0.23963          0.23963
        30          0.88445          0.90200          0.90200
        40          0.23023          0.82493          0.82493
        50          0.76124          0.53714          0.53714
        60          0.72004          0.66878          0.66879
        70          0.89952          0.53189          0.53203
        80          0.61599          0.93573          0.93275
        90          0.84367          0.68312          0.79884
       100          0.94858          0.65620          0.23138

10 0.18510 0.18510 0.18510

20 0.23951 0.23963 0.23963

30 0.88445 0.90200 0.90200

40 0.23023 0.82493 0.82493

50 0.76124 0.53714 0.53714

60 0.72004 0.66878 0.66879

70 0.89952 0.53189 0.53203

80 0.61599 0.93573 0.93275

90 0.84367 0.68312 0.79884

100 0.94858 0.65620 0.23138

Jak widać rezultaty różnią się w sposób znaczący. Jedynym wytłumaczeniem jest różna liczba bitów uwzględnianych w obliczeniach. Zatem przypuszczać można że wersja 128-bitowa jest najbliższa rzeczywistości.

Jeżeli wierzyć rozważaniom z [1] „poprawne wyniki” (uzyskane w arytmetyce o 1000 cyfr) są następujące:

           10            0.18509
           20            0.23963
           30            0.90200
           40            0.82492
           50            0.53713
           60            0.66878
           70            0.53202
           80            0.93275
           90            0.79885
          100            0.23161

10 0.18509

20 0.23963

30 0.90200

40 0.82492

50 0.53713

60 0.66878

70 0.53202

80 0.93275

90 0.79885

100 0.23161

i jak widać odrobinę różnią się od najlepszego uzyskanego wyniku.

Okazuje się, że dla $0\le \alpha < 3$ zachowanie obliczeń jest deterministyczne, a dla $3 \le \alpha <4$ — chaotyczne. Takie systemy (chaotyczne) są bardzo czułe na dokładność obliczeń.

Mathematica

Mathematica dla każdego obliczenia potrafi podać precyzję wyniku. Służy do tego funkcja Precision[]. Zazwyczaj obliczenia wykonywane są z „maszynową precyzją“ (co oznacza obliczenia na liczbach podwójnej precyzji (64 bity).

Dodatkowo można zadeklarować w jakiej precyzji mają być wykonywane obliczenia.

Block[{$MinPrecision = 5,$MaxPrecision = 5},x = a~* x * (1 - x)]

1	Block[{$MinPrecision = 5,$MaxPrecision = 5},x = a~* x * (1 - x)]

Znalazłem też instrukcję (której do końca nie rozumiem) powodującą, że wszystkie obliczenia w notatniku będą odbywać się z zadaną precyzją. Wygląda ona jakoś tak:

$PreRead=(#/.s_String/; StringMatchQ[s, NumberString] &amp;&amp; 
   Precision @ ToExpression @ s==MachinePrecision:&gt; s&lt;&gt;"`1000."&amp;);

1 2	$PreRead=(#/.s_String/; StringMatchQ[s, NumberString] && Precision @ ToExpression @ s==MachinePrecision:> s<>"`1000."&);

i powinna być umieszczona na początku notatnika². Precyzja określona jest stałą tekstową na samym końcu polecenia ("`1000."). Demonstruje to prosty przykład (notatnik Mathematici).

Powinno to wystarczyć do zaprogramowania opisywanego przykładu.

Mathematica to pełnowymiarowy język programowania, wyposażony w instrukcje warunkowe (If czy służące do tworzenia pętli (For, Do, While). Prosty przykład programiku z tymi instrukcjami zawiera notatnik.

Jak ktoś ciągle ma kłopoty może skorzystać z kolejnego przykładowego programu.

Pakiet Computer Arithmetics

Mathematica wyposażona jest w pakiet Computer Aruithmetics służący do symulowania obliczeń korzystając z jakiejś „wymyślonej” arytmetyki komputerowe. Pozwala to badać jaki wpływ na na precyzję obliczeń może mieć liczba dostępnych bitów czy zakres zmian wartości wykładnika.

W szczególności pakiet pozwala na symulowanie komputerów o arytmetyce dziesiętnej (nie binarnej). Wydaje się że do prowadzenia złożonych obliczeń inżynierskich taki rodzaj arytmetyki może być lepszy niż stosowana dziś arytmetyka binarna.

Aby z pakietu skorzystać, trzeba go załadować. Służy d tego polecenie iNeeds.

Needs["ComputerArithmetic`"].

1	Needs["ComputerArithmetic`"].

Następnie definiujemy rodzaj arytmetyki poleceniem SetArithmetic[d,b] gdzie d to liczba używanych cyfr (niestety z zakresu 1 do 10), a b to podstawa systemu liczenia z zakresu od 2 do 16. Nie można zatem poszaleć z liczbą używanych w obliczeniach cyfr (co utrudnia nieco nasze zadanie).

Kolejne polecenie to ComputerNumber definiujące obiekt w wybranej arytmetyce:

x = ComputerNumber[0.5]
a = ComputerNumber[3.8]

1 2	x = ComputerNumber[0.5] a = ComputerNumber[3.8]

(Standardowo nie można wykonywać obliczeń na liczbach mieszanych więc trzeba uważać na zapis:

x = x * a * (1 - x)

1	x = x * a * (1 - x)

gdyż jedynka jest z innego świata (obiektem innego typu). Czyli powyższe trzeba zapisać jako:

x = x * a * (ComputerNumber[1] - x)

1	x = x * a * (ComputerNumber[1] - x)

albo

one = ComputerNumber[1]
x = x * a * (one - x)

1 2	one = ComputerNumber[1] x = x * a * (one - x)

Takie obliczenia również można wykonać aby sprawdzić zachowanie naszego problemu. Od razu podpowiadam, że można wykorzystać arytmetykę szesnastkową o 10 cyfrach (uzyskamy większą precyzję niż w przypadku arytmetyki 10 o 10 cyfrach). Ale i tak to za mało.

Python

Język programowania python może być wyposażony w bibliotekę mpmath ([2]), pozwalającą na obliczenia w dowolnej precyzji.

Korzystanie z niej jest stosunkowo proste. Trzeba ją najpierw „załadować”

from mpmath import *

1	from mpmath import *

Polecenia mp.prec i mp.dps definiują precyzję obliczeń. Pierwsze określa ją w bitach, drugie w cyfrach dziesiętnych. Wystarczy korzystać tylko z jednego, gdyż ich wartości zależne są od siebie: $\text{prec}\sim 3.33 \text{ dps}$

>>> mp.dps=100
>>> mp.prec
336

>>> mp.dps=100

>>> mp.prec

336

Największy kłopot jest z deklaracją stałych, które nie mogą być standardowymi literałami.

>>> x=mpf(4)
>>> x
mpf('4.0')
>>> x/mpf('3.')
mpf('1.333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333338')
>>> mp.dps=30
>>> x/mpf('3.')
mpf('1.33333333333333333333333333333327')

>>> x=mpf(4)

>>> x

mpf('4.0')

>>> x/mpf('3.')

mpf('1.333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333338')

>>> mp.dps=30

>>> x/mpf('3.')

mpf('1.33333333333333333333333333333327')

Wprowadzając dane trzeba pamiętać o błędach konwersji z układu dziesiętnego do dwójkowego. W związku z tym spodziewać można się dziwnych zjawisk

>>> mpf(.1)
mpf('0.100000000000000005551115123125783')
>>> mpf('.1')
mpf('0.0999999999999999999999999999999951')
>>> mpf(1)/mpf(10)
mpf('0.0999999999999999999999999999999951')

>>> mpf(.1)

mpf('0.100000000000000005551115123125783')

>>> mpf('.1')

mpf('0.0999999999999999999999999999999951')

>>> mpf(1)/mpf(10)

mpf('0.0999999999999999999999999999999951')

w pierwszym przypadku najpierw robiona jest konwersja z liczby dziesiętnej do binarnej, a później wynik traktowany jest jako liczba wielokrotnej precyzji. W drugim wypadku do funkcji mpf przekazujemy tekst, który jest interpretowany i konwertowany od razu na liczbę wielokrotnej precyzji.

Reszta jest już stosunkowo prosta.

Dodawanie: fadd()
Odejmowanie: fsub()
Negacja (zmiana znaku): fneg()
Iloczyn: fmul()
Iloraz: fdiv()

Chociaż jak już zmienne zostały „zadeklarowane” jako dużej precyzji można korzystać ze standardowych operatorów +-/*.

>>> mp.dps = 100
>>> x = mpf('1.')
>>> y = mpf('3.')
>>> z~= x / y
>>> nprint(z,20)
0.33333333333333333333

>>> mp.dps = 100

>>> x = mpf('1.')

>>> y = mpf('3.')

>>> z~= x / y

>>> nprint(z,20)

0.33333333333333333333

Idle

Można ułatwić sobie pracę w pythonie uruchamiając proste środowisko graficzne zwane idle. Pozwala ono zapisać program korzystając z prostego edytora i bardzo łatwo uruchamiać go.

Jupyter

Znacznie bardziej wygodne może być skorzystanie z notatnika Jupyter. Przygotowałem odpowiedni przykład, który można
pobrać (wraz z innymi) stąd.

Zadania do wykonania

Zaproponować eksperyment obliczeniowy wyznaczający liczbę bitów i dokładność prowadzonych obliczeń.
Przetestować go na matlabie i Mathematici, w arkuszu kalkulacyjnym LibreOffice calc, Blockly i — jeżeli ktoś konto Google posiada dla arkusza kalkulacyjnego Google (i ewentualnie dla Microsoft Office Online).
Pobawić się w pythonie³ i Mathematici prowadzeniem obliczeń w dużej precyzji programując podany na początku przykład i porównując (jakieś wykresy?) wyniki uzyskane dla różnych precyzji obliczeń.

Instrukcja w postaci jednego pliku…

…jest również dostępna.

Problem pojawił się po raz pierwszy podczas prób symulacji izolowanej populacji insektów.↩
Ustalona precyzja będzie obowiązywała we wszystkich równocześnie otwartych notebookach↩
Programowanie w pythonie jest stosunkowo proste ([3], [4], [5], [6]). Można również zerknąć do instrukcji laboratoryjnych do zajęć z Technologii Informacyjnych. Na stronach TI dla informatykji znajduje się również obszerny spis literatury. Natomiast, jeżeli ktoś boi się zaczynać, powinien zerknąc na przygotowaną stronę w Blockly zawierający program i pozwalającą na podejrzenie wersji skonwertowanej do pythona. Skopiować z tego okienka jest trudno, ale się da. Klikamy w nie i w ciemno Ctrl-A, Ctrl-C i zawartość wklejamy do edytora. Później trzeba tylko jeszcze podzielić na linie (i zadbać o właściwe wcięcia, ale jak robić to ostrożne – udaje się).↩

A. Ward, „Some issues on floating-point precision under linux,” Linux gazette, iss. 53, 2000.
[BibTeX] [Download PDF]

@Article{Ward2000,
author = {Alan Ward},
title = {Some issues on floating-point precision under Linux},
journal = {Linux Gazette},
year = {2000},
number = {53},
url = {https://linuxgazette.net/53/ward.html},
month = {May},
owner = {myszka},
timestamp = {2012.10.15},
}

F. Johansson, „mpmath’s documentation,” , 2015.
[BibTeX] [Download PDF]

@Electronic{Johansson2015,
author = {Fredrik Johansson},
title = {{mpmath}’s documentation},
year = {2015},
url = {http://mpmath.org/doc/current/index.html},
owner = {myszka},
timestamp = {2015.11.05},
}

M. Pilgrim, Zanurkuj w pythonie, WikiBooks, 2010.
[BibTeX] [Abstract] [Download PDF]

Niniejszy podręcznik powstaje na podstawie książki Dive into Python Marka Pilgrima (w większości jest to tłumaczenie) udostępnionej na licencji {GNU} Free Documentation License. Kody wszystkich przykładów można pobrać stąd.

@Book{mark-pilgrim-zanurkuj-2010,
author = {Mark Pilgrim},
title = {Zanurkuj w Pythonie},
year = {2010},
publisher = {{WikiBooks}},
url = {http://pl.wikibooks.org/wiki/Python},
abstract = {Niniejszy podręcznik powstaje na podstawie książki Dive into Python Marka Pilgrima (w większości jest to tłumaczenie) udostępnionej na licencji {GNU} Free Documentation License. Kody wszystkich przykładów można pobrać stąd.},
owner = {myszka},
timestamp = {2010.10.27},
}

Python documentation index, 2010.
[BibTeX] [Download PDF]

@Misc{python-2010,
title = {Python documentation index},
year = {2010},
howpublished = {\url{http://www.python.org/doc/}},
month = sep,
url = {http://www.python.org/doc/},
key = {Python},
owner = {myszka},
timestamp = {2015.11.06},
}

Z. A. Shaw, Learn python the hard way, , 2010.
[BibTeX] [Download PDF]

@Book{Shaw2010,
author = {Zed A. Shaw},
title = {Learn Python The Hard Way},
year = {2010},
url = {http://learnpythonthehardway.org/book/},
owner = {myszka},
timestamp = {2010.10.04},
}

Python programming, WikiBooks, 2010.
[BibTeX] [Abstract] [Download PDF]

Python is a general-purpose interpreted programming language. It currently has distributions available for Microsoft Windows, Apple Mac {OS} X, {GNU/Linux,} {BSD,} and many other platforms. There are currently three major implementations: the standard implementation written in C, Jython written in Java, and {IronPython} written in C\# for the {MS} {.NET} environment.

@Book{wiki-python-2010,
title = {Python Programming},
year = {2010},
publisher = {{WikiBooks}},
url = {http://en.wikibooks.org/wiki/Python_Programming},
abstract = {Python is a general-purpose interpreted programming language. It currently has distributions available for Microsoft Windows, Apple Mac {OS} X, {GNU/Linux,} {BSD,} and many other platforms. There are currently three major implementations: the standard implementation written in C, Jython written in Java, and {IronPython} written in C\# for the {MS} {.NET} environment.},
key = {Python},
owner = {myszka},
timestamp = {2010.10.27},
}