6  Maschinenzahlen

for x ∈ ( 3, 3.3e4,  Int16(20), Float32(3.3e4), UInt16(9) )
    @show x sizeof(x) typeof(x)
    println()
end
x = 3
sizeof(x) = 8
typeof(x) = Int64

x = 33000.0
sizeof(x) = 8
typeof(x) = Float64

x = 20
sizeof(x) = 2
typeof(x) = Int16

x = 33000.0f0
sizeof(x) = 4
typeof(x) = Float32

x = 0x0009
sizeof(x) = 2
typeof(x) = UInt16

6.1 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.

6.1.1 Unsigned integers

subtypes(Unsigned)
5-element Vector{Any}:
 UInt128
 UInt16
 UInt32
 UInt64
 UInt8

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.

x = 0x0033efef
@show x typeof(x) Int(x)

z = UInt(32)
@show z typeof(z);
x = 0x0033efef
typeof(x) = UInt32
Int(x) = 3403759
z = 0x0000000000000020
typeof(z) = UInt64

6.1.2 Signed Integers

subtypes(Signed)
6-element Vector{Any}:
 BigInt
 Int128
 Int16
 Int32
 Int64
 Int8

Integers haben den Wertebereich \[ -2^{n-1} \le x < 2^{n-1} \]

In Julia sind ganze Zahlen standardmäßig 64 Bit groß:

x = 42
typeof(x)
Int64

Sie haben daher den Wertebereich: \[ -9.223.372.036.854.775.808 \le x \le 9.223.372.036.854.775.807 \]

32-Bit-Integers haben den Wertebereich \[ -2.147.483.648 \le x \le 2.147.483.647 \] Der Maximalwert \(2^{31}-1\) is zufällig gerade eine Mersenne-Primzahl:

using Primes 
isprime(2^31-1)
true

Negative Zahlen werden im Zweierkomplement dargestellt:

\(x \Rightarrow -x\) entspricht: flip all bits, then add 1

Das sieht so aus:

Eine Darstellung des fiktiven Datentyps Int4

Das höchstwertige Bit zeigt das Vorzeichen an. Sein „Umkippen“ entspricht allerdings nicht der Operation x → -x.

Wichtig

Die gesamte Ganzzahlarithmetik ist eine Arithmetik modulo \(2^n\)

x = 2^62 - 10 + 2^62
9223372036854775798
x + 20
-9223372036854775798

Keine Fehlermeldung, keine Warnung! Ganzzahlen fester Länge liegen nicht auf einer Geraden, sondern auf einem Kreis!

Schauen wir uns ein paar Integers an:

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
                   1 Int64 ⇒  0000000000000000000000000000000000000000000000000000000000000001
                   2 Int64 ⇒  0000000000000000000000000000000000000000000000000000000000000010
                   7 Int64 ⇒  0000000000000000000000000000000000000000000000000000000000000111
                  -7 Int64 ⇒  1111111111111111111111111111111111111111111111111111111111111001
    1125899906842624 Int64 ⇒  0000000000000100000000000000000000000000000000000000000000000000
                  -1 Int64 ⇒  1111111111111111111111111111111111111111111111111111111111111111
                   7 Int16 ⇒  0000000000000111
                   7 Int8  ⇒  00000111

Julia hat Funktionen, die über die Datentypen informieren (introspection):

typemin(Int64), typemax(Int64)
(-9223372036854775808, 9223372036854775807)
typemin(UInt64), typemax(UInt64), BigInt(typemax(UInt64))
(0x0000000000000000, 0xffffffffffffffff, 18446744073709551615)
typemin(Int8), typemax(Int8)
(-128, 127)

6.2 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 gleich 1.
(-2)^63, 2^64, 3^(-3), 0^0
(-9223372036854775808, 0, 0.037037037037037035, 1)
  • Für natürliche Exponenten wird exponentiation by squaring 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.
x = 40/5
8.0

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:
@show divrem( 27, 4)
@show ( 27 ÷  4,  27 %  4)   
@show (-27 ÷  4, -27 %  4 )
@show ( 27 ÷ -4,  27 % -4);
divrem(27, 4) = (6, 3)
(27 ÷ 4, 27 % 4) = (6, 3)
(-27 ÷ 4, -27 % 4) = (-6, -3)
(27 ÷ -4, 27 % -4) = (-6, 3)
  • 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):
@show divrem(-27, 4, RoundDown)
@show (fld(-27,  4), mod(-27,  4))
@show (fld( 27, -4), mod( 27, -4));
divrem(-27, 4, RoundDown) = (-7, 1)
(fld(-27, 4), mod(-27, 4)) = (-7, 1)
(fld(27, -4), mod(27, -4)) = (-7, -1)

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:

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);
typeof(z) = Int64
typeof(z) = Int64
typeof(z) = Int128
typeof(z) = BigInt

Meist wird man allerdings den Datentyp BigInt explizit anfordern müssen, damit nicht modulo \(2^{64}\) gerechnet wird:

@show 3^300        BigInt(3)^300;
3 ^ 300 = 4157753088978724465
BigInt(3) ^ 300 = 136891479058588375991326027382088315966463695625337436471480190078368997177499076593800206155688941388250484440597994042813512732765695774566001

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.

# 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)
10000000-element Vector{BigInt}:
  1254509
  2789949
  -179922
 -6472591
 -9800554
 -4320501
   738617
 -9238805
  9482724
 -7719974
        ⋮
 -2153395
 -7458756
  7107169
  6789723
  9531434
 -4123834
  9023344
  3436240
    39559

Einen ersten Eindruck vom Zeit- und Speicherbedarf gibt das @time-Macro:

@time x = sum(vec_int)
@show x typeof(x)
  0.007893 seconds (491 allocations: 22.984 KiB, 64.47% compilation time)
x = -32235833853
typeof(x) = Int64
Int64
@time x = sum(vec_bigint)
@show x typeof(x);
  0.091471 seconds (53.71 k allocations: 2.653 MiB, 11.84% compilation time)
x = -32235833853
typeof(x) = BigInt

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.

using BenchmarkTools

@benchmark sum($vec_int)
BenchmarkTools.Trial: 1953 samples with 1 evaluation per sample.
 Range (min … max):  2.422 ms …  3.381 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.524 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.550 ms ± 87.910 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

        ▃▅▆▇▆███▄▄▁▂▁                                         
  ▂▂▄▅▆▇█████████████▇██▆▇▇▆▆▆▆▅▆▅▄▃▄▅▄▄▃▄▄▃▃▃▃▃▂▃▃▃▂▂▂▁▂▁▂▂ ▄
  2.42 ms        Histogram: frequency by time        2.83 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.
@benchmark sum($vec_bigint)
BenchmarkTools.Trial: 61 samples with 1 evaluation per sample.
 Range (min … max):  80.020 ms … 96.473 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     81.507 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   82.641 ms ±  3.988 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

    █▅ ▃▂▃                                                     
  ▅▇██▃███▃▇▅▁▁▁▃▁▁▁▁▁▁▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▁▃▃▁▁▃ ▁
  80 ms           Histogram: frequency by time          96 ms <

 Memory estimate: 48 bytes, allocs estimate: 3.

Die BigInt-Addition ist mehr als 30 mal langsamer.

6.3 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.

6.3.1 Grundidee

  • Eine „feste Anzahl von Vor- und Nachkommastellen“ ist für viele Probleme ungeeignet.
  • Eine Trennung zwischen „gültigen Ziffern“ und Größenordnung (Mantisse und Exponent) wie in der wissenschaftlichen Notation ist wesentlich flexibler. \[ 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\]
  • 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\): \[ 1 > m \ge (0.1)_b = b^{-1}, \]
    wobei \(b\) die gewählte Basis des Zahlensystems bezeichnet.
  • Wir wählen im Folgenden die Form, die der tatsächlichen Implementation auf dem Computer entspricht und legen fest: Die Darstellung mit genau einer Ziffer ungleich Null vor dem Komma ist die Normalisierte Darstellung. Damit gilt \[ (10)_b = b> m \ge 1. \]
  • 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 \[m = 1 + M\] gilt.
Maschinenzahlen

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.

Wenn \(b=2\) ist, braucht man daher nur \(p-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 ).\)

6.3.2 Aufbau von Float64 nach Standard IEEE 754

Aufbau einer Float64 (Quelle:1)
  • 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 (gekürzte) 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:

x = 27.56640625
bitstring(x)
"0100000000111011100100010000000000000000000000000000000000000000"

Das geht auch schöner:

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

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:

x = (1 + 1/2 + 1/8 + 1/16 + 1/32 + 1/256 + 1/4096) * 2^4
27.56640625
  • 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:
ux = parse(UInt64, "0100000000111011100100010000000000000000000000000000000000000001", base=2)
reinterpret(Float64, ux)
27.566406250000004

Das kann man in Julia allerdings auch einfacher haben:

y = nextfloat(x)
27.566406250000004

Der Vorgänger von x ist:

z = prevfloat(x)
println(z)
printbitsf64(z)
27.566406249999996
0100000000111011100100001111111111111111111111111111111111111111

6.4 Maschinenepsilon

  • Den Abstand zwischen 1 und dem Nachfolger nextfloat(1) nennt man Maschinenepsilon.
  • Für Float64 mit einer Mantissenlänge von 52 Bit ist \(\epsilon=2^{-52}\).
@show nextfloat(1.) - 1   2^-52   eps(Float64);
nextfloat(1.0) - 1 = 2.220446049250313e-16
2 ^ -52 = 2.220446049250313e-16
eps(Float64) = 2.220446049250313e-16
  • 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 positive Gleitkommazahl:
floatmin(Float64)
2.2250738585072014e-308
  • 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:

Eps = 1
while(1 != 1 + Eps)
    Eps /= 2
    println(1+Eps)
end
Eps
1.5
1.25
1.125
1.0625
1.03125
1.015625
1.0078125
1.00390625
1.001953125
1.0009765625
1.00048828125
1.000244140625
1.0001220703125
1.00006103515625
1.000030517578125
1.0000152587890625
1.0000076293945312
1.0000038146972656
1.0000019073486328
1.0000009536743164
1.0000004768371582
1.000000238418579
1.0000001192092896
1.0000000596046448
1.0000000298023224
1.0000000149011612
1.0000000074505806
1.0000000037252903
1.0000000018626451
1.0000000009313226
1.0000000004656613
1.0000000002328306
1.0000000001164153
1.0000000000582077
1.0000000000291038
1.000000000014552
1.000000000007276
1.000000000003638
1.000000000001819
1.0000000000009095
1.0000000000004547
1.0000000000002274
1.0000000000001137
1.0000000000000568
1.0000000000000284
1.0000000000000142
1.000000000000007
1.0000000000000036
1.0000000000000018
1.0000000000000009
1.0000000000000004
1.0000000000000002
1.0
1.1102230246251565e-16

oder als Bitmuster:

Eps=1
while 1 != 1 + Eps
    Eps/= 2
    printbitsf64(1+Eps)
end
Eps
0011111111111000000000000000000000000000000000000000000000000000
0011111111110100000000000000000000000000000000000000000000000000
0011111111110010000000000000000000000000000000000000000000000000
0011111111110001000000000000000000000000000000000000000000000000
0011111111110000100000000000000000000000000000000000000000000000
0011111111110000010000000000000000000000000000000000000000000000
0011111111110000001000000000000000000000000000000000000000000000
0011111111110000000100000000000000000000000000000000000000000000
0011111111110000000010000000000000000000000000000000000000000000
0011111111110000000001000000000000000000000000000000000000000000
0011111111110000000000100000000000000000000000000000000000000000
0011111111110000000000010000000000000000000000000000000000000000
0011111111110000000000001000000000000000000000000000000000000000
0011111111110000000000000100000000000000000000000000000000000000
0011111111110000000000000010000000000000000000000000000000000000
0011111111110000000000000001000000000000000000000000000000000000
0011111111110000000000000000100000000000000000000000000000000000
0011111111110000000000000000010000000000000000000000000000000000
0011111111110000000000000000001000000000000000000000000000000000
0011111111110000000000000000000100000000000000000000000000000000
0011111111110000000000000000000010000000000000000000000000000000
0011111111110000000000000000000001000000000000000000000000000000
0011111111110000000000000000000000100000000000000000000000000000
0011111111110000000000000000000000010000000000000000000000000000
0011111111110000000000000000000000001000000000000000000000000000
0011111111110000000000000000000000000100000000000000000000000000
0011111111110000000000000000000000000010000000000000000000000000
0011111111110000000000000000000000000001000000000000000000000000
0011111111110000000000000000000000000000100000000000000000000000
0011111111110000000000000000000000000000010000000000000000000000
0011111111110000000000000000000000000000001000000000000000000000
0011111111110000000000000000000000000000000100000000000000000000
0011111111110000000000000000000000000000000010000000000000000000
0011111111110000000000000000000000000000000001000000000000000000
0011111111110000000000000000000000000000000000100000000000000000
0011111111110000000000000000000000000000000000010000000000000000
0011111111110000000000000000000000000000000000001000000000000000
0011111111110000000000000000000000000000000000000100000000000000
0011111111110000000000000000000000000000000000000010000000000000
0011111111110000000000000000000000000000000000000001000000000000
0011111111110000000000000000000000000000000000000000100000000000
0011111111110000000000000000000000000000000000000000010000000000
0011111111110000000000000000000000000000000000000000001000000000
0011111111110000000000000000000000000000000000000000000100000000
0011111111110000000000000000000000000000000000000000000010000000
0011111111110000000000000000000000000000000000000000000001000000
0011111111110000000000000000000000000000000000000000000000100000
0011111111110000000000000000000000000000000000000000000000010000
0011111111110000000000000000000000000000000000000000000000001000
0011111111110000000000000000000000000000000000000000000000000100
0011111111110000000000000000000000000000000000000000000000000010
0011111111110000000000000000000000000000000000000000000000000001
0011111111110000000000000000000000000000000000000000000000000000
1.1102230246251565e-16
Die Menge der (normalisierten) Maschinenzahlen
  • 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 \[ 𝕄 = -𝕄_+ \cup \{0\} \cup 𝕄_+ \]

Die größte und die kleinste positive normalisiert darstellbare Gleitkommazahl eines Gleitkommatyps kann man abfragen:

@show floatmax(Float64)
printbitsf64(floatmax(Float64))

@show floatmin(Float64)
printbitsf64(floatmin(Float64))
floatmax(Float64) = 1.7976931348623157e308
0111111111101111111111111111111111111111111111111111111111111111
floatmin(Float64) = 2.2250738585072014e-308
0000000000010000000000000000000000000000000000000000000000000000

6.5 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.
  • Es gilt: \[ \frac{|x-\text{rd}(x)|}{|x|} \le \frac{1}{2} \epsilon \]

6.6 Maschinenzahlarithmetik

Die Maschinenzahlen sind als Untermenge von ℝ nicht algebraisch abgeschlossen. Schon die Summe zweier Maschinenzahlen wird in der Regel keine Maschinenzahl sein.

Wichtig

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:

1 +  10^-16 + 10^-16
1.0
1 + (10^-16 + 10^-16)
1.0000000000000002

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} \]

printbitsf64(0.1)
printbitsf64(0.3)
0011111110111001100110011001100110011001100110011001100110011010
0011111111010011001100110011001100110011001100110011001100110011

Folge:

0.1 + 0.1 == 0.2
true
0.2 + 0.1 == 0.3
false
0.2 + 0.1
0.30000000000000004

Bei der Ausgabe einer Maschinenzahl muss der Binärbruch in einen Dezimalbruch entwickelt werden. Man kann sich auch mehr Stellen dieser Dezimalbruchentwicklung anzeigen lassen:

using Printf
@printf("%.30f", 0.1)
0.100000000000000005551115123126
@printf("%.30f", 0.3)
0.299999999999999988897769753748

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!

Wichtig

Moral: wenn man Floats auf Gleichheit testen will, sollte man fast immer eine dem Problem angemessene realistische Genauigkeit epsilon festlegen und darauf testen:

epsilon = 1.e-10

if abs(x-y) < epsilon
   # ...
end

6.7 Normalisierte und Denormalisierte Maschinenzahlen

Die Lücke zwischen Null und der kleinsten normalisierten Maschinenzahl \(2^{-1022} \approx 2.22\times 10^{-308}\) ist mit denormalisierten Maschinenzahlen besiedelt.

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:

using Printf

x = 2 * floatmin(Float64)   # 2*kleinste normalisierte Gleitkommazahl > 0

while x != 0
    x /= 2
    @printf "%14.6g  " x
    printbitsf64(x)
end
  2.22507e-308  0000000000010000000000000000000000000000000000000000000000000000
  1.11254e-308  0000000000001000000000000000000000000000000000000000000000000000
  5.56268e-309  0000000000000100000000000000000000000000000000000000000000000000
  2.78134e-309  0000000000000010000000000000000000000000000000000000000000000000
  1.39067e-309  0000000000000001000000000000000000000000000000000000000000000000
  6.95336e-310  0000000000000000100000000000000000000000000000000000000000000000
  3.47668e-310  0000000000000000010000000000000000000000000000000000000000000000
  1.73834e-310  0000000000000000001000000000000000000000000000000000000000000000
  8.69169e-311  0000000000000000000100000000000000000000000000000000000000000000
  4.34585e-311  0000000000000000000010000000000000000000000000000000000000000000
  2.17292e-311  0000000000000000000001000000000000000000000000000000000000000000
  1.08646e-311  0000000000000000000000100000000000000000000000000000000000000000
  5.43231e-312  0000000000000000000000010000000000000000000000000000000000000000
  2.71615e-312  0000000000000000000000001000000000000000000000000000000000000000
  1.35808e-312  0000000000000000000000000100000000000000000000000000000000000000
  6.79039e-313  0000000000000000000000000010000000000000000000000000000000000000
  3.39519e-313  0000000000000000000000000001000000000000000000000000000000000000
   1.6976e-313  0000000000000000000000000000100000000000000000000000000000000000
  8.48798e-314  0000000000000000000000000000010000000000000000000000000000000000
  4.24399e-314  0000000000000000000000000000001000000000000000000000000000000000
    2.122e-314  0000000000000000000000000000000100000000000000000000000000000000
    1.061e-314  0000000000000000000000000000000010000000000000000000000000000000
  5.30499e-315  0000000000000000000000000000000001000000000000000000000000000000
  2.65249e-315  0000000000000000000000000000000000100000000000000000000000000000
  1.32625e-315  0000000000000000000000000000000000010000000000000000000000000000
  6.63124e-316  0000000000000000000000000000000000001000000000000000000000000000
  3.31562e-316  0000000000000000000000000000000000000100000000000000000000000000
  1.65781e-316  0000000000000000000000000000000000000010000000000000000000000000
  8.28905e-317  0000000000000000000000000000000000000001000000000000000000000000
  4.14452e-317  0000000000000000000000000000000000000000100000000000000000000000
  2.07226e-317  0000000000000000000000000000000000000000010000000000000000000000
  1.03613e-317  0000000000000000000000000000000000000000001000000000000000000000
  5.18065e-318  0000000000000000000000000000000000000000000100000000000000000000
  2.59033e-318  0000000000000000000000000000000000000000000010000000000000000000
  1.29516e-318  0000000000000000000000000000000000000000000001000000000000000000
  6.47582e-319  0000000000000000000000000000000000000000000000100000000000000000
  3.23791e-319  0000000000000000000000000000000000000000000000010000000000000000
  1.61895e-319  0000000000000000000000000000000000000000000000001000000000000000
  8.09477e-320  0000000000000000000000000000000000000000000000000100000000000000
  4.04739e-320  0000000000000000000000000000000000000000000000000010000000000000
  2.02369e-320  0000000000000000000000000000000000000000000000000001000000000000
  1.01185e-320  0000000000000000000000000000000000000000000000000000100000000000
  5.05923e-321  0000000000000000000000000000000000000000000000000000010000000000
  2.52962e-321  0000000000000000000000000000000000000000000000000000001000000000
  1.26481e-321  0000000000000000000000000000000000000000000000000000000100000000
  6.32404e-322  0000000000000000000000000000000000000000000000000000000010000000
  3.16202e-322  0000000000000000000000000000000000000000000000000000000001000000
  1.58101e-322  0000000000000000000000000000000000000000000000000000000000100000
  7.90505e-323  0000000000000000000000000000000000000000000000000000000000010000
  3.95253e-323  0000000000000000000000000000000000000000000000000000000000001000
  1.97626e-323  0000000000000000000000000000000000000000000000000000000000000100
  9.88131e-324  0000000000000000000000000000000000000000000000000000000000000010
  4.94066e-324  0000000000000000000000000000000000000000000000000000000000000001
             0  0000000000000000000000000000000000000000000000000000000000000000

6.8 Spezielle Werte

Die Gleitkommaarithmetik kennt einige spezielle Werte, z.B.

nextfloat(floatmax(Float64))
Inf
for x ∈ (NaN, Inf, -Inf, -0.0)
    @printf "%6g  " x
    printbitsf64(x)
end
   NaN  0111111111111000000000000000000000000000000000000000000000000000
   Inf  0111111111110000000000000000000000000000000000000000000000000000
  -Inf  1111111111110000000000000000000000000000000000000000000000000000
    -0  1000000000000000000000000000000000000000000000000000000000000000
  • Ein Exponentenüberlauf (overflow) führt zum Ergebnis Inf oder -Inf.
2/0, -3/0, floatmax(Float64) * 1.01, exp(1300)
(Inf, -Inf, Inf, Inf)
  • Damit kann weitergerechnet werden:
-Inf + 20, Inf/30, 23/-Inf, sqrt(Inf),  Inf * 0, Inf - Inf
(-Inf, Inf, -0.0, Inf, NaN, NaN)
  • NaN (Not a Number) steht für das Resultat einer Operation, das undefiniert ist. Alle weiteren Operationen mit NaN ergeben ebenfalls NaN.
0/0, Inf - Inf, 2.3NaN, sqrt(NaN)
(NaN, NaN, NaN, 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().
x = 0/0
y = Inf - Inf
@show x==y   NaN==NaN  isfinite(NaN) isinf(NaN) isnan(x) isnan(y);
x == y = false
NaN == NaN = false
isfinite(NaN) = false
isinf(NaN) = false
isnan(x) = true
isnan(y) = true
  • Es gibt eine „minus Null“. Sie signalisiert einen Exponentenunterlauf (underflow) einer betragsmäßig zu klein gewordenen negativen Größe.
@show 23/-Inf   -2/exp(1200)    -0.0==0.0;
23 / -Inf = -0.0
-2 / exp(1200) = -0.0
-0.0 == 0.0 = true

6.9 Mathematische Funktionen

Julia verfügt über die üblichen mathematischen Funktionen
sqrt, exp, log, log2, log10, sin, cos,..., asin, acos,..., sinh,..., gcd, lcm, factorial,...,abs, max, min,...,

darunter z.B. die Rundungsfunktionen

  • floor(T,x) = \(\lfloor x \rfloor\)
  • ceil(T,x) = \(\lceil x \rceil\)
floor(3.4),   floor(Int64, 3.5),   floor(Int64, -3.5)
(3.0, 3, -4)
ceil(3.4),    ceil(Int64, 3.5),    ceil(Int64, -3.5)
(4.0, 4, -3)

Es sei noch hingewiesen auf atan(y, x), den Arkustangens mit 2 Argumenten, 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)
atan(3, -2),    atan(-3, 2),    atan(-3/2)
(2.1587989303424644, -0.982793723247329, -0.982793723247329)

6.10 Umwandlung Strings \(\Longleftrightarrow\) Zahlen

Die Umwandlung ist mit den Funktionen parse() und string() möglich.

parse(Int64, "1101", base=2)
13
string(13, base=2)
"1101"
string(1/7)
"0.14285714285714285"
string(77, base=16)
"4d"

Zur Umwandlung der numerischen Typen ineinander kann man die Typnamen verwenden. Typenamen sind auch Konstruktoren:

x = Int8(67)
@show x   typeof(x);
x = 67
typeof(x) = Int8
z = UInt64(3459)
0x0000000000000d83
y = Float64(z)
3459.0

6.11 Literatur

  1. Quelle: Codekaizen, CC BY-SA 4.0, via Wikimedia Commons ↩︎