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;
}
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 notatnika2. 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.
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 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
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 pythonie3 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ń.