commit 8a5c33459d6368b2699d38993e3076066362ed32 Author: jacekpoz Date: Sun Oct 27 21:56:24 2024 +0100 l1 diff --git a/.envrc b/.envrc new file mode 100644 index 0000000..3550a30 --- /dev/null +++ b/.envrc @@ -0,0 +1 @@ +use flake diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..be5f547 --- /dev/null +++ b/.gitignore @@ -0,0 +1,4 @@ +**/*.aux +**/*.log +**/*.pdf +.direnv/ diff --git a/flake.lock b/flake.lock new file mode 100644 index 0000000..aa911b0 --- /dev/null +++ b/flake.lock @@ -0,0 +1,43 @@ +{ + "nodes": { + "nixpkgs": { + "locked": { + "lastModified": 1729850857, + "narHash": "sha256-WvLXzNNnnw+qpFOmgaM3JUlNEH+T4s22b5i2oyyCpXE=", + "owner": "NixOS", + "repo": "nixpkgs", + "rev": "41dea55321e5a999b17033296ac05fe8a8b5a257", + "type": "github" + }, + "original": { + "owner": "NixOS", + "ref": "nixpkgs-unstable", + "repo": "nixpkgs", + "type": "github" + } + }, + "root": { + "inputs": { + "nixpkgs": "nixpkgs", + "systems": "systems" + } + }, + "systems": { + "locked": { + "lastModified": 1681028828, + "narHash": "sha256-Vy1rq5AaRuLzOxct8nz4T6wlgyUR7zLU309k9mBC768=", + "owner": "nix-systems", + "repo": "default", + "rev": "da67096a3b9bf56a91d16901293e51ba5b49a27e", + "type": "github" + }, + "original": { + "owner": "nix-systems", + "repo": "default", + "type": "github" + } + } + }, + "root": "root", + "version": 7 +} diff --git a/flake.nix b/flake.nix new file mode 100644 index 0000000..b2a4954 --- /dev/null +++ b/flake.nix @@ -0,0 +1,30 @@ +{ + description = "on"; + + inputs = { + nixpkgs.url = "github:NixOS/nixpkgs/nixpkgs-unstable"; + systems.url = "github:nix-systems/default"; + }; + + outputs = { nixpkgs, systems, ... }: let + name = "on"; + + forEachSystem = nixpkgs.lib.genAttrs (import systems); + pkgsForEach = nixpkgs.legacyPackages; + in { + devShells = forEachSystem ( + system: let + pkgs = pkgsForEach.${system}; + in { + default = pkgs.mkShell { + inherit name; + + packages = with pkgs; [ + julia + texliveFull + ]; + }; + } + ); + }; +} diff --git a/l1/1.jl b/l1/1.jl new file mode 100755 index 0000000..69b15d8 --- /dev/null +++ b/l1/1.jl @@ -0,0 +1,87 @@ +#!/usr/bin/env julia + +# Jacek Poziemski 272389 + +""" + macheps(T) + +Calculate the machine epsilon iteratively +for the given IEEE 754 floating point type. + +# Arguments +- `T`: a floating point number type <: AbstractFloat +""" +function macheps(T) + one_T = one(T) + ret = one_T + while one_T + (ret / T(2)) > one_T + ret /= T(2) + end + + return ret +end + +""" + eta(T) + +Calculate the eta number iteratively +for the given IEEE 754 floating point type. + +# Arguments +- `T`: a floating point number type <: AbstractFloat +""" +function eta(T) + ret = one(T) + + while ret / T(2) > zero(T) + ret /= T(2) + end + + return ret +end + +""" + max(T) + +Calculate the largest possible value +for the given IEEE 754 floating type. + +# Arguments +- `T`: a floating point number type <: AbstractFloat +""" +function max(T) + ret = one(T) + + while isfinite(ret * T(2)) + ret *= T(2) + end + + x = ret / T(2) + while isfinite(ret + x) + ret += x + x /= T(2) + end + + return ret +end + +for T in [Float16, Float32, Float64] + println("me - macheps($T): $(macheps(T))") + println("julia - eps($T): $(eps(T))") + + println() + + println("me - eta($T): $(eta(T))") + println("julia - nextfloat($T(0.0)): $(nextfloat(T(0.0)))") + + println() + + println("me - max($T): $(max(T))") + println("julia - floatmax($T): $(floatmax(T))") + + println() + + println("julia - floatmin($T): $(floatmin(T))") + + println() +end diff --git a/l1/2.jl b/l1/2.jl new file mode 100755 index 0000000..3c282a8 --- /dev/null +++ b/l1/2.jl @@ -0,0 +1,23 @@ +#!/usr/bin/env julia + +# Jacek Poziemski 272389 + +""" + kahanmacheps(T) + +Calculate the Kahan machine epsilon +using the following expression: `3(4/3 - 1) - 1`. + +# Arguments +- `T`: a floating point number type <: AbstractFloat +""" +function kahanmacheps(T) + return T(3) * ((T(4) / T(3)) - one(T)) - one(T) +end + +for T in [Float16, Float32, Float64] + println("kahanmacheps($T): $(kahanmacheps(T))") + println("eps($T): $(eps(T))") + + println() +end diff --git a/l1/3.jl b/l1/3.jl new file mode 100755 index 0000000..0a6b201 --- /dev/null +++ b/l1/3.jl @@ -0,0 +1,47 @@ +#!/usr/bin/env julia + +# Jacek Poziemski 272389 + +""" + exponentsequal(x::Float64, y::Float64) + +Check if the exponents of 2 64-bit IEEE 754 +numbers are equal. + +# Arguments +- `x::Float64`: first of the 2 numbers to check +- `y::Float64`: second of the 2 numbers to check +""" +function exponentsequal(x::Float64, y::Float64) + xexp = SubString(bitstring(x), 2:12) + yexp = SubString(bitstring(y), 2:12) + + return xexp == yexp +end + +""" + numberdistribution(start::Float64, finish::Float64) + +Calculate the step between 2 neighbouring +64-bit IEEE 754 numbers in a given range. + +# Arguments +- `start::Float64`: bottom bound of the range +- `finish::Float64`: upper bound of the range +""" +function numberdistribution(start::Float64, finish::Float64) + if !exponentsequal(nextfloat(start), prevfloat(finish)) + return missing + end + + exp = parse(Int, SubString(bitstring(start), 2:12), base = 2) + + return 2.0^(exp - 1023) * 2.0^(-52) +end + +println("[1.0, 2.0]: $(numberdistribution(1.0, 2.0))") +println("2^(-52): $(2.0^(-52))") +println("[0.5, 1.0]: $(numberdistribution(0.5, 1.0))") +println("2^(-53): $(2.0^(-53))") +println("[2.0, 4.0]: $(numberdistribution(2.0, 4.0))") +println("2^(-51): $(2.0^(-51))") diff --git a/l1/4.jl b/l1/4.jl new file mode 100755 index 0000000..17dec99 --- /dev/null +++ b/l1/4.jl @@ -0,0 +1,37 @@ +#!/usr/bin/env julia + +# Jacek Poziemski 272389 + +""" + findsmallest(start::Float64, finish::Float64) + +Find the smallest Float64 larger than `start` +and smaller than `finish` that satisfies the +following inequality `x * (1/x) ≠ 1`. + +# Arguments +- `start::Float64`: bottom limit for the wanted number +- `finish::Float64`: upper limit for the wanted number +""" +function findsmallest(start::Float64, finish::Float64) + x = nextfloat(Float64(start)) + res = zero(Float64) + + while x < Float64(finish) + res = x * (one(Float64) / x) + if res != one(Float64) + return x + end + x = nextfloat(x) + end + + return missing +end + +smallest = findsmallest(1.0, 2.0) + +if ismissing(smallest) + println("smallest not found") +else + println("smallest: $smallest") +end diff --git a/l1/5.jl b/l1/5.jl new file mode 100755 index 0000000..7e6a980 --- /dev/null +++ b/l1/5.jl @@ -0,0 +1,151 @@ +#!/usr/bin/env julia + +# Jacek Poziemski 272389 + +""" + sumvectorsforwards(x::Vector{<: AbstractFloat}, y::Vector{<: AbstractFloat}) + +Calculate the scalar sum of `x` and `y`, +going from the first element to the last. + +# Arguments +- `x::Vector{<: AbstractFloat}`: first of the 2 vectors +- `y::Vector{<: AbstractFloat}`: second of the 2 vectors +""" +function sumvectorsforwards(x::Vector{<: AbstractFloat}, y::Vector{<: AbstractFloat}) + s = 0.0 + + for i in 1:(length(x)) + s += x[i] * y[i] + end + + return s +end + +""" + sumvectorsbackwards(x::Vector{<: AbstractFloat}, y::Vector{<: AbstractFloat}) + +Calculate the scalar sum of `x` and `y`, +going from the last element to the first. + +# Arguments +- `x::Vector{<: AbstractFloat}`: first of the 2 vectors +- `y::Vector{<: AbstractFloat}`: second of the 2 vectors +""" +function sumvectorsbackwards(x::Vector{<: AbstractFloat}, y::Vector{<: AbstractFloat}) + s = 0.0 + + for i in (length(x)):-1:1 + s += x[i] * y[i] + end + + return s +end + +""" + sumvectorsdecreasing(x::Vector{<: AbstractFloat}, y::Vector{<: AbstractFloat}) + +Calculate the scalar sum of `x` and `y` by calculating +the partial sums of positive and negative numbers, +summing up the partial sums using a decreasing order +according to their absolute values. + +# Arguments +- `x::Vector{<: AbstractFloat}`: first of the 2 vectors +- `y::Vector{<: AbstractFloat}`: second of the 2 vectors +""" +function sumvectorsdecreasing(x::Vector{<: AbstractFloat}, y::Vector{<: AbstractFloat}) + s = zeros(length(x)) + + for i in 1:(length(x)) + s[i] = x[i] * y[i] + end + + spos = s[s .> 0] + sort!(spos, rev = true) + + sneg = s[s .<= 0] + sort!(sneg) + + neg = 0 + for i in sneg + neg += i + end + + pos = 0 + for i in spos + pos += i + end + + return pos + neg +end + +""" + sumvectorsincreasing(x::Vector{<: AbstractFloat}, y::Vector{<: AbstractFloat}) + +Calculate the scalar sum of `x` and `y` by calculating +the partial sums of positive and negative numbers, +summing up the partial sums using an increasing order +according to their absolute values. + +# Arguments +- `x::Vector{<: AbstractFloat}`: first of the 2 vectors +- `y::Vector{<: AbstractFloat}`: second of the 2 vectors +""" +function sumvectorsincreasing(x::Vector{<: AbstractFloat}, y::Vector{<: AbstractFloat}) + s = zeros(length(x)) + + for i in 1:(length(x)) + s[i] = x[i] * y[i] + end + + spos = s[s .> 0] + sort!(spos) + + sneg = s[s .<= 0] + sort!(sneg, rev = true) + + neg = 0 + for i in sneg + neg += i + end + + pos = 0 + for i in spos + pos += i + end + + return pos + neg +end + +x32::Vector{Float32} = [2.718281828, −3.141592654, 1.414213562, 0.5772156649, 0.3010299957] +y32::Vector{Float32} = [1486.2497, 878366.9879, −22.37492, 4773714.647, 0.000185049] + +x64::Vector{Float64} = [2.718281828, −3.141592654, 1.414213562, 0.5772156649, 0.3010299957] +y64::Vector{Float64} = [1486.2497, 878366.9879, −22.37492, 4773714.647, 0.000185049] + +expected = −1.00657107000000 * 10^(−11) + +res = sumvectorsforwards(x32, y32) +println("summing x32 and y32 forwards: $res; correct: $(res == expected)") + +res = sumvectorsbackwards(x32, y32) +println("summing x32 and y32 backwards: $res; correct: $(res == expected)") + +res = sumvectorsdecreasing(x32, y32) +println("summing x32 and y32 decreasing: $res; correct: $(res == expected)") + +res = sumvectorsincreasing(x32, y32) +println("summing x32 and y32 increasing: $res; correct: $(res == expected)") + +res = sumvectorsforwards(x64, y64) +println("summing x64 and y64 forwards: $res; correct: $(res == expected)") + +res = sumvectorsbackwards(x64, y64) +println("summing x64 and y64 backwards: $res; correct: $(res == expected)") + +res = sumvectorsdecreasing(x64, y64) +println("summing x64 and y64 decreasing: $res; correct: $(res == expected)") + +res = sumvectorsincreasing(x64, y64) +println("summing x64 and y64 increasing: $res; correct: $(res == expected)") diff --git a/l1/6.jl b/l1/6.jl new file mode 100755 index 0000000..172e40d --- /dev/null +++ b/l1/6.jl @@ -0,0 +1,33 @@ +#!/usr/bin/env julia + +# Jacek Poziemski 272389 + +""" + f(x::Float64) + +Calculate `sqrt(x^2 + 1) - 1`. + +# Arguments +- `x::Float64`: function argument +""" +function f(x::Float64) + return sqrt(x^2 + one(Float64)) - one(Float64) +end + +""" + g(x::Float64) + +Calculate `x^2 / (sqrt(x^2 + 1) + 1)`. + +# Arguments +- `x::Float64`: function argument +""" +function g(x::Float64) + return x^2 / (sqrt(x^2 + one(Float64)) + one(Float64)) +end + +for i in 1:20 + x::Float64 = 8.0^(-i) + println("f(8^-$i) = $(f(x))") + println("g(8^-$i) = $(g(x))") +end diff --git a/l1/7.jl b/l1/7.jl new file mode 100755 index 0000000..15a02f3 --- /dev/null +++ b/l1/7.jl @@ -0,0 +1,53 @@ +#!/usr/bin/env julia + +# Jacek Poziemski 272389 + +""" + f(x::Float64) + +Calculate `sin(x) + cos(3x)`. + +# Arguments +- `x::Float64`: function argument +""" +function f(x::Float64) + return sin(x) + cos(Float64(3) * x) +end + +""" + fderivative(x::Float64) + +Calculate the derivative of `sin(x) + cos(3x)`, `cos(x) - 3sin(3x)`. + +# Arguments +- `x::Float64`: function argument +""" +function fderivative(x::Float64) + return cos(x) - Float64(3) * sin(Float64(3) * x) +end + +""" + derivative(f::Function, x0::Float64, h::Float64) + +Calculate an approximation of the derivative of `f` at `x0`, +using `h` in the following formula: `(f(x0 + h) - f(x0)) / h`. + +# Arguments +- `f::Function`: function to calculate the derivative of +- `x0::Float64`: function argument +- `h::Float64`: value used in the approximation formula +""" +function derivative(f::Function, x0::Float64, h::Float64) + return (f(x0 + h) - f(x0)) / h +end + +x = one(Float64) + +for n in 0:54 + h = 2.0^(-n) + d = derivative(f, x, h) + e = abs(fderivative(x) - d) + println("n: $n; h: $h; h + 1: $(h + 1); derivative: $d; error: $e") +end + +println("exact value: $(fderivative(1.0))") diff --git a/l1/sprawozdanie.tex b/l1/sprawozdanie.tex new file mode 100644 index 0000000..de0a162 --- /dev/null +++ b/l1/sprawozdanie.tex @@ -0,0 +1,404 @@ +\documentclass[a4paper]{article} +\usepackage[T1]{fontenc} +\usepackage[polish]{babel} +\usepackage[utf8]{inputenc} +\usepackage{amsmath} +\usepackage{adjustbox} +\usepackage{longtable} +\usepackage{siunitx} +\title{Obliczenia naukowe - Lista 1} +\author{Jacek Poziemski 272389} +\date{27.10.2024} +\begin{document} +\maketitle + +\textbf{Zad. 1} + +Wyznaczanie iteracyjnie epsilonów maszynowych, liczb maszynowych $eta$ i maksymalnych wartości dla wszystkich dostępnych typów zmiennopozycyjnych zgodnych ze standardem IEEE 754 + +\vspace{0.5cm} + +\textbf{Epsilon maszynowy} + +\vspace{0.25cm} + +Epsilon maszynowy to najmniejsza wartość dla danego typu zmiennopozycyjnego, dla której spełniona jest następująca nierówność: $1.0 + \varepsilon > 1.0$. Innymi słowami, jest to wartość, która po dodaniu do liczby da nam następną liczbę po niej. Epsilon maszynowy jest odwrotnie proporcjonalny do precyzji typu zmiennopozycyjnego, w którym operujemy. + +\vspace{0.25cm} + +Epsilon liczę zaczynając od wartości $1.0$ dla danego typu i dzieląc ją przez $2.0$ dopóki $1.0 + \frac{\varepsilon}{2.0} > 1.0$. Podzielenie $\varepsilon$ przez $2.0$ w warunku pętli powoduje wyjście z pętli po otrzymaniu najmniejszej wartości, dla której warunek jest spełniony. + +\vspace{0.25cm} + +Porównanie wartości $\varepsilon$ liczonych iteracyjnie, za pomocą funkcji $eps()$ z Julii oraz z headera $\langle float.h \rangle$ biblioteki standardowej języka C: + +\begin{table}[h] + \begin{adjustbox}{max width=\textwidth} + \begin{tabular}{|c|c|c|c|} + \hline + Typ zmiennopozycyjny & Wartość wyznaczona iteracyjnie & $eps()$ & $\langle float.h \rangle$ \\ \hline + Float16 & 0.000977 & 0.000977 & — \\ + Float32 & 1.1920929e-7 & 1.1920929e-7 & 1.1920929e-07 \\ + Float64 & 2.220446049250313e-16 & 2.220446049250313e-16 & 2.220446049250313e-16 \\ + \hline + \end{tabular} + \end{adjustbox} +\end{table} + +Brak wartości Float16 dla headera $\langle float.h \rangle$ wynika z braku 16-bitowego typu zmiennopozycyjnego w standardzie języka C. + +\vspace{0.25cm} + +Wartości wyliczone iteracyjnie są identyczne do tych z bibliotek standardowych języków Julia i C, więc jeśli ufamy twórcom tych języków, metoda iteracyjna jest wystarczająco dokładna. + +\vspace{0.5cm} + +\textbf{Liczba maszynowa eta} + +\vspace{0.25cm} + +Liczba maszynowa eta to najmniejsza zdenormalizowana wartość dla danego typu zmiennopozycyjnego większa od $0.0$. + +\vspace{0.25cm} + +Liczbę eta liczę podobnie do liczby maszynowej epsilon – zaczynam od wartości $1.0$ i dzielę ją przez $2.0$ dopóki warunek $\frac{\eta}{2.0} > 0.0$ jest spełniony. Tak samo jak poprzednio podzielenie $\eta$ przez $2.0$ powoduje wyjście z pętli po otrzymaniu najmniejszej wartości, dla której warunek jest spełniony. + +Porównanie wartości $\eta$ liczonych iteracyjnie, za pomocą $nextfloat(0.0)$ z Julii oraz wartościami zwróconymi przez funkcję $floatmin()$ dla danego typu zmiennopozycyjnego: + +\begin{table}[h] + \begin{adjustbox}{max width=\textwidth} + \begin{tabular}{|c|c|c|c|} + \hline + Typ zmiennopozycyjny & Wartość wyznaczona iteracyjnie & $nextfloat(0.0)$ & $floatmin()$ \\ \hline + Float16 & 6.0e-8 & 6.0e-8 & 6.104e-5 \\ + Float32 & 1.0e-45 & 1.0e-45 & 1.1754944e-38 \\ + Float64 & 5.0e-324 & 5.0e-324 & 2.2250738585072014e-308 \\ + \hline + \end{tabular} + \end{adjustbox} +\end{table} + +Wartości wyliczone iteracyjnie są identyczne do tych zwróconych przez funkcję $nextfloat(0.0)$, ale $floatmin()$ zwraca wartości o wiele większe. Jest tak dlatego, że ta funkcja zwraca wartości znormalizowane, a metoda iteracyjna i funkcja $nextfloat(0.0)$ zwracają wartości zdenormalizowane - takie, których cecha składa się z samych bitów 0. + +\vspace{0.5cm} + +\textbf{Maksymalne wartości} + +\vspace{0.5cm} + +Wartość maksymalną dla danego typu zmiennopozycyjnego liczę zaczynając od $1.0$ dla tego typu. Następnie mnożę tą wartość przez $2.0$ dopóki $max \cdot 2.0 < \infty$ (sprawdzam tę nierówność korzystając z funkcji $isfinite()$). Osiągnięta w ten sposób wartość niekoniecznie jest maksymalna, żeby to zmienić tworzę pomocniczą zmienną $x$, której na początku przypisuję wartość $\frac{max}{2.0}$. Następnie, dopóki $max + x < \infty$ dodaję $x$ do $max$ i dzielę $x$ przez $2.0$. + +\vspace{0.25cm} + +W przypadku gdyby w ostatniej iteracji pierwszej pętli $max$ nie był maksymalny, ale $max \cdot 2.0$ był nieskończony, zostało nam co najwyżej $max \cdot 2.0 - \varepsilon$ do "`dopełnienia"'. Zaczynamy to "`dopełnianie"' od $\frac{max}{2.0}$ – wiemy, że $max + max = max \cdot 2.0$ będzie nieskończone. W każdej iteracji dodajemy $x$ do $max$ i dzielimy $x$ przez $2.0$. W ten sposób "`dopełnimy"' $max$ do końca i otrzymamy rzeczywistą wartość maksymalną. + +\vspace{0.25cm} + +Porównanie wartości maksymalnych liczonych iteracyjnie, za pomocą funkcji $floatmax()$ z Julii oraz z headera $\langle float.h \rangle$ biblioteki standardowej języka C: + +\begin{table}[h] + \begin{adjustbox}{max width=\textwidth} + \begin{tabular}{|c|c|c|c|} + \hline + Typ zmiennopozycyjny & Wartość wyznaczona iteracyjnie & $floatmax()$ & $\langle float.h \rangle$ \\ \hline + Float16 & 6.55e4 & 6.55e4 & — \\ + Float32 & 3.4028235e38 & 3.4028235e38 & 3.4028235e+38 \\ + Float64 & 1.7976931348623157e308 & 1.7976931348623157e308 & 1.7976931348623157e+308 \\ + \hline + \end{tabular} + \end{adjustbox} +\end{table} + +Brak wartości Float16 dla headera $\langle float.h \rangle$ wynika z braku 16-bitowego typu zmiennopozycyjnego w standardzie języka C. + +\vspace{0.25cm} + +Wartości wyznaczone iteracyjnie są identyczne do tych zwróconych przez $floatmax()$ jak i tych z headera $\langle float.h \rangle$ z biblioteki standardowej języka C. + +\vspace{0.5cm} + +\textbf{Zad. 2} + +Eksperymentalne sprawdzenie słuszności twierdzenia Kahana dla wszystkich typów zmiennopozycyjnych dostępnych w języku Julia. + +\vspace{0.25cm} + +Twierdzenie Kahana mówi, że epsilon maszynowy można uzyskać w następujący sposób: $\varepsilon = 3(4/3 - 1) / 1$. Aby sprawdzić wyniki korzystam z funkcji $eps()$ z biblioteki standardowej Julii. + +\vspace{0.25cm} + +Porównanie wartości epsilonu maszynowego liczonego za pomocą twierdzenia Kahana oraz za pomocą funkcji $eps()$: + +\begin{table}[h] + \begin{adjustbox}{max width=\textwidth} + \begin{tabular}{|c|c|c|} + \hline + Typ zmiennopozycyjny & Twierdzenie Kahana & $eps()$ \\ \hline + Float16 & -0.000977 & 0.000977 \\ + Float32 & 1.1920929e-7 & 1.1920929e-7 \\ + Float64 & -2.220446049250313e-16 & 2.220446049250313e-16 \\ + \hline + \end{tabular} + \end{adjustbox} +\end{table} + +Dla typu $Float32$ wartości są identyczne, dla typów $Float16$ i $Float64$ wartości różnią się znakiem - twierdzenie Kahana daje nam wartości ujemne. Wynika to z definicji typów zmiennoprzecinkowych i rozwinięcia liczby $4/3$ binarnie – $1.(10)$. Typy $Float16$, $Float32$ i $Float64$ mają odpowiednio 10, 23 i 52 bity mantysy. Rozwijając $4/3$, na indeksach parzystych dostajemy bit 0, na nieparzystych bit 1. W typach o parzystej ilości bitów w mantysie ($Float16$ i $Float64$) na ostatnim indeksie będzie 0. Mnożąc wartość $4/3 - 1$ przez $3$ otrzymamy wartość minimalnie mniejszą od $1$ Przez co wynik jest ujemny. Z kolei dla typu $Float32$ ostatnim bitem będzie 1, co da nam wartość dodatnią. + +\vspace{0.25cm} + +Twierdzenie Kahana poprawnie liczy epsilon maszynowy jednak trzeba pamiętać o wzięciu modułu z wyniku. + +\vspace{0.5cm} + +\textbf{Zad. 3} + +Eksperymentalne sprawdzenie rozmieszczenia liczb zmiennopozycyjnych w arytmetyce $Float64$ w języku Julia. + +\vspace{0.25cm} + +Aby sprawdzić, czy liczby w przedziale $[1, 2]$ są równomiernie rozmieszczone można zrobić pętlę, w której od $1.0$ do $2.0$ jedną liczbę iterujemy dodając do niej krok $\delta = 2^{-52}$, drugą korzystając z funkcji $nextfloat()$ jednak jest to bardzo nieefektywna metoda. + +\vspace{0.25cm} + +Lepszą metodą jest sprawdzenie pierwszej i ostatniej liczby z przedziału, w naszym przypadku $nextfloat(1.0)$ i $nextfloat(2.0)$. Sprawdzamy cechę obu – jeśli są różne, wyklucza to równomierne rozmieszczenie liczb. W przeciwnym przypadku możemy obliczyć $\delta$ w następujący sposób: + +\[2^{cecha - 1023} \cdot 2^{-52}\] + +gdzie 1023 to przesunięcie cechy typu $Float64$, a 52 to rozmiar mantysy. + +\begin{table}[h] + \begin{adjustbox}{max width=\textwidth} + $\begin{array}{|c|c|} + \hline + \text{Przedział} & \delta \\ \hline + \text{[1.0, 2.0]} & \num{2.220446049250313e-16} = 2^{-52} \\ + \text{[0.5, 1.0]} & \num{1.1102230246251565e-16} = 2^{-53} \\ + \text{[2.0, 4.0]} & \num{4.440892098500626e-16} = 2^{-51} \\ + \hline + \end{array}$ + \end{adjustbox} +\end{table} + +Dla przedziału [1.0, 2.0] rzeczywiście $\delta = 2^{-52}$, z kolei $\delta$ dla przedziałów [0.5, 1.0] i [2.0, 4.0] pokazuje, że jest ona zależna od cechy liczb z zakresu nad którym się skupiamy. Im bliżej zera jesteśmy tym mniejsza jest $\delta$ – zwiększa się dokładność liczb. Analogicznie oddalając się od zera tracimy precyzję części ułamkowej – przy coraz większych liczbach prawdopodobnie bardziej obchodzi nas część całkowita. + +\vspace{0.5cm} + +\textbf{Zad. 4} + +Eksperymentalne znalezienie liczby w arytmetyce $Float64$ w przedziale $(1, 2)$ takiej, że $x \cdot (1/x) \neq 1$ oraz znalezienie najmniejszej takiej liczby. + +\vspace{0.25cm} + +Implementacja zadania opiera się na funkcji $nextfloat()$ z języka Julia. Zaczynamy od $nextfloat(1.0)$, iterujemy do $2.0$. W każdej iteracji liczymy wynik $x \cdot (1/x)$, jeśli jest różny od $1.0$ zwracamy $x$ jako wynik. Oczywiście jest to również najmniejsza taka liczba. + +\vspace{0.25cm} + +Najmniejsza taka liczba to $1.000000057228997$. Błędny wynik wyrażenia, które zawsze powinno zwracać $1.0$ jest spowodowany oczywiście skończoną mantysą w liczbach zmiennoprzecinkowych – przez co tracimy precyzję – ale również kolejnością wykonywania działań. Gdyby zamienić $x \cdot (1/x)$ na $(x \cdot 1) / x$ zawsze dostalibyśmy odpowiedni wynik – $x \cdot 1$ nie zmieniłoby wartości liczby, a $x / x$ zakładając dobrze napisany interpreter lub kompilator zawsze zwróci $1.0$ + +\vspace{0.5cm} + +\textbf{Zad. 5} + +Eksperymentalne obliczanie iloczynu skalarnego dwóch wektorów: + +\begin{gather*} + x = [2.718281828, -3.141592654, 1.414213562, 0.5772156649, 0.3010299957] \\ + y = [1486.2497, 878366.9879, -22.37492, 4773714.647, 0.000185049] +\end{gather*} + +na 4 sposoby: + +\vspace{0.25cm} + +(a) w przód – od pierwszego do ostatniego elementu + +\vspace{0.25cm} + +(b) w tył – od ostatniego do pierwszego elementu + +\vspace{0.25cm} + +(c) od największego do najmniejszego – liczenie częściowych sum elementów dodatnich i ujemnych w porządku malejącym względem wartości absolutnej liczb, następnie dodanie do siebie sum częściowych + +\vspace{0.25cm} + +(d) od największego do najmniejszego – liczenie częściowych sum elementów dodatnich i ujemnych w porządku rosnącym względem wartości absolutnej liczb, następnie dodanie do siebie sum częściowych + +\vspace{0.25cm} + +używając typów zmiennopozycyjnych $Float32$ i $Float64$. + +\vspace{0.25cm} + +Oczekiwanym wynikiem iloczynu skalarnego na tych dwóch wektorach jest $-1.00657107000000 \cdot 10^{-11}$. + +\vspace{0.25cm} + +Porównanie wyników iloczynu skalarnego przeprowadzonego za pomocą powyższych czterech sposobów: + +\begin{table}[h] + \begin{adjustbox}{max width=\textwidth} + \begin{tabular}{|c|c|c|} + \hline + Sposób & Wynik dla $Float32$ & Wynik dla $Float64$ \\ \hline + (a) & -0.3472038161853561 & 1.0251881368296672e-10 \\ + (b) & -0.3472038162872195 & -1.5643308870494366e-10 \\ + (c) & -0.3472038162872195 & 0.0 \\ + (d) & -0.3472038162872195 & 0.0 \\ + \hline + \end{tabular} + \end{adjustbox} +\end{table} + +Żaden ze sposobów nie dał wyniku bliskiego oczekiwanemu, choć sposób (b) w arytmetyce $Float64$ dał wynik dość bliski prawdzie. Patrząc na wyniki korzystające z różnych sposobów widzimy również, że kolejność wykonywania działań może znacznie wpłynąć na poprawność wyniku. + +\vspace{0.5cm} + +\textbf{Zad. 6} + +Policzenie w języku Julia w arytmetyce $Float64$ wartości funkcji: + +\[f(x) = \sqrt{x^2 + 1} - 1\] + +\[g(x) = x^2 / (\sqrt{x^2 + 1} + 1)\] + +dla argumentów $x = 8^{-1}, 8^{-2}, 8^{-3}, ...$. + +\vspace{0.25cm} + +Porównanie wyników dla obu funkcji: + +\newpage + +\begin{table}[h] + \begin{adjustbox}{max width=\textwidth} + \begin{tabular}{|c|c|c|} + \hline + $x$ & $f(x)$ & $g(x)$ \\ \hline + $8^{-1}$ & 0.0077822185373186414 & 0.0077822185373187065 \\ + $8^{-2}$ & 0.00012206286282867573 & 0.00012206286282875901 \\ + $8^{-3}$ & 1.9073468138230965e-6 & 1.907346813826566e-6 \\ + $8^{-4}$ & 2.9802321943606103e-8 & 2.9802321943606116e-8 \\ + $8^{-5}$ & 4.656612873077393e-10 & 4.6566128719931904e-10 \\ + $8^{-6}$ & 7.275957614183426e-12 & 7.275957614156956e-12 \\ + $8^{-7}$ & 1.1368683772161603e-13 & 1.1368683772160957e-13 \\ + $8^{-8}$ & 1.7763568394002505e-15 & 1.7763568394002489e-15 \\ + $8^{-9}$ & 0.0 & 2.7755575615628914e-17 \\ + $8^{-10}$ & 0.0 & 4.336808689942018e-19 \\ + $8^{-11}$ & 0.0 & 6.776263578034403e-21 \\ + $8^{-12}$ & 0.0 & 1.0587911840678754e-22 \\ + $8^{-13}$ & 0.0 & 1.6543612251060553e-24 \\ + $8^{-14}$ & 0.0 & 2.5849394142282115e-26 \\ + $8^{-15}$ & 0.0 & 4.0389678347315804e-28 \\ + $8^{-16}$ & 0.0 & 6.310887241768095e-30 \\ + $8^{-17}$ & 0.0 & 9.860761315262648e-32 \\ + $8^{-18}$ & 0.0 & 1.5407439555097887e-33 \\ + $8^{-19}$ & 0.0 & 2.407412430484045e-35 \\ + $8^{-20}$ & 0.0 & 3.76158192263132e-37 \\ + \hline + \end{tabular} + \end{adjustbox} +\end{table} + +Już przy $x = 8^{-9}$, $f(x)$ przyjmuje wartość $0.0$, którą zachowuje przy mniejszych wartościach. Przy wartościach większych od $8^{-9}$ obie funkcje mają podobne wartości. $f(x)$ szybciej traci precyzję i spada do $0.0$ przez odejmowanie $1.0$ od pierwiastka – im mniejsza wartość $x$ tym bliżej $x^2 + 1$ jest $1.0$. Przez pierwiastek całe wyrażenie zbliża się jeszcze bardziej do $1.0$, od którego również odejmujemy $1.0$. Przekształcenie $f(x)$ w $g(x)$ zwalcza ten problem poprzez pozbycie się problematycznego odejmowania, dzięki czemu tracimy mniej precyzji wykonując działania. + +\vspace{0.5cm} + +\textbf{Zad. 7} + +Obliczanie przybliżonej wartości pochodnej funkcji $f(x) = \sin{x} + \cos{3x}$ w punkcie $x_0 = 1$ oraz błędów $|f^{'}(x_0) - \widetilde{f^{'}}(x_0)|$ korzystając ze wzoru + +\[f^{'}(x_0) \approx \widetilde{f^{'}}(x_0) = \frac{f(x_0 + h) - f(x_0)}{h}\] + +dla $h = 2^{-n}$, gdzie $n = 0, 1, 2, ..., 54$. + +\vspace{0.25cm} + +Dokładna pochodna $f(x)$ to $f^{'}(x) = \cos{x} - 3\sin{3x}$ + +\vspace{0.25cm} + +Dokładna wartość pochodnej dla $x_0 = 1$ to $0.11694228168853815$. + +\begin{longtable}{|c|c|c|c|c|} + \hline + $n$ & $h$ & $1 + h$ & $f^{'}(x_0)$ & błąd \\ \hline + $0$ & 1.0 & 2.0 & 2.0179892252685967 & 1.9010469435800585 \\ + $1$ & 0.5 & 1.5 & 1.8704413979316472 & 1.753499116243109 \\ + $2$ & 0.25 & 1.25 & 1.1077870952342974 & 0.9908448135457593 \\ + $3$ & 0.125 & 1.125 & 0.6232412792975817 & 0.5062989976090435 \\ + $4$ & 0.0625 & 1.0625 & 0.3704000662035192 & 0.253457784514981 \\ + $5$ & 0.03125 & 1.03125 & 0.24344307439754687 & 0.1265007927090087 \\ + $6$ & 0.015625 & 1.015625 & 0.18009756330732785 & 0.0631552816187897 \\ + $7$ & 0.0078125 & 1.0078125 & 0.1484913953710958 & 0.03154911368255764 \\ + $8$ & 0.00390625 & 1.00390625 & 0.1327091142805159 & 0.015766832591977753 \\ + $9$ & 0.001953125 & 1.001953125 & 0.1248236929407085 & 0.007881411252170345 \\ + $10$ & 0.0009765625 & 1.0009765625 & 0.12088247681106168 & 0.0039401951225235265 \\ + $11$ & 0.00048828125 & 1.00048828125 & 0.11891225046883847 & 0.001969968780300313 \\ + $12$ & 0.000244140625 & 1.000244140625 & 0.11792723373901026 & 0.0009849520504721099 \\ + $13$ & 0.0001220703125 & 1.0001220703125 & 0.11743474961076572 & 0.0004924679222275685 \\ + $14$ & 6.103515625e-5 & 1.00006103515625 & 0.11718851362093119 & 0.0002462319323930373 \\ + $15$ & 3.0517578125e-5 & 1.000030517578125 & 0.11706539714577957 & 0.00012311545724141837 \\ + $16$ & 1.52587890625e-5 & 1.0000152587890625 & 0.11700383928837255 & 6.155759983439424e-5 \\ + $17$ & 7.62939453125e-6 & 1.0000076293945312 & 0.11697306045971345 & 3.077877117529937e-5 \\ + $18$ & 3.814697265625e-6 & 1.0000038146972656 & 0.11695767106721178 & 1.5389378673624776e-5 \\ + $19$ & 1.9073486328125e-6 & 1.0000019073486328 & 0.11694997636368498 & 7.694675146829866e-6 \\ + $20$ & 9.5367431640625e-7 & 1.0000009536743164 & 0.11694612901192158 & 3.8473233834324105e-6 \\ + $21$ & 4.76837158203125e-7 & 1.0000004768371582 & 0.1169442052487284 & 1.9235601902423127e-6 \\ + $22$ & 2.384185791015625e-7 & 1.000000238418579 & 0.11694324295967817 & 9.612711400208696e-7 \\ + $23$ & 1.1920928955078125e-7 & 1.0000001192092896 & 0.11694276239722967 & 4.807086915192826e-7 \\ + $24$ & 5.960464477539063e-8 & 1.0000000596046448 & 0.11694252118468285 & 2.394961446938737e-7 \\ + $25$ & 2.9802322387695312e-8 & 1.0000000298023224 & 0.116942398250103 & 1.1656156484463054e-7 \\ + $26$ & 1.4901161193847656e-8 & 1.0000000149011612 & 0.11694233864545822 & 5.6956920069239914e-8 \\ + $27$ & 7.450580596923828e-9 & 1.0000000074505806 & 0.11694231629371643 & 3.460517827846843e-8 \\ + $28$ & 3.725290298461914e-9 & 1.0000000037252903 & 0.11694228649139404 & 4.802855890773117e-9 \\ + $29$ & 1.862645149230957e-9 & 1.0000000018626451 & 0.11694222688674927 & 5.480178888461751e-8 \\ + $30$ & 9.313225746154785e-10 & 1.0000000009313226 & 0.11694216728210449 & 1.1440643366000813e-7 \\ + $31$ & 4.656612873077393e-10 & 1.0000000004656613 & 0.11694216728210449 & 1.1440643366000813e-7 \\ + $32$ & 2.3283064365386963e-10 & 1.0000000002328306 & 0.11694192886352539 & 3.5282501276157063e-7 \\ + $33$ & 1.1641532182693481e-10 & 1.0000000001164153 & 0.11694145202636719 & 8.296621709646956e-7 \\ + $34$ & 5.820766091346741e-11 & 1.0000000000582077 & 0.11694145202636719 & 8.296621709646956e-7 \\ + $35$ & 2.9103830456733704e-11 & 1.0000000000291038 & 0.11693954467773438 & 2.7370108037771956e-6 \\ + $36$ & 1.4551915228366852e-11 & 1.000000000014552 & 0.116943359375 & 1.0776864618478044e-6 \\ + $37$ & 7.275957614183426e-12 & 1.000000000007276 & 0.1169281005859375 & 1.4181102600652196e-5 \\ + $38$ & 3.637978807091713e-12 & 1.000000000003638 & 0.116943359375 & 1.0776864618478044e-6 \\ + $39$ & 1.8189894035458565e-12 & 1.000000000001819 & 0.11688232421875 & 5.9957469788152196e-5 \\ + $40$ & 9.094947017729282e-13 & 1.0000000000009095 & 0.1168212890625 & 0.0001209926260381522 \\ + $41$ & 4.547473508864641e-13 & 1.0000000000004547 & 0.116943359375 & 1.0776864618478044e-6 \\ + $42$ & 2.2737367544323206e-13 & 1.0000000000002274 & 0.11669921875 & 0.0002430629385381522 \\ + $43$ & 1.1368683772161603e-13 & 1.0000000000001137 & 0.1162109375 & 0.0007313441885381522 \\ + $44$ & 5.684341886080802e-14 & 1.0000000000000568 & 0.1171875 & 0.0002452183114618478 \\ + $45$ & 2.842170943040401e-14 & 1.0000000000000284 & 0.11328125 & 0.003661031688538152 \\ + $46$ & 1.4210854715202004e-14 & 1.0000000000000142 & 0.109375 & 0.007567281688538152 \\ + $47$ & 7.105427357601002e-15 & 1.000000000000007 & 0.109375 & 0.007567281688538152 \\ + $48$ & 3.552713678800501e-15 & 1.0000000000000036 & 0.09375 & 0.023192281688538152 \\ + $49$ & 1.7763568394002505e-15 & 1.0000000000000018 & 0.125 & 0.008057718311461848 \\ + $50$ & 8.881784197001252e-16 & 1.0000000000000009 & 0.0 & 0.11694228168853815 \\ + $51$ & 4.440892098500626e-16 & 1.0000000000000004 & 0.0 & 0.11694228168853815 \\ + $52$ & 2.220446049250313e-16 & 1.0000000000000002 & -0.5 & 0.6169422816885382 \\ + $53$ & 1.1102230246251565e-16 & 1.0 & 0.0 & 0.11694228168853815 \\ + $54$ & 5.551115123125783e-17 & 1.0 & 0.0 & 0.11694228168853815 \\ + \hline +\end{longtable} + +Najdokładniejszy wynik jest dla $n = 28$, $h = 2^{-28}$, dokładnie po środku wartości testowanych (jeśli zignorujemy 54, dla którego wynik jest identyczny do 53). Idąc w górę lub w dół błąd jest coraz większy, dla $n = 53$ dochodząc do wartości pochodnej. + +\vspace{0.25cm} + +Wzięcie zbyt dużego $h$ (zbyt małego $n$) powoduje ogromny błąd z definicji pochodnej, którą jest granica wzoru wykorzystanego w tym zadaniu przy h dążącym do 0. Przy coraz mniejszych $h$ (coraz większych $n$) na papierze zbliżamy się do poprawnej wartości, ale w rzeczywistości ograniczona precyzja liczb zmiennopozycyjnych w standardzie IEEE 754 powoduje zwiększanie się błędu, lecz nie do takich rozmiarów jak dla pierwszych kilku $n$. Maksymalny błąd po przekroczeniu $n = 28$ otrzymujemy przy $n = 53$, gdzie $h < \varepsilon$ w precyzji $Float64$. Wtedy + +\[x_0 + h = x_0 = 1\] + +a we wzorze na pochodną + +\[\frac{f(x_0 + h) - f(x_0)}{h} = \frac{f(x_0) - f(x_0)}{h} = \frac{0}{h} = 0\] + +zatem dla $n \geq 53$ błąd będzie zawsze równy dokładnej wartości pochodnej. + +\vspace{0.5cm} + +\textbf{Wnioski z listy} + +Liczby zmiennopozycyjne w standardzie IEEE 754 mają skończoną precyzję, co przy liczbach bliskich zeru powoduje błędy. Pisząc kod operujący na tych liczbach musimy pamiętać o tym fakcie i wybrać odpowiedni typ do trzymania danych lub zastanowić się nad trzymaniem ich w typach całkowitych, np. waluty zamiast we floacie trzymać w dwóch integerach. + +\end{document} diff --git a/l1/values.c b/l1/values.c new file mode 100644 index 0000000..72e255d --- /dev/null +++ b/l1/values.c @@ -0,0 +1,16 @@ +#include +#include + +#define PRINT(var, format) printf(#var " = " format "\n", var) + +// https://en.cppreference.com/w/c/language/arithmetic_types#Real_floating_types +// https://en.cppreference.com/w/c/types/limits#Limits_of_floating-point_types +int main(void) { + PRINT(FLT_EPSILON, "%.9g"); + PRINT(DBL_EPSILON, "%.16lg"); + + PRINT(FLT_MAX, "%.8g"); + PRINT(DBL_MAX, "%.17lg"); + + return 0; +}