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
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.
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.
- Zur Eindeutigkeit muss man eine dieser Formen auswählen. In der mathematischen Analyse von Maschinenzahlen wählt man oft die Form, bei der die erste Nachkommastelle ungleich Null ist. Damit gilt für die Mantisse $m$:
- Bei Binärzahlen `1.01101`: ist diese Ziffer immer gleich Eins und man kann auf das Abspeichern dieser Ziffer verzichten. Diese tatsächlich abgespeicherte (gekürzte) Mantisse bezeichnen wir mit $M$, so dass
Die Menge der Maschinenzahlen $𝕄(b, p, e_{min}, e_{max})$ ist charakterisiert durch die verwendete Basis $b$, die Mantissenlänge $p$ 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 $p-1$ Nachkommastellen.
- 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:
- 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.“
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:
- Im Intervall $[1,2)$ liegen $2^{52}$ äquidistante Maschinenzahlen.
- Danach erhöht sich der Exponent um 1 und die Mantisse $M$ wird auf 0 zurückgesetzt. Damit enthält das Intervall $[2,4)$ wiederum $2^{52}$ äquidistante Maschinenzahlen, ebenso das Intervall $[4,8)$ bis hin zu $[2^{1023}, 2^{1024})$.
- Ebenso liegen in den Intervallen $\ [\frac{1}{2},1), \ [\frac{1}{4},\frac{1}{2}),...$ je $2^{52}$ äquidistante Maschinenzahlen, bis hinunter zu $[2^{-1022}, 2^{-1021})$.
- Dies bildet die Menge $𝕄_+$ der positiven Maschinenzahlen und es ist
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.
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:
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
- 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:
- `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()`.
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)`
