Reprezentacja zmiennoprzecinkowa liczb

W poniższej prezentacji zostanie przedstawiona w działaniu reprezentacja zmiennoprzecinkowa liczb w arytmetyce komputerowej. Ogólnie daną liczbę będziemy zapisywać w postaci x = ϵ M (b^w), gdzie M jest mantysą liczby, b - bazą systemu liczbowego, w - wykładnikiem , a ϵ znakiem liczby

Clear[Evaluate[Context[]<>"*"]];
<<ComputerArithmetic`

Arithmetic[]

4 oznacza liczbę cyfr mantyssy
RoundToEven: zaokrąglanie do najbliższej liczby maszynowej, w przypadku "remisu" do parzystej mantyssy.
RoundToInfinity: zaokrąglanie do najbliższej liczby maszynowej, w przypadku "remisu" do liczby dalej od zera.
Truncation: po prostu obciąć nadmierne liczby mantyssy.

Poniżej wypisujemy liczbę 2 w obecnej arytmetyce komputerowej:

ComputerNumber[2]

Poniżej wypisujemy liczbę pi w obecnej arytmetyce komputerowej:

ComputerNumber[Pi]

Można też podać znak, mantysę (musi mieć tyle cyfr ile zadane w arytmetyce, tj. aktualnie 4) i wykładnik:

ComputerNumber[1, 3142, -3]

Pomiżej jest wypisana lista wszystkich liczb maszynowych o długości mantysy m=2, bazie b=2, minimalnym wykładniku minexp=-4 i maksymalnym wykładniku maxexp=1 (każda liczba ze względów technicznych jest podana w postaci pary {liczba, 0.1})

repliczby[m_,b_,minexp_,maxexp_]:=Module[{list={}},
For[i=minexp,i≤maxexp,
For[j=b^(m-1),j≤b^m-1,
list=Append[list,{j*(b^i),0.1}];list=Append[list,{-j*(b^i),0.1}];j++];i++];
list
];repliczby[2,2,-4,1]

Na poniższym wykresie możemy zobaczyć jak rozkładają się na osi x wszystkie możliwe do zapisania liczby w systemie zmiennoprzecinkowym (dodatnie) przy ustalonej mantysie, bazie i wykładniku (m=długość mantysy, b=baza, minexp=minimalny wykładnik, maxexp=maksymalny wykładnik r= zakres wykresu od -r do r).

Wyraźnie widać, że wokół zera jest "dziura", a liczby te rozkładają się nierównomiernie, co jest powodem wielu błędów obliczeniowych w tej arytmetyce! Liczby maszynowe rozmieszczone są dyskretnie a odległość między dwoma sąsiednimi liczbami to Ulp (zmienny)

Manipulate[ListPlot[(repliczby[m,b,minexp,maxexp]),PlotStyle->PointSize[0.01], PlotRange->{{-r,r},{-1,1}},PlotLabel->Style["wszystkie liczby maszynowe o danych parametrach",Bold,16]],{{m,2,"długość mantysy"},2,5,1},{{b,2, "baza"},2,10,1},{{minexp,-2,"najmniejszy wykładnik"},-6,-2,1},{{maxexp,2,"największy wykładnik"},2,6,1},{{r,1,"zakres wykresu"},0.1,10,0.1}]

Funkcja SetArithmetic[b] pozwala ustawić bazę dla danego systemu, przykładowo zobaczmy jak wygląda liczba π w bazie dwójkowej:

SetArithmetic[10,2];
ComputerNumber[π]

Poniższy przykład demonstruje, że dodawanie liczb w arytmetyce komputerowej jest operacją niełączną. Liczba (1+(1/4)+(1/9)+...+(1/2000)) ma inną wartość w zależności od tego jaki 
przyjmiemy porządek sumowania! (obliczenia wykonujemy z dokładnością do 4 miejsc po przecinku w systemie dziesiętnym)

Suma w kolejności naturalnej:

SetArithmetic[4,10]

sum=0;Do[sum+=(1/((ComputerNumber[i])^2)),{i,2000}];sum

Sumowanie od końca:

sum=0;Do[sum+=(1/((ComputerNumber[i])^2)),{i,2000,1,-1}];sum

Różnica wyników wynosi:

"1.644000000000000000"-"1.625000000000000000"

Błąd rzędu 1 % dla 2000 wyrazów !

0.019`16.718928450563464/"1.644000000000000000"

Możemy przebadać jak różnica wyników zależy od ilości sumowanych wyrazów:

RoznicaWynikow[n_,precision_,base_]:=Module[{r=0},SetArithmetic[precision,base];sum=0;(Do[sum+=(1/((ComputerNumber[i])^2)),{i,n,1,-1}]);
r=sum;sum=0;
(Do[sum+=(1/((ComputerNumber[i])^2)),{i,n}]);
r=r-sum;
r]

RoznicaWynikow[2000,4,10]

RoznicaWynikow[1001,5,8]

SetArithmetic[4,10]

Poniżej przedstawiono porównanie wartości symbolu nieoznaczonego (1/(2-((Sqrt[2]))^2)) ,który w arytmetykach o różnej bazie wcale nie przyjmuje wartości bliskiej nieskończoności !

SetArithmetic[5,8];
1/(2-ComputerNumber[Sqrt[2]]*ComputerNumber[Sqrt[2]])//N

ComputerArithmetic[(1/(2.`+(Subscript["-2.0001", 8])))]

ComputerNumber[1/(2-ComputerNumber[Sqrt[2]]*ComputerNumber[Sqrt[2]])]

ComputerAritmetic[(1/(2+(Subscript["-2.00010000000000000000", 8])))]

(1/(2+(Subscript["-2.00000010000000000", 11])))//N

(1/(2+(Subscript["-1.55555550000000000000000", 6])))//N

SetArithmetic[4,10];

Teraz przebadamy arytmetyke na tej konkretnej maszynie. Poniższy wykres przedstawia wartość względną błędu obliczenia wartości funkcji f(x) = (x^2) w małym sąsiedztwie argumentu x = 100. Widać, że wartość błędu zmienia się nieregularnie. Jednostki na osi pionowej to Ulp (tj. odległość dwóch sąsiednich liczb maszynowych w otoczeniu wartości).

MicroscopicErrorPlot[x^2,{x,100}]

Poniższy wykres pokazuje z kolei wygląd funkcji w arytmetyce komputerowej, ukazując "schodkową" naturę tych funkcji. Liczbę x podnosimy do potęgi a/b i pokazujemy wykres w okolicy liczby x=c.

Manipulate[MicroscopePlot[x^(a/b),{x,c}],{{a,3},1,10,1},{{b,5},1,10,1},{{c,7},0,10}]

Poniższy wykres pokazuje wartośc Ulp na przedziale (a,b).

Manipulate[Plot[Ulp[x],{x,a,b}],{{a,0.001},0.00001,1},{{b,5},1,10,1}]