718 lines
21 KiB
Plaintext
718 lines
21 KiB
Plaintext
# Maschinenzahlen
|
||
|
||
|
||
|
||
```{julia}
|
||
for x ∈ ( 3, 3.3e4, Int16(20), Float32(3.3e4), UInt16(9) )
|
||
@show x sizeof(x) typeof(x)
|
||
println()
|
||
end
|
||
```
|
||
|
||
## Ganze Zahlen *(integers)*
|
||
|
||
Ganze Zahlen werden grundsätzlich als Bitmuster fester Länge gespeichert. Damit ist der Wertebereich endlich.
|
||
**Innerhalb dieses Wertebereichs** sind Addition, Subtraktion, Multiplikation und die Ganzzahldivision mit Rest
|
||
exakte Operationen ohne Rundungsfehler.
|
||
|
||
Ganze Zahlen gibt es in den zwei Sorten `Signed` (mit Vorzeichen) und `Unsigned`, die man als Maschinenmodelle für ℤ bzw. ℕ auffassen kann.
|
||
|
||
|
||
### *Unsigned integers*
|
||
```{julia}
|
||
subtypes(Unsigned)
|
||
```
|
||
UInts sind Binärzahlen mit n=8, 16, 32, 64 oder 128 Bits Länge und einem entsprechenden Wertebereich von
|
||
$$
|
||
0 \le x < 2^n
|
||
$$
|
||
Sie werden im *scientific computing* eher selten verwendet. Bei hardwarenaher Programmierung dienen sie z.B. dem Umgang mit Binärdaten und Speicheradressen. Deshalb werden sie von Julia standardmäßig auch als Hexadezimalzahlen (mit Präfix `0x` und Ziffern `0-9a-f`) dargestellt.
|
||
|
||
```{julia}
|
||
x = 0x0033efef
|
||
@show x typeof(x) Int(x)
|
||
|
||
z = UInt(32)
|
||
@show z typeof(z);
|
||
|
||
```
|
||
|
||
### *Signed Integers*
|
||
```{julia}
|
||
subtypes(Signed)
|
||
```
|
||
*Integers* haben den Wertebereich
|
||
$$
|
||
-2^{n-1} \le x < 2^{n-1}
|
||
$$
|
||
|
||
|
||
In Julia sind ganze Zahlen standardmäßig 64 Bit groß:
|
||
```{julia}
|
||
x = 42
|
||
typeof(x)
|
||
```
|
||
Sie haben daher den Wertebereich:
|
||
$$
|
||
-9.223.372.036.854.775.808 \le x \le 9.223.372.036.854.775.807
|
||
$$
|
||
|
||
Negative Zahlen werden im Zweierkomplement dargestellt:
|
||
|
||
$x \Rightarrow -x$ entspricht: _flip all bits, then add 1_
|
||
|
||
Das sieht so aus:
|
||
|
||
::: {.content-visible when-format="html"}
|
||
{width=50%}
|
||
:::
|
||
|
||
::: {.content-visible when-format="pdf"}
|
||
{width=50%}
|
||
:::
|
||
|
||
Das höchstwertige Bit zeigt das Vorzeichen an. Sein „Umkippen“ entspricht allerdings nicht der Operation `x → -x`.
|
||
|
||
:::{.callout-important}
|
||
Die gesamte Ganzzahlarithmetik ist eine __Arithmetik modulo $2^n$__
|
||
|
||
```{julia}
|
||
x = 2^62 - 10 + 2^62
|
||
```
|
||
```{julia}
|
||
x + 20
|
||
```
|
||
Keine Fehlermeldung, keine Warnung! Ganzzahlen fester Länge liegen nicht auf einer Geraden, sondern auf einem Kreis!
|
||
|
||
:::
|
||
|
||
Schauen wir uns ein paar *Integers* mal an:
|
||
|
||
|
||
```{julia}
|
||
using Printf
|
||
|
||
for i in (1, 2, 7, -7, 2^50, -1, Int16(7), Int8(7))
|
||
s = bitstring(i)
|
||
t = typeof(i)
|
||
@printf "%20i %-5s ⇒ %s\n" i t s
|
||
end
|
||
```
|
||
|
||
... und noch etwas mehr *introspection*:
|
||
|
||
|
||
```{julia}
|
||
typemin(Int64), typemax(Int64)
|
||
```
|
||
|
||
```{julia}
|
||
typemin(UInt64), typemax(UInt64)
|
||
```
|
||
|
||
```{julia}
|
||
typemin(Int8), typemax(Int8)
|
||
```
|
||
|
||
## Arithmetik ganzer Zahlen
|
||
|
||
#### Addition, Multiplikation
|
||
Die Operationen `+`,`-`,`*` haben die übliche exakte Arithmetik **modulo $2^{64}$**.
|
||
|
||
#### Potenzen `a^b`
|
||
|
||
- Potenzen `a^n` werden für natürliche Exponenten `n` ebenfalls modulo $2^{64}$ exakt berechnet.
|
||
- Für negative Exponenten ist das Ergebnis eine Gleitkommazahl.
|
||
- `0^0` ist [selbstverständlich](https://en.wikipedia.org/wiki/Zero_to_the_power_of_zero#cite_note-T4n3B-4) gleich 1.
|
||
```{julia}
|
||
(-2)^63, 2^64, 3^(-3), 0^0
|
||
```
|
||
- Für natürliche Exponenten wird [*exponentiation by squaring*](https://de.wikipedia.org/wiki/Bin%C3%A4re_Exponentiation) verwendet, so dass z.B. `x^23` nur 7 Multiplikationen benötigt:
|
||
$$
|
||
x^{23} = \left( \left( (x^2)^2 \cdot x \right)^2 \cdot x \right)^2 \cdot x
|
||
$$
|
||
|
||
#### Division
|
||
|
||
- Division `/` erzeugt eine Gleitkommazahl.
|
||
|
||
```{julia}
|
||
x = 40/5
|
||
```
|
||
#### Ganzzahldivision und Rest
|
||
|
||
- Die Funktionen `div(a,b)`, `rem(a,b)` und `divrem(a,b)` berechnen den Quotient der Ganzzahldivision, den dazugehörigen Rest *(remainder)* bzw. beides als Tupel.
|
||
- Für `div(a,b)` gibt es die Operatorform `a ÷ b` (Eingabe: `\div<TAB>`) und für `rem(a,b)` die Operatorform `a % b`.
|
||
- Standardmäßig wird bei der Division „zur Null hin gerundet“, wodurch der dazugehörige Rest dasselbe Vorzeichen wie der Dividend `a` trägt:
|
||
|
||
```{julia}
|
||
@show divrem( 27, 4)
|
||
@show ( 27 ÷ 4, 27 % 4)
|
||
@show (-27 ÷ 4, -27 % 4 )
|
||
@show ( 27 ÷ -4, 27 % -4);
|
||
```
|
||
|
||
- Eine von `RoundToZero` abweichende Rundungsregel kann bei den Funktionen als optionales 3. Argument angegeben werden.
|
||
- `?RoundingMode` zeigt die möglichen Rundungsregeln an.
|
||
- Für die Rundungsregel `RoundDown` („in Richtung minus unendlich"), wodurch der dazugehörige Rest dasselbe Vorzeichen wie der Divisor `b` bekommt, gibt es auch die Funktionen `fld(a,b)` *(floored division)* und `mod(a,b)`:
|
||
|
||
```{julia}
|
||
@show divrem(-27, 4, RoundDown)
|
||
@show (fld(-27, 4), mod(-27, 4))
|
||
@show (fld( 27, -4), mod( 27, -4));
|
||
```
|
||
|
||
Für alle Rundungsregeln gilt:
|
||
```
|
||
div(a, b, RoundingMode) * b + rem(a, b, RoundingMode) = a
|
||
```
|
||
|
||
#### Der Datentyp `BigInt`
|
||
|
||
Der Datentyp `BigInt` ermöglicht Ganzzahlen beliebiger Länge. Der benötigte Speicher wird dynamisch allokiert.
|
||
|
||
Numerische Konstanten haben automatisch einen ausreichend großen Typ:
|
||
|
||
```{julia}
|
||
z = 10
|
||
@show typeof(z)
|
||
z = 10_000_000_000_000_000 # 10 Billiarden
|
||
@show typeof(z)
|
||
z = 10_000_000_000_000_000_000 # 10 Trillionen
|
||
@show typeof(z)
|
||
z = 10_000_000_000_000_000_000_000_000_000_000_000_000_000 # 10 Sextilliarden
|
||
@show typeof(z);
|
||
```
|
||
|
||
Meist wird man allerdings den Datentyp `BigInt` explizit anfordern müssen, damit nicht modulo $2^{64}$ gerechnet wird:
|
||
|
||
```{julia}
|
||
@show 3^300 BigInt(3)^300;
|
||
```
|
||
|
||
*Arbitrary precision arithmetic* kostet einiges an Speicherplatz und Rechenzeit.
|
||
|
||
Wir vergleichen den Zeit- und Speicherbedarf bei der Aufsummation von 10 Millionen Ganzzahlen als `Int64` und als `BigInt`.
|
||
|
||
```{julia}
|
||
# 10^7 Zufallszahlen, gleichverteilt zwischen -10^7 und 10^7
|
||
vec_int = rand(-10^7:10^7, 10^7)
|
||
|
||
# Dieselben Zahlen als BigInts
|
||
vec_bigint = BigInt.(vec_int)
|
||
```
|
||
|
||
Einen ersten Eindruck vom Zeit- und Speicherbedarf gibt das `@time`-Macro:
|
||
|
||
```{julia}
|
||
@time x = sum(vec_int)
|
||
@show x typeof(x)
|
||
```
|
||
```{julia}
|
||
@time x = sum(vec_bigint)
|
||
@show x typeof(x);
|
||
```
|
||
|
||
Durch die Just-in-Time-Compilation von Julia ist die einmalige Ausführung einer Funktion wenig aussagekräftig. Das Paket `BenchmarkTools` stellt u.a. das Macro `@benchmark` bereit, das eine Funktion mehrfach aufruft und die Ausführungszeiten als Histogramm darstellt.
|
||
|
||
:::{.ansitight}
|
||
```{julia}
|
||
using BenchmarkTools
|
||
|
||
@benchmark sum($vec_int)
|
||
```
|
||
|
||
|
||
```{julia}
|
||
@benchmark sum($vec_bigint)
|
||
```
|
||
:::
|
||
Die `BigInt`-Addition ist mehr als 30 mal langsamer.
|
||
|
||
|
||
|
||
:::{.content-hidden unless-format="xxx"}
|
||
Die folgende Funktion soll die Summe aller Zahlen von 1 bis n mit der Arithmetik des Datentyps T berechnen.
|
||
Auf Grund der *type promotion rules* reicht es für `T ≥ Int64` dazu aus, die Akkumulatorvariable mit einer Zahl vom Typ T zu initialisieren.
|
||
```{julia}
|
||
function mysum(n, T)
|
||
s = T(0)
|
||
for i = 1:n
|
||
s += i
|
||
end
|
||
return s
|
||
end
|
||
```
|
||
|
||
Einen ersten Eindruck vom Zeit- und Speicherbedarf gibt das `@time`-Macro:
|
||
|
||
```{julia}
|
||
@time x = mysum(10_000_000, Int64)
|
||
@show x typeof(x);
|
||
```
|
||
|
||
```{julia}
|
||
@time x = mysum(10_000_000, BigInt)
|
||
@show x typeof(x);
|
||
```
|
||
|
||
Durch die Just-in-Time-Compilation von Julia ist die einmalige Ausführung einer Funktion wenig aussagekräftig. Das Paket `BenchmarkTools` stellt u.a. das Macro `@benchmark` bereit, das eine Funktion mehrfach aufruft und die Ausführungszeiten als Histogramm darstellt.
|
||
|
||
:::{.ansitight}
|
||
```{julia}
|
||
using BenchmarkTools
|
||
|
||
@benchmark mysum(10_000_000, Int64)
|
||
```
|
||
|
||
|
||
```{julia .myclass}
|
||
|
||
@benchmark mysum(10_000_000, BigInt)
|
||
```
|
||
|
||
|
||
Die Berechnung von $\sum_{n=1}^{10000000} n$ dauert auf meinem PC mit den Standard-64bit-Integern im Schnitt 2 Nanosekunden, in *arbitrary precision arithmetic* über eine Sekunde, wobei auch noch fast 500MB Speicher allokiert werden.
|
||
|
||
:::
|
||
:::
|
||
|
||
|
||
## Gleitkommazahlen
|
||
|
||
Aus _floating point numbers_ kann man im Deutschen **[Gleit|Fließ]--[Komma|Punkt]--Zahlen** machen
|
||
und tatsächlich kommen alle 4 Varianten in der Literatur vor.
|
||
|
||
In der Numerischen Mathematik spricht man auch gerne von **Maschinenzahlen**.
|
||
|
||
|
||
### Grundidee
|
||
|
||
- Eine „feste Anzahl von Vor- und Nachkommastellen“ ist für viele Probleme ungeeignet.
|
||
- Eine Trennung zwischen „gültigen Ziffern“ und Größenordnung wie in der wissenschaftlichen Notation
|
||
$$ 345.2467 \times 10^3\qquad 34.52467\times 10^4\qquad \underline{3.452467\times 10^5}\qquad 0.3452467\times 10^6\qquad 0.03452467\times 10^7$$
|
||
ist wesentlich flexibler.
|
||
- Zur Eindeutigkeit legt man fest: Die Darstellung mit genau einer Ziffer vor dem Komma ist die __Normalisierte Darstellung__.
|
||
- Bei Binärzahlen `1.01101`: ist diese Ziffer immer gleich Eins und man kann auf das Abspeichern dieser Ziffer verzichten.
|
||
|
||
|
||
:::{.callout-note }
|
||
|
||
## Maschinenzahlen
|
||
|
||
Die Menge der Maschinenzahlen $𝕄(b, m, e_{min}, e_{max})$ ist charakterisiert durch die verwendete Basis $b$, die Mantissenlänge $m$ und den Wertebereich des Exponenten $\{e_{min}, ... ,e_{max}\}$.
|
||
|
||
In unserer Konvention hat die Mantisse einer normalisierten Maschinenzahl eine Ziffer (der Basis $b$) ungleich Null vor dem Komma und $m-1$ Nachkommastellen.
|
||
|
||
Wenn $b=2$ ist, braucht man daher nur $m-1$ Bits zur Speicherung der Mantisse normalisierter Gleitkommazahlen.
|
||
|
||
Der Standard IEEE 754, der von der Mahrzahl der modernen Prozessoren und Programmiersprachen implementiert wird, definiert
|
||
|
||
- `Float32` als $𝕄(2, 24, -126, 127 )$ und
|
||
- `Float64` als $𝕄(2, 53, -1022, 1023 ).$
|
||
|
||
:::
|
||
|
||
### Aufbau von `Float64` nach [Standard IEEE 754](https://de.wikipedia.org/wiki/IEEE_754)
|
||
|
||
|
||
::: {.content-visible when-format="html"}
|
||
![Aufbau einer `Float64` (Quelle:^[Quelle: <a href="https://commons.wikimedia.org/wiki/File:IEEE_754_Double_Floating_Point_Format.svg">Codekaizen</a>, <a href="https://creativecommons.org/licenses/by-sa/4.0">CC BY-SA 4.0</a>, via Wikimedia Commons ])
|
||
](../images/1024px-IEEE_754_Double_Floating_Point_Format.png)
|
||
:::
|
||
|
||
::: {.content-visible when-format="pdf"}
|
||
{width="70%"}
|
||
:::
|
||
|
||
|
||
- 1 Vorzeichenbit $S$
|
||
- 11 Bits für den Exponenten, also $0\le E \le 2047$
|
||
- die Werte $E=0$ und $E=(11111111111)_2=2047$ sind reserviert für die Codierung von Spezialwerten wie
|
||
$\pm0, \pm\infty$, NaN _(Not a Number)_ und denormalisierte Zahlen.
|
||
- 52 Bits für die Mantisse $M,\quad 0\le M<1$, das entspricht etwa 16 Dezimalstellen
|
||
- Damit wird folgende Zahl dargestellt:
|
||
$$ x=(-1)^S \cdot(1+M)\cdot 2^{E-1023}$$
|
||
|
||
Ein Beispiel:
|
||
```{julia}
|
||
x = 27.56640625
|
||
bitstring(x)
|
||
```
|
||
Das geht auch schöner:
|
||
|
||
```{julia}
|
||
function printbitsf64(x::Float64)
|
||
s = bitstring(x)
|
||
printstyled(s[1], color = :blue, reverse=true)
|
||
printstyled(s[2:12], color = :green, reverse=true)
|
||
printstyled(s[13:end], color=:red, bold=true, reverse=true)
|
||
print("\n")
|
||
end
|
||
|
||
printbitsf64(x)
|
||
```
|
||
und wir können S (in blau),E (grün) und M (rot) ablesen:
|
||
$$
|
||
\begin{aligned}
|
||
S &= 0\\
|
||
E &= (10000000011)_2 = 2^{10} + 2^1 + 2^0 = 1027\\
|
||
M &= (.1011100100001)_2 = \frac{1}{2} + \frac{1}{2^3} + \frac{1}{2^4} + \frac{1}{2^5} + \frac{1}{2^8} + \frac{1}{2^{12}}\\
|
||
x &=(-1)^S \cdot(1+M)\cdot 2^{E-1023}
|
||
\end{aligned}
|
||
$$
|
||
|
||
|
||
... und damit die Zahl rekonstruieren:
|
||
|
||
```{julia}
|
||
x = (1 + 1/2 + 1/8 + 1/16 + 1/32 + 1/256 + 1/4096) * 2^4
|
||
```
|
||
|
||
- Die Maschinenzahlen 𝕄 bilden eine endliche, diskrete Untermenge von ℝ. Es gibt eine kleinste und eine größte Maschinenzahl und abgesehen davon haben alle x∈𝕄 einen Vorgänger und Nachfolger in 𝕄.
|
||
- Was ist der Nachfolger von x in 𝕄? Dazu setzen wir das kleinste Mantissenbit von 0 auf 1.
|
||
- Die Umwandlung eines Strings aus Nullen und Einsen in die entsprechende Maschinenzahl ist z.B. so möglich:
|
||
|
||
|
||
```{julia}
|
||
ux = parse(UInt64, "0100000000111011100100010000000000000000000000000000000000000001", base=2)
|
||
reinterpret(Float64, ux)
|
||
```
|
||
|
||
Das kann man in Julia allerdings auch einfacher haben:
|
||
|
||
```{julia}
|
||
y = nextfloat(x)
|
||
```
|
||
|
||
Der Vorgänger von x ist:
|
||
|
||
```{julia}
|
||
z = prevfloat(x)
|
||
println(z)
|
||
printbitsf64(z)
|
||
```
|
||
|
||
|
||
## Maschinenepsilon
|
||
|
||
- Den Abstand zwischen `1` und dem Nachfolger `nextfloat(1)` nennt man [**Maschinenepsilon**](https://en.wikipedia.org/wiki/Machine_epsilon).
|
||
- Für `Float64` mit einer Mantissenlänge von 52 Bit, ist $\epsilon=2^{-52}$.
|
||
|
||
```{julia}
|
||
@show nextfloat(1.) - 1 2^-52 eps(Float64);
|
||
```
|
||
|
||
- Das Maschinenepsilon ist ein Maß für den relativen Abstand zwischen den Maschinenzahlen und quantifiziert die Aussage: „64-Bit-Gleitkommazahlen haben eine Genauigkeit von etwa 16 Dezimalstellen.“
|
||
- Das Maschinenepsilon ist etwas völlig anderes als die kleinste von Null verschiedene Gleitkommazahl:
|
||
|
||
```{julia}
|
||
floatmin(Float64)
|
||
```
|
||
|
||
- Ein Teil der Literatur verwendet eine andere Definition des Maschinenepsilons, die halb so groß ist.
|
||
$$
|
||
\epsilon' = \frac{\epsilon}{2}\approx 1.1\times 10^{-16}
|
||
$$
|
||
ist der maximale relative Fehler, der beim Runden einer reellen Zahl auf die nächste Maschinenzahl entstehen kann.
|
||
- Da Zahlen aus dem Intervall $(1-\epsilon',1+\epsilon']$ auf die Maschinenzahl $1$ gerundet werden, kann man $\epsilon'$ auch definieren als: *die größte Zahl, für die in der Maschinenzahlarithmetik noch gilt: $1+\epsilon' = 1$.*
|
||
|
||
Auf diese Weise kann man das Maschinenepsilon auch berechnen:
|
||
|
||
:::{.ansitight}
|
||
|
||
```{julia}
|
||
Eps = 1
|
||
while(1 != 1 + Eps)
|
||
Eps /= 2
|
||
println(1+Eps)
|
||
end
|
||
Eps
|
||
```
|
||
|
||
oder als Bitmuster:
|
||
|
||
```{julia}
|
||
Eps=1
|
||
while 1 != 1 + Eps
|
||
Eps/= 2
|
||
printbitsf64(1+Eps)
|
||
end
|
||
Eps
|
||
```
|
||
|
||
:::
|
||
|
||
Die größte und die kleinste positive normalisiert darstellbare Gleitkommazahl eines Gleitkommatyps kann man abfragen:
|
||
|
||
```{julia}
|
||
@show floatmax(Float64)
|
||
printbitsf64(floatmax(Float64))
|
||
|
||
@show floatmin(Float64)
|
||
printbitsf64(floatmin(Float64))
|
||
```
|
||
|
||
|
||
|
||
|
||
|
||
## Runden auf Maschinenzahlen
|
||
|
||
- Die Abbildung rd: ℝ $\rightarrow$ 𝕄 soll zur nächstgelegenen darstellbaren Zahl runden.
|
||
- Standardrundungsregel: _round to nearest, ties to even_
|
||
Wenn man genau die Mitte zwischen zwei Maschinenzahlen trifft *(tie)*, wählt man die, deren letztes Mantissenbit 0 ist.
|
||
- Begründung: damit wird statistisch in 50% der Fälle auf- und in 50% der Fälle abgerundet und so ein „statistischer Drift“ bei längeren Rechnungen vermieden.
|
||
|
||
## Maschinenzahlarithmetik
|
||
|
||
Die Maschinenzahlen sind als Untermenge von ℝ nicht algebraisch abgeschlossen. Schon die Summe zweier Maschinenzahlen wird in der Regel keine Maschinenzahl sein.
|
||
|
||
:::{.callout-important}
|
||
Der Standard IEEE 754 fordert, dass die Maschinenzahlarithmetik das *gerundete exakte Ergebnis* liefert:
|
||
|
||
Das Resultat muss gleich demjenigen sein, das bei einer exakten Ausführung der entsprechenden Operation mit anschließender Rundung entsteht.
|
||
$$
|
||
a \oplus b = \text{rd}(a + b)
|
||
$$
|
||
Analoges muss für die Implemetierung der Standardfunktionen wie
|
||
wie `sqrt()`, `log()`, `sin()` ... gelten: Sie liefern ebenfalls die Maschinenzahl, die dem exakten Ergebnis am nächsten kommt.
|
||
:::
|
||
|
||
|
||
Die Arithmetik ist *nicht assoziativ*:
|
||
|
||
```{julia}
|
||
1 + 10^-16 + 10^-16
|
||
```
|
||
|
||
```{julia}
|
||
1 + (10^-16 + 10^-16)
|
||
```
|
||
|
||
Im ersten Fall (ohne Klammern) wird von links nach rechts ausgewertet:
|
||
$$
|
||
\begin{aligned}
|
||
1 \oplus 10^{-16} \oplus 10^{-16} &=
|
||
(1 \oplus 10^{-16}) \oplus 10^{-16} \\ &=
|
||
\text{rd}\left(\text{rd}(1 + 10^{-16}) + 10^{-16} \right)\\
|
||
&= \text{rd}(1 + 10^{-16})\\
|
||
&= 1
|
||
\end{aligned}
|
||
$$
|
||
|
||
Im zweiten Fall erhält man:
|
||
$$
|
||
\begin{aligned}
|
||
1 \oplus \left(10^{-16} \oplus 10^{-16}\right) &=
|
||
\text{rd}\left(1 + \text{rd}(10^{-16} + 10^{-16}) \right)\\
|
||
&= \text{rd}(1 + 2*10^{-16})\\
|
||
&= 1.0000000000000002
|
||
\end{aligned}
|
||
$$
|
||
|
||
|
||
Es sei auch daran erinnert, dass sich selbst „einfache“ Dezimalbrüche häufig nicht exakt als Maschinenzahlen darstellen lassen:
|
||
|
||
$$
|
||
\begin{aligned}
|
||
(0.1)_{10} &= (0.000110011001100110011001100...)_2 = (0.000\overline{1100})_2 \\
|
||
(0.3)_{10} &= (0.0100110011001100110011001100..)_2 = (0.0\overline{1001})_2
|
||
\end{aligned}
|
||
$$
|
||
|
||
|
||
```{julia}
|
||
printbitsf64(0.1)
|
||
printbitsf64(0.3)
|
||
```
|
||
|
||
Folge:
|
||
|
||
```{julia}
|
||
0.1 + 0.1 == 0.2
|
||
```
|
||
|
||
```{julia}
|
||
0.2 + 0.1 == 0.3
|
||
```
|
||
|
||
```{julia}
|
||
0.2 + 0.1
|
||
```
|
||
|
||
Bei der Ausgabe einer Maschinenzahl muss der Binärbruch in einen Dezimalbruch entwickelt werden. Man kann sich auch mehr Stellen dieser Dezimalbruchentwicklung anzeigen lassen:
|
||
```{julia}
|
||
using Printf
|
||
@printf("%.30f", 0.1)
|
||
```
|
||
|
||
```{julia}
|
||
@printf("%.30f", 0.3)
|
||
```
|
||
Die Binärbruch-Mantisse einer Maschinenzahl kann eine lange oder sogar unendlich-periodische Dezimalbruchentwicklung haben. Dadurch
|
||
sollte man sich nicht eine „höheren Genauigkeit“ suggerieren lassen!
|
||
|
||
:::{.callout-important}
|
||
Moral: wenn man `Float`s auf Gleichheit testen will, sollte man fast immer eine dem Problem angemessene realistische Genauigkeit `epsilon` festlegen und
|
||
darauf testen:
|
||
|
||
```julia
|
||
epsilon = 1.e-10
|
||
|
||
if abs(x-y) < epsilon
|
||
# ...
|
||
end
|
||
```
|
||
:::
|
||
|
||
## Normalisierte und Denormalisierte Maschinenzahlen
|
||
|
||
Zum Verständnis nehmen wir ein einfaches Modell:
|
||
|
||
- Sei 𝕄(10,4,±5) die Menge der Maschinenzahlen zur Basis 10 mit 4 Mantissenstellen (eine vor dem Komma, 3 Nachkommastellen) und dem Exponentenbereich -5 ≤ E ≤ 5.
|
||
- Dann ist die normalisierte Darstellung (Vorkommastelle ungleich 0)
|
||
von z.B. 1234.0 gleich 1.234e3 und von 0.00789 gleich 7.890e-3.
|
||
- Es ist wichtig, dass die Maschinenzahlen bei jedem Rechenschritt normalisiert gehalten werden. Nur so wird die Mantissenlänge voll ausgenutzt und die Genauigkeit ist maximal.
|
||
- Die kleinste positive normalisierte Zahl in unserem Modell ist `x = 1.000e-5`. Schon `x/2` müsste auf 0 gerundet werden.
|
||
- Hier erweist es sich für viele Anwendungen als günstiger, auch denormalisierte *(subnormal)* Zahlen zuzulassen und `x/2` als `0.500e-5` oder `x/20` als `0.050e-5` darzustellen.
|
||
- Dieser *gradual underflow* ist natürlich mit einem Verlust an gültigen Stellen und damit Genauigkeit verbunden.
|
||
|
||
|
||
Im `Float`-Datentyp werden solche *subnormal values* dargestellt durch ein Exponentenfeld, in dem alle Bits gleich Null sind:
|
||
|
||
```{julia}
|
||
#| echo: false
|
||
#| output: false
|
||
flush(stdout)
|
||
```
|
||
|
||
|
||
:::{.ansitight}
|
||
```{julia}
|
||
using Printf
|
||
|
||
x = 2 * floatmin(Float64) # 2*kleinste normalisierte Gleitkommazahl > 0
|
||
|
||
while x != 0
|
||
x /= 2
|
||
@printf "%14.6g " x
|
||
printbitsf64(x)
|
||
end
|
||
```
|
||
:::
|
||
|
||
## Spezielle Werte
|
||
|
||
Die Gleitkommaarithmetik kennt einige spezielle Werte:
|
||
|
||
```{julia}
|
||
for x ∈ (NaN, Inf, -Inf, -0.0)
|
||
@printf "%6g " x
|
||
printbitsf64(x)
|
||
end
|
||
```
|
||
|
||
- Ein Exponentenüberlauf *(overflow)* führt zum Ergebnis `Inf` oder `-Inf`.
|
||
```{julia}
|
||
2/0, -3/0, floatmax(Float64) * 1.01, exp(1300)
|
||
```
|
||
- Damit kann weitergerechnet werden:
|
||
|
||
```{julia}
|
||
-Inf + 20, Inf/30, 23/-Inf, sqrt(Inf), Inf * 0, Inf - Inf
|
||
```
|
||
- `NaN` *(Not a Number)* steht für das Resultat einer Operation, das undefiniert ist. Alle weiteren Operationen mit `NaN` ergeben ebenfalls `NaN`.
|
||
|
||
```{julia}
|
||
0/0, Inf - Inf, 2.3NaN, sqrt(NaN)
|
||
```
|
||
|
||
- Da `NaN` einen undefinierten Wert repräsentiert, ist es zu nichts gleich, nichtmal zu sich selbst. Das ist sinnvoll, denn wenn zwei Variablen `x` und `y` als `NaN` berechnet wurden, sollte man nicht schlußfolgern, dass sie gleich sind.
|
||
- Zum Testen auf `NaN` gibt es daher die boolsche Funktion `isnan()`.
|
||
|
||
```{julia}
|
||
x = 0/0
|
||
y = Inf - Inf
|
||
@show x==y NaN==NaN isfinite(NaN) isinf(NaN) isnan(x) isnan(y);
|
||
```
|
||
|
||
- Es gibt eine „minus Null“. Sie signalisiert einen Exponentenunterlauf *(underflow)* einer betragsmäßig zu klein gewordenen *negativen* Größe.
|
||
|
||
```{julia}
|
||
@show 23/-Inf -2/exp(1200) -0.0==0.0;
|
||
```
|
||
|
||
## Mathematische Funktionen
|
||
|
||
Julia verfügt über die [üblichen mathematischen Funktionen](https://docs.julialang.org/en/v1/manual/mathematical-operations/#Rounding-functions)
|
||
`sqrt, exp, log, log2, log10, sin, cos,..., asin, acos,..., sinh,..., gcd, lcm, factorial,...,abs, max, min,...`,
|
||
|
||
darunter z.B. die [Rundungsfunktionen](https://de.wikipedia.org/wiki/Abrundungsfunktion_und_Aufrundungsfunktion)
|
||
|
||
- `floor(T,x)` = $\lfloor x \rfloor$
|
||
- `ceil(T,x)` = $\lceil x \rceil$
|
||
|
||
```{julia}
|
||
floor(3.4), floor(Int64, 3.5), floor(Int64, -3.5)
|
||
```
|
||
|
||
```{julia}
|
||
ceil(3.4), ceil(Int64, 3.5), ceil(Int64, -3.5)
|
||
```
|
||
|
||
Es sei noch hingewiesen auf `atan(y, x)`, den [Arkustangens mit 2 Argumenten](https://de.wikipedia.org/wiki/Arctan2), Er ist in anderen Programmiersprachen oft als Funktion mit eigenem Namen *atan2* implementiert.
|
||
Dieser löst das Problem der Umrechnung von kartesischen in Polarkoordinaten ohne lästige Fallunterscheidung.
|
||
|
||
- `atan(y,x)` ist Winkel der Polarkoordinaten von (x,y) im Intervall $(-\pi,\pi]$. Im 1. und 4. Quadranten ist er also gleich `atan(y/x)`
|
||
|
||
```{julia}
|
||
atan(3, -2), atan(-3, 2), atan(-3/2)
|
||
```
|
||
|
||
|
||
## Umwandlung Strings $\Longleftrightarrow$ Zahlen
|
||
|
||
Die Umwandlung ist mit den Funktionen `parse()` und `string()` möglich.
|
||
|
||
```{julia}
|
||
parse(Int64, "1101", base=2)
|
||
```
|
||
|
||
```{julia}
|
||
string(13, base=2)
|
||
```
|
||
|
||
```{julia}
|
||
string(1/7)
|
||
```
|
||
|
||
```{julia}
|
||
string(77, base=16)
|
||
```
|
||
|
||
Zur Umwandlung der numerischen Typen ineinander kann man die Typnamen verwenden. Typenamen sind auch Konstruktoren:
|
||
|
||
|
||
```{julia}
|
||
x = Int8(67)
|
||
@show x typeof(x);
|
||
```
|
||
|
||
```{julia}
|
||
z = UInt64(3459)
|
||
```
|
||
|
||
|
||
```{julia}
|
||
y = Float64(z)
|
||
```
|
||
|
||
|
||
|
||
## Literatur
|
||
|
||
- D. Goldberg, [What Every Computer Scientist Should Know About Floating-Point Arithmetic](http://www.validlab.com/goldberg/paper.pdf)
|
||
- C. Vuik, [Some Disasters caused by numerical errors](http://ta.twi.tudelft.nl/users/vuik/wi211/disasters.html)
|