namespace cpp

C++ lernen, kennen, anwenden

Benutzer-Werkzeuge

Webseiten-Werkzeuge


numerik:maschinenzahlen

Maschinenzahlen

The aim of computing is insight, not numbers.
— R.W. Hamming

Aus der Schulmathematik sind natürliche, ganze, gebrochene, reelle und komplexe Zahlen bekannt. In der üblichen Schreibweise gilt: $\mathbb{N} \subset \mathbb{Z} \subset \mathbb{Q} \subset \mathbb{R} \subset \mathbb{C}$. Der Übergang zu immer umfassenderen Zahlbereichen wurzelt in dem Wunsch, im kleineren Zahlbereich unlösbare Aufgaben lösen zu können:

$x+5 = 2$ unlösbar in $\mathbb{N} \leadsto \mathbb{Z}$
$x*5 = 2$ unlösbar in $\mathbb{Z} \leadsto \mathbb{Q}$
$x^2 = 2$ unlösbar in $\mathbb{Q} \leadsto \mathbb{R}$
$x^2 = -1$ unlösbar in $\mathbb{R} \leadsto \mathbb{C}$

Umso erstaunlicher ist für den Neuling, dass Computer nur mit einem Teilbereich der ganzen bzw. gebrochenen Zahlen umgehen können. Und auch das nur eingeschränkt. Warum? Nun, Speicher ist endlich.

Ganzzahlen

Gott schuf die ganzen Zahlen, der Rest ist Menschenwerk.
— Leopold Kronecker

Die Menge natürlicher Zahlen $\mathbb{N}$ ist unbegrenzt, es gibt unendlich viele. Computerspeicher ist aber endlich, daher stehen zur Zahldarstellung nur begrenzt viele Stellen ($b$ Bits) zur Verfügung. Da Berechnungen mehr als eine Zahl erfordern und für viele Berechnungen Zahlen mit kleinen Beträgen genügen, ist der Wertebereich von Ganzzahlen im Computer eingeschränkt:

kronecker.cpp
#include <iostream>
 
int main()
{
  unsigned int n = 0;        // Peano-Axiom 1: 0 ist eine natürliche Zahl.
  do { ++n; } while (n > 0); // Peano-Axiom 2: Jede natürliche Zahl hat einen Nachfolger.
 
  auto min = n;
  auto max = n-1;
  std::cout << min << "..." << max << '\n';
}

Je nach Maschine dauert es eine kleine Weile, diese "Endlosschleife" zu durchlaufen. Abhängig von Prozessor, Betriebssystem, Compiler ist max = 65767, ca. 4 Mrd oder ein anderer Wert $2^b-1$ mit $b%$=16, 32, 64, 128 o.a. Rechenoperationen können nicht aus Menge maschinell darstellbarer Ganzzahlen $\mathbb{G} \subseteq \mathbb{Z}$ herausführen. So kommt es zu befremdlichen Ergebnissen. Die Summe positiver ganzer Zahlen, z.B. n+1, kann eigentlich nicht kleiner sein als die größere der beiden Zahlen. Dieser Überlauf des Wertebereichs wird in obigem Programm ausgenutzt, genauso wie der Unterlauf n-1 beim Unterschreiten der kleinsten darstellbaren Zahl. Willkommen in der Welt der Modulo-Arithmetik. Wird das höchstwertige Bit als Vorzeichen interpretiert (Datentyp int mit Zweierkomplement-Darstellung), so ist die Hälfte der Zahlen negativ: $-2^{b-1}\ldots -1$, der Rest $0\ldots 2^{b-1}-1$ nichtnegativ.1)

Gleitkommazahlen

Nobody wants to be a zero… Almost everybody wants to be a number one…
But the problem with these two number is, that they're just too close, leaves very little room in there for everybody else.
—- Laurie Anderson

Durch Abspalten von Zehnerpotenzen lässt sich die gleiche Zahl auf verschiedene Weisen \begin{align*} & 299792,458 \\ & 299,792458\cdot 10^3 \\ & 2,99792458\cdot 10^5 \approx 3,00\cdot 10^5 \end{align*} darstellen, bei denen das Komma (Gleitkomma, engl. floating point). Mit nur einer Ziffer vor dem Komma heißt die Darstellung normalisiert. Bei sehr großen und sehr kleinen Zahlen ist das vorteilhaft: $3,00\cdot10^5$ statt $300000$ mit 5 Stellen nach und $1,602\cdot 10^{-19}$ statt $0,0\ldots01602$ mit 19 Nullen vor der höchstwertigen von Null verschiedenen Ziffer. Neben der Größenordnung ist die Anzahl der Bedeutung tragenden (signifikanten, "zählenden" oder "gültigen") Ziffern sofort erkennbar. Wie viele Stellen sinnvoll oder "ausreichend" sind, hängt vom Anwendungsgebiet ab.

Im englischen Sprachraum wird die Schreibweise 3.00E5 "scientific notation" genannt. Dabei steht E oder e für Exponent zur Basis 10.

Eine Definition

Eine Gleitkommazahl $x$ zur Basis $B$ mit dem Wert $x = \pm \sum_{i=0}^n z_i B^{E-i}$ hat die Darstellung $\pm (z_0, z_1 z_2 \ldots z_n)_B \cdot B^E$ mit $n+1$ zählenden Stellen $z_i \in \{0, \ldots, B-1 \}$. $(z_0, z_1 z_2 \ldots z_n)_B$ wird Mantisse (oder Signifikand) genannt. Falls $z_0 \neq 0$, heißt $x$ normalisiert.

Endliche Mantissen und Exponenten

Für das Vorzeichen $s$ muss ein Bit zur Verfügung stehen: $\pm 1 = (-1)^s$.

Bei begrenzter Anzahl von Nachkommastellen $n$ erlaubt diese Darstellung nur endliche Brüche zur Basis $B$. Irrationale Zahlen und periodische Brüche wie $\pi = 3,14159\ldots_{10}$ und $\frac{1}{7} = 0,\overline{142857}_{10}$ sind ausgeschlossen. Ein gemeiner Bruch hat nur dann endlich viele Stellen, wenn sein Nenner durch Erweitern oder Kürzen in eine Potenz der Basis $B$ überführt werden kann. Manchmal läuft ein Basiswechsel unserer gewohnten Vorstellung zuwider: $\frac{1}{10} = 0,1_{10} = 0,0\overline{0011}_2$ ist nicht als endliche Binärzahl darstellbar (http://www.floating-point-gui.de). Dennoch ist dieser Basiswechsel nötig, wenn Computer Zahlen zur Basis 2 speichern — und die gebräuchlichsten tun dies.

Für normalisierte binäre Zahlen ($B=2$) gilt $z_0 = 1$. Als Mantisse $1,m = (1,z_1 z_2 \ldots z_n)_2$ dient dann die Bitfolge der Nachkommaziffern. Die führende Eins muss nicht extra gespeichert werden.

Der Wertebereich des Exponenten $E$ muss für die Speicherung ebenfalls eingeschränkt werden: $E_{min} \leq E \leq E_{max}$. (Hardware-)Implementierer haben es einfacher, wenn negative Zahlen im Exponenten vermieden werden können. Daher arbeiten sie mit einem verschobenen Exponenten (biased exponent) $e = E+v$ und einem festgelegten Versatz (engl. bias) $v = E_{max}$ statt $E = e -v$.2)

Bitmuster mit dem niedrigsten $e = 0$ stellen den Zahlbereichs-Unterlauf (Null und subnormale Zahlen mit $z_0 = 0$) dar. Das erlaubt einen sanfteren Übergang zur Null. Bitmuster mit dem höchstmöglichen $e = 2v+1$ kodieren den Überlauf ($\pm \infty$) und nicht als Zahl deutbare Bitmuster (NaN = not a number), die der Fehlerdiagnose dienen. So können Situationen wie 0/0, $\infty-\infty$, $\infty/\infty$, $0\cdot\infty$ oder unzulässige Funktionsargumente wie $\arcsin(2)$ abgebildet werden.

Zahlenbasis, Mantissenlänge und Exponentenbereich der Gleitkommazahlen können bei jeder Maschine anders sein. Das daraus entstehende Durcheinander führte ab 1977 zum Entwurf eines Standards.

IEEE 754-2008

Der Standard IEEE 754 in der aktuellen Fassung schlägt einige Gleitkommaformate vor. Konforme Implementierungen stellen wenigstens eine davon bereit:

  • einfach genaue Gleitkommazahl (binary32, single, float in C++) mit 8-Bit-Exponent und 23 Nachkommabits (siehe Tabelle unten),
  • doppelt genaue Gleitkommazahl (binary64,double, double in C++) mit 11-Bit-Exponent und 52 Nachkommabits
  • sowie weitere Formate: binary16, binary128, decimal32, decimal64, decimal128 u.a.
Bitmuster s[...e..][..........m..........] Wert
NaN 011111111xxxxxxxxxxxxxxxxxxxxxxx "not a number"
Überlauf 01111111100000000000000000000000 $+\infty$
$\nearrow$ max. 01111111011111111111111111111111 $1,11111111111111111111111_2\cdot2^{254-127}$
normalisiert 00111111100000000000000000000000 $1,m\cdot2^{e-127} = 1.0$
$\searrow$ min. 00000000100000000000000000000000 $1,00000000000000000000000_2\cdot2^{1-127}$
$\nearrow$ max. 00000000011111111111111111111111 $0,11111111111111111111111_2\cdot2^{-126}$
subnormal 000000000xxxxxxxxxxxxxxxxxxxxxxx $0,m\cdot2^{-126}$
$\searrow$ min. 00000000000000000000000000000001 $0,00000000000000000000001_2\cdot2^{-126}$
+Null 00000000000000000000000000000000 $0,0$

Mit höchstwertigem Bit 1 werden entsprechend negative Zahlen dargestellt. Als Kuriosität existieren zwei Bitmuster +0 und -0 für die Null, die aber als wertgleich behandelt werden. NaNs werden als ungeordnet betrachtet.

Rundung

Nicht als Maschinenzahl darstellbare reelle Zahlen $x$ müssen gerundet werden. Die Art der Rundung ist einstellbar (zur nächstliegenden Maschinenzahl, bei "Komma 5" wahlweise zur nächsten geraden oder zur größeren Ziffer; in Richtung $+\infty$, in Richtung $-\infty$ oder zur Null hin).

Wird diese Rundung $\text{rd}(x)$ zur nächstliegenden Maschinenzahl hin ausgeführt, beträgt die maximale absolute Abweichung $|x- \text{rd}(x)| = \frac{1}{2} \text{ulp}(x)$ die Hälfte des Abstandes zweier Maschinenzahlen in der Umgebung von $x$. Die Auflösung, d.h. der Unterschied zwischen $z_n = 0$ und $z_n = 1$ (engl. unit of least precision) $\text{ulp}(x) = B^{E-n}$ hängt von der Größenordnung der Zahl ab. Die Gleitkommazahlen sind nicht gleichmäßig über den gesamten darstellbaren Bereich verteilt. Betragsmäßige größere Zahlen haben einen größeren Abstand voneinander. Die relative Auflösung einer normalisierten Zahl zur Basis 2 bleibt dadurch im Bereich $2^{-(n+1)} < \frac{\text{ulp(x)}}{x} = \frac{2^{E-n}}{1,m\cdot2^E} \leq 2^{-n}$. Bei einfacher "Genauigkeit" sind das $2^{-23} \approx 1,2\cdot10^{-7}$, das entspricht etwa 6 bis 7 Dezimalstellen, bei doppelter "Genauigkeit" $2^{-52} \approx 2,2\cdot10^{-16}$, etwa 15 bis 16 Dezimalen.

Das Wort "Genauigkeit" (engl. precision) wurde hier bewusst in Anführungsstriche gesetzt, da es hier irreführend benutzt wird. Auflösung und Genauigkeit sollten nicht verwechselt werden. Die Auflösung begrenzt die maximal erreichbare Genauigkeit einer Maschinenzahl. Die Genauigkeit einer nur näherungsweise angegebenen Zahl $\tilde{x}$ kann jedoch viel geringer sein, wenn wir uns z.B. mit der Näherung 3,14 für die Kreiszahl begnügen. Ebenso führen Rechenoperationen aus dem Bereich der Maschinenzahlen heraus. Die Ergebnisse sind wieder zu einer Maschinenzahl zu runden, so dass sich Rundungsfehler akkumulieren können. Die Auflösung täuscht dann eine nicht vorhandene Genauigkeit vor.

Als Folge der zwischen den einzelnen Operationen ausgeführten Rundungen gehorcht die Maschinenarithmetik nicht dem Assoziativgesetz, wie ein Beispiel mit nur zwei signifikanten Dezimalen verdeutlicht: \begin{align*} (0,17 + 0,17) + 3,0 &= 0,34 + 3,0 \approx 3,3 \\ 0,17 + (0,17 + 3,0) &\approx 0,17 + 3,2 \approx 3,4 \end{align*}

Fehlerfortpflanzung

Seien $\varepsilon_x = (x-\tilde{x})/x$ und $\varepsilon_y = (y-\tilde{y})/y$ die relativen Fehler der Maschinenzahlen $\tilde{x} = (1-\varepsilon_x)x$ und $\tilde{x} = (1-\varepsilon_y)y$ zu den "exakten" Werten $x$ und $y$. Damit lassen sich die absoluten und relativen Fehler für die Ergebnisse der Grundrechenarten abschätzen.3)

Addition (Subtraktion als Addition von Zahlen mit unterschiedlichen Vorzeichen): \begin{align*} \text{rd}(\tilde{x}+\tilde{y}) &= x + y - \varepsilon_x x - \varepsilon_y y \pm \frac{1}{2} \text{ulp}(x+y) \\ \varepsilon_{x+y} = \frac{x+y - \text{rd}(\tilde{x}+\tilde{y})}{x+y} &= \varepsilon_x \frac{x}{x+y} + \varepsilon_y \frac{y}{x+y} \mp \frac{1}{2} \frac{\text{ulp}(x+y)}{x+y} = \varepsilon_x \frac{x}{x+y} + \varepsilon_y \frac{y}{x+y} \mp \frac{1}{2} B^{-n} \end{align*}

Multiplikation: \begin{align*} \text{rd}(\tilde{x}\cdot\tilde{y}) &= (1 - \varepsilon_x - \varepsilon_y + \varepsilon_x \varepsilon_y)xy \pm \frac{1}{2} \text{ulp}(x\cdot y) \\ \varepsilon_{x\cdot y} = \frac{x\cdot y - \text{rd}(\tilde{x}\cdot \tilde{y})}{x\cdot y} &= \varepsilon_x + \varepsilon_y - \varepsilon_x \varepsilon_y \mp \frac{1}{2} \frac{\text{ulp}(x\cdot y)}{x\cdot y} \approx \varepsilon_x + \varepsilon_y \mp \frac{1}{2} B^{-n} \end{align*}

Division: \begin{align*} \text{rd}(\tilde{x}/\tilde{y}) &= \frac{1 - \varepsilon_x}{1 - \varepsilon_y} \frac{x}{y} \pm \frac{1}{2} \text{ulp}(x/y) \\ \varepsilon_{x/y} = \frac{x/y - \text{rd}(\tilde{x}/\tilde{y})}{x/y} &= \frac{\varepsilon_x- \varepsilon_y}{1 - \varepsilon_y} \mp \frac{1}{2} \frac{\text{ulp}(x/y)}{x/y} \approx \varepsilon_x- \varepsilon_y \mp \frac{1}{2} B^{-n} \end{align*}

Bei kleinen relativen Fehlern sind die Vereinfachungen bei Multiplikation und Division gerechtfertigt. Beide Operationen sind gutartig in dem Sinne, dass kleine Fehler der Eingabewerte nicht zu großen Fehlern im Ergebnis führen. Addition von Zahlen $x + y \approx 0$ und Subtraktion von Zahlen $x - y \approx 0$ führen dagegen zu großen relativen Fehlern (Auslöschungseffekt) bis hin zu völlig sinnlosen Resultaten.

Für die häufig zusammen vorkommende Multiplikation und Addition steht die Funktion fma(a,x,y) = ax+y (fused multiply-add) bereit, die eine Rundung einspart:

\begin{align*} \varepsilon_{ax + y} &\approx (\varepsilon_a + \varepsilon_x) \frac{ax}{ax+y} + \varepsilon_y \frac{y}{ax+y} \mp \frac{1}{2} B^{-n} \end{align*}

Kondition

Der absolute Fehler einer (als exakt vorausgesetzten) Funktionswertberechnung $f(x) - f(\tilde{x}) \approx f'(x) \cdot(x-\tilde{x})$ wird vor allem durch den Anstieg $f'$ in der Nähe der Werte $x$ und $\tilde{x}$ bestimmt. Für den relativen Fehler \begin{align*} \varepsilon_{f(x)} &= \frac{f(x) - f(\tilde{x})}{f(x)} = \frac{f(x) - f(x -\varepsilon_x x)}{\varepsilon_x x} \cdot \frac{\varepsilon_x x}{f(x)} \\ &= \frac{f'(x) \cdot x}{f(x)} \cdot \varepsilon_x = K_{f(x)} \cdot \varepsilon_x \end{align*} ist die Konditionszahl $K_{f(x)} = \frac{f'(x) \cdot x}{f(x)}$ maßgeblich. Sie gibt an, um welchen Faktor der relative Fehler der Eingangsgröße $x$ verstärkt wird.4) Eine Funktion wie $\sqrt{x}$ verhält sich überall gutartig: $K_{\sqrt{x}} = \frac{x}{\sqrt{x} \cdot 2 \sqrt{x}} = 1/2$. Sie halbiert den Fehler. Dagegen führt die Berechnung von $f(x) = (x-1)^3$ in der Nähe der Nullstelle $x = 1 + \varepsilon$ wegen $K_{f(x)} = \frac{3(x-1)^2\cdot(1+\varepsilon)}{(x-1)^3} \approx 3/\varepsilon$ bei $\varepsilon = 10^{-5}$ zu einem 300000fachen Fehler, d.h. mehr als 5 Dezimalstellen gehen verloren. Aufgaben mit sehr großen Konditionszahlen heißen schlecht gestellt.

Auswege

Für ganzzahlige Berechnungen mit großer Stellenanzahl, z.B. in der Kryptarithmetik, gibt es Bibliotheken.

Die Nachverfolgung der Rundungsfehler bei Gleitkommaberechnungen per Hand gestaltet sich schnell unübersichtlich. Aufgabe der numerischen Mathematik ist, trotz dieser Unzulänglichkeiten stabile Rechenvorschriften zu finden, die den unvermeidbaren Fehler nicht unnötig vergrößern.

Obwohl schon seit längerem Vorschläge existieren, die Ergebnisverifikation von Gleitkommaberechnungen zu automatisieren (Einschluss des korrekten Ergebnisses durch Intervallarithmetik), werden diese nicht flächendeckend eingesetzt. Offensichtlich ist die Mehrzahl der Computernutzer mehr an schnellen als an richtigen Ergebnissen interessiert.

1)
Der Überlauf vorzeichenbehafteter Ganzzahlen ist laut C++-Standard undefiniertes Verhalten. Weglassen von signed in obigem Programm und Übersetzen mit g++ -O2 bzw. g++ -O3 führt tatsächlich zu einer Endlosschleife, wenn beim GNU-Compiler nicht die Option -fwrapv angegeben wird.
2)
Deutschsprachige Lehrbücher bezeichen den verschobenen Exponenten zuweilen als Charakteristik.
3)
Zur Vereinfachung wird hier der schlimmste Fall angenommen, dass jede Rundung den maximalen Fehler erzeugt. Ist der Rundungsfehler klein gegenüber den relativen Fehlern der Ausgangswerte, kann dieser vernachlässigt werden. Dann ergeben sich die Formeln aus [Gerhard Opfer: Numerische Mathematik für Anfänger. Vieweg, Braunschweig (2002), S. 7f]
4)
Wirtschaftler nennen diese Größe Elastizität $\frac{dy/y}{dx/x}$: Um wieviel Prozent ändert sich die Ausgangsgröße bei einer Änderung der Eingangsgröße um ein Prozent?
numerik/maschinenzahlen.txt · Zuletzt geändert: 2014-12-29 13:27 (Externe Bearbeitung)