Laboratorium 5: 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 ¹

#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/80 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

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  „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

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)]

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] && 
   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 do tego poniższe polecenie:

<<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]

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

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)

albo

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 (), pozwalającą na obliczenia w dowolnej precyzji.

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

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

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')

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')

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.

1. Dodawanie: fadd()

2. Odejmowanie: fsub()

3. Negacja (zmiana znaku): fneg()

4. Iloczyn: fmul()

5. 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 s= 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

1. Zaproponować eksperyment obliczeniowy wyznaczający liczbę bitów i dokładność prowadzonych obliczeń.

2. 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).

3. 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ń.