l1

.envrc

@@ -0,0 +1 @@
+use flake

.gitignore

@@ -0,0 +1,4 @@
+**/*.aux
+**/*.log
+**/*.pdf
+.direnv/

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
+}

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
+          ];
+        };
+      }
+    );
+  };
+}

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

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

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

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

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

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

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

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}

l1/values.c

@@ -0,0 +1,16 @@
+#include <float.h>
+#include <stdio.h>
+
+#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;
+}