1 Theorie | |
→ | 1.1 Die Programmiersprache C |
gespeichert.
Siehe Anhang Gleitkommazahlen für Details.
Sowohl das Exponentenfeld als auch das Mantissenfeld sind in der
Bitanzahl begrenzt.
Die Begrenzung des Exponentenfeldes führt dazu, dass die
darstellbaren Zahlen betragsmäßig nach oben, aber auch nach unten
begrenzt sind.
Die Begrenzung des Mantissenfeldes führt dazu, dass die Genauigkeit
begrenzt ist.
Der kontinuierliche Zahlenbereich der reellen Zahlen wird rechnerintern durch einen Satz diskreter Werte abgebildet. Rechnerintern wird also ein Näherungswert des "richtigen" Wertes gespeichert.
Für Gleitkommazahlen sind spezielle Werte definiert:
Gleitkommaoperationen werden nicht in der ALU (arithmetic logic
unit) des Prozessors sondern in der FPU (floating point unit)
durchgeführt.
Diese war bei früheren Prozessoren ein Zusatz-IC (mathematischer
Coprozessor), bei modernen Prozessoren ist die FPU mit in den
Prozessor integriert.
Bei Division durch 0 oder betragsmäßig sehr kleine Werte erfolgt
kein sofortiger Programmabbruch, stattdessen ergibt sich das
Ergebnis INF oder NAN.
Als normalisierte Darstellung wird bezeichnet, wenn das gesamte Mantissenfeld benutzt wird.
Sind Zahlen betragsmäßig zu klein, um das
gesamte Mantissenfeld zu benutzen, spricht man von
nicht-normalisierter Darstellung.
Durch die Verwendung von weniger signifikanten Bits ist die
Darstellung u.U. erheblich ungenauer als in der normalisierten
Darstellung.
In C stehen folgende Datentypen für Gleitkommazahlen zur Verfügung
Im Beispiel soll die Zahl 4,75 (dezimal) als double dargestellt
werden.
Die Umwandlung in das Dualsystem ergibt 100,11
(1*22+0*21+0*20+1*2-1+1*2-2).
Die normalisierte Darstellung (Komma verrutscht und Exponent zur
Basis 2 angepasst, so dass das Komma nach der führenden 1 steht)
ergibt 1,0011*22.
Da es sich um eine nichtnegative Zahl handelt, hat das
Vorzeichenfeld den Wert 0.
Zum Exponenten 2 wird der BIAS 1023 hinzuaddiert, daraus wird ein
11 Bit breites Exponentenfeld 10000000001 gebildet.
Bei der Bildung des 52 Bit breiten Mantissenfeldes wird die
führende 1 weggelassen, die auf das Komma folgenden Bits ergeben
den Anfang des Mantissenfeldes. Hinten wird mit Nullen aufgefüllt.
Somit ergibt sich das Mantissenfeld
0011000000000000000000000000000000000000000000000000.
Insgesamt ergibt sich der double-Wert zu
0100000000010011000000000000000000000000000000000000000000000000.
Unter Windows mit Visual Studio kann mit
int _finite(double x);
getestet werden, ob eine Zahl x weder INF noch NAN ist. Falls ja, handelt es sich um eine gültige Zahl.
Auf POSIX-Systemen wird hierzu
int isfinite(x)
verwendet, dies ist meist als Makro implementiert.
Unter Windows kann mit der Funktion _fpclass() ein Gleitkommawert klassifiziert werden. Die Ergebniswerte _FPCLASS_ND oder _FPCLASS_PD zeigen nicht-normalisierte negative bzw. positive Werte an.
Auf POSIX-Systemen kann das Makro
int isnormal(x)
benutzt werden, um zu prüfen, ob normalisierte Darstellung vorliegt.
Liefert eine Berechnung mehrere Ergebnisse (z.B. bei einer Matrizenmultiplikation), so möchte man nicht jeden Ergebniswert einzeln mit _finite() bzw. isfinite() überprüfen.
Treten bei einer Berechnung Probleme auf, wird eine Gleitkomma-Exception ausgelöst. Dabei wird in einem speziellen Statuswort der Gleitkomma-Einheit ein Bit gesetzt. Dieses Statuswort kann abgefragt und auch zurückgesetzt werden.
Gleitkomma-Exceptions haben nichts mit den aus Java und C++ bekannten Exceptions zu tun, die Abbruch und Ausnahmebehandlung auslösen, es wird hier lediglich derselbe Begriff verwendet.
Unter Windows mit Visual Studio stehen die Funktionen _clearfp() (Gleitkomma-Exceptions im Statuswort löschen) und_statusfp() (Statuswort auslesen) aus der Header-Datei float.h zur Verfügung.
Auf POSIX-Systemen werden stattdessen die Funktionen feclearexcept() (Statuswort löschen) und fetestexcept() (Statuswort auslesen) aus der Header-Datei fenv.h verwendet. Durch Optimierungen des Compilers können die Ergebnisse dieser Tests verfälscht werden. Mit
#pragma STDC FENV_ACCESS ON
werden derartige Optimierungen verhindert.
Sind Teile der Berechnung bzw. die gesamte Berechnung in andere
Module ausgelagert, sollte die #pragma-Anweisung in jedem an der
Berechnung beteiligten Modul stehen.
Vor einer evtl. auch umfangreichen Berechnung werden die
Exceptions zurückgesetzt. Nach der Berechnung kann dann geprüft
werden, ob in der Berechnung Probleme auftraten.
Folgende Gleitkomma-Exceptions können auftreten:
Name (Windows) |
Name (POSIX) |
Bedeutung |
---|---|---|
_SW_ZERODIVIDE | FE_DIVBYZERO | Division durch 0 |
_SW_INVALID | FE_INVALID | Unbestimmtes Ergebnis, z.B. 0/0 oder ∞-∞ |
_SW_OVERFLOW | FE_OVERFLOW | Zahlenbereichsüberlauf, Ergebnis ist betragsmäßig zu groß für Datentyp |
_SW_UNDERFLOW | FE_UNDERFLOW | Zahlenbereichsunterlauf, Ergebnis ist betragsmäßig zu klein für Datentyp. Diese Exception wird auch ausgelöst, wenn alle Operanden normalisiert darstellbar sind, das Ergebnis einer Operation aber nur nichtnormalisiert darstellbar ist. |
_SW_INEXACT | FE_INEXACT | Ungenaues Ergebnis, wird meist gemeinsam
mit overflow oder underflow gesetzt. Diese Exception wird auch
gesetzt, wenn ein Ergebnis mehr Nachkommabits aufweist, als im
Mantissenfeld gespeichert werden können. Da diese Exception bei nichttrivialen Berechnungen wie z.B. "1/3" oder "1/10" (liefert im Dualsystem 0,000110011...0011...) häufig vorkommt, wird im Normalfall nicht auf diese Exception getestet. |
_SW_DENORMAL | Ein Ergebnis einer Operation ist
betragsmäßig so klein, dass keine normalisierte Darstellung mehr
möglich ist. Dies geht meist mit dem Abschneiden einer erheblichen
Bitanzahl daher und führt somit zu einem ungenauen Ergebnis. Hinweis: Durchlauf eines Testprogrammes zeigte, dass diese Exception nicht ausgelöst wurde, auch wenn nicht normalisierte Ergebnisse bei einer Berechnung auftraten. |
#include <float.h> #include <math.h> /** Gleitkomma-Exceptions. */ unsigned int exc = 0; /** Bit-Maske zum Pruefen der Exceptions. */ const unsigned int mask = _SW_ZERODIVIDE | _SW_INVALID | _SW_OVERFLOW | _SW_UNDERFLOW | _SW_DENORMAL; ... /* Alle Exceptions zuruecksetzen */ (void)_clearfp(); /* ... Berechnung durchfuehren */ /* Exceptions in exc speichern, alle Exceptions zuruecksetzen */ exc = _clearfp(); if (0 != (mask & exc)) { /* Es ist eine Exception aufgetreten */ /* Bei Bedarf genauer pruefen, welche Exception(s) */ if (0 != (_SW_ZERODIVIDE & exc)) { /* Division durch 0 */ } if (0 != (_SW_INVALID & exc)) { /* Ungueltiger Operand/Argument, z.B. Wurzel aus negativer Zahl */ } if (0 != (_SW_OVERFLOW & exc)) { /* Zahlenbereichsueberlauf */ } if (0 != (_SW_UNDERFLOW & exc)) { /* Zahlenbereichsunterlauf */ } if (0 != (_SW_DENORMAL & exc)) { /* Denormalisierte Darstellung notwendig, Genauigkeitsverlust */ } }
Die Beispieldatei ex072.c erweitert das Programm zur
Quadratflächenberechnung so, dass für Windows eine Abfrage nach
Gleitkomma-Exceptions stattfindet.
Zum Test können neben "normalen" auch extrem große Zahlen wie
1.5e200 oder extrem kleine Zahlen wie 1.5e-200 verwendet
werden.
#include <fenv.h> #include <float.h> #include <math.h> #pragma STDC FENV_ACCESS on /** Gleitkomma-Exceptions. */ int exc = 0; ... /* Alle Exceptions zuruecksetzen */ feclearexcept(FE_ALL_EXCEPT); /* ... Berechnung durchfuehren */ /* Exceptions in exc speichern, alle Exceptions zuruecksetzen */ exc = fetestexcept(FE_ALL_EXCEPT & (~(FE_INEXACT))); feclearexcept(FE_ALL_EXCEPT); if (0 != (exc & FE_ALL_EXCEPT)) { /* Es ist eine Exception aufgetreten */ /* Bei Bedarf genauer pruefen, welche Exceptions */ if (0 != (exc & FE_DIVBYZERO)) { /* Division durch 0 */ } if (0 != (exc & FE_INVALID)) { /* Ungueltiger Operand/Argument, z.B. Wurzel aus negativer Zahl */ } if (0 != (exc & FE_OVERFLOW)) { /* Zahlenbereichsueberlauf */ } if (0 != (exc & FE_UNDERFLOW)) { /* Zahlenbereichsunterlauf */ } }
Die Beispieldatei ex073.c erweitert das Programm zur Quadratflächenberechnung so, dass für POSIX eine Abfrage nach Gleitkomma-Exceptions stattfindet. Zum Test können neben "normalen" auch extrem große Zahlen wie 1.5e200 oder extrem kleine Zahlen wie 1.5e-200 verwendet werden.
Mit
#define _USE_MATH_DEFINES 1 #include <math.h>
wird die Header-Datei math.h eingeschlossen, die Konstanten,
Makros und Funktionsprototypen für double-Berechnungen enthält.
Die erste Zeile (#define...) wird nur unter Windows benötigt.
Die folgenden Konstanten, Makros und Funktionen sind verfügbar (hier eine unvollständige Darstellung):
Nicht alle der hier genannten Konstanten, Makros und Funktionen sind auf allen Systemen zu finden.
Unter Windows mit Visual Studio wird die Mathematik-Bibliothek bei Bedarf automatisch mit zum Projekt gelinkt.
Unter Linux/UNIX muss dem Linker die Option "-lm" übergeben werden.
Dieser Abschnitt basiert zum großen Teil auf dem Artikel "The Perils of Floating Point" von Bruce M. Bush, [BMB1996].
Wie oben bereits dargestellt, ist die Darstellung von
Gleitkommazahlen auf dem Computer sowohl im Wertebereich als auch
in der Genauigkeit beschränkt, da sowohl für das Exponentenfeld als
auch für das Mantissenfeld jeweils eine exakt festgelegte Anzahl
von Bits benutzt wird.
Der kontinuierliche Zahlenbereich der reellen Zahlen wird
rechnerintern durch einen Satz diskreter Werte abgebildet. Wird das
Ergebnis einer Berechnung in einer Gleitkommavariable gespeichert,
wird der diskrete Wert in der Variable abgelegt, der dem Ergebnis
am nächsten kommt.
Die Zahlenwerte werden im Dualsystem gespeichert, dies betrifft sowohl Vor- als auch die Nachkommastellen. Einige Brüche, die im Dezimalsystem endlich sind (z.B. ⅒=0,1), sind im Dualsystem unendlich (0,000110011).
In einem Beispielprogramm setzen wir eine double-Variable auf den Wert 4,2, ziehen anschließend 4 ab und vergleichen das Ergebnis mit 0,2.
#include <stdio.h>
int main(void)
{
double a, y;
y = 4.2;
a = y - 4.0;
printf("a = %g\n", a);
if (0.2 != a) {
printf("a ungleich 0.2\n");
}
else {
printf("a gleich 0.2\n");
}
return 0;
}
Der Wert 4,2 hat unendlich viele Nachkommastellen, die Darstellung im Dualsystem ergibt 100,00110011. In normierter Darstellung ergibt dies 1,0000110011*22. Im Mantissenfeld wird die führende "1," weggelassen, das Exponentenfeld entsteht aus Summe des Exponenten 2 und BIAS 1023. Für die Darstellung mit Vorzeichen-, Expontenten- und Mantissenfeld ergibt sich:
0 10000000001 0000110011001100110011001100110011001100110011001101
Die 1 im letzten Bit ergibt sich, weil in der unendlich langen Dualdarstellung an dieser Stelle eine 0 gefolgt von 11 steht, die Darstellung aber wegen der begrenzten Länge des Mantissenfeldes an dieser Stelle abbricht und Aufrunden an dieser Stelle einen geringeren Fehler produziert als Abrunden.
Der Wert 4 wird dual als 100,000... dargestellt, in normierter Darstellung ergibt dies 1,000...*22. Für Vorzeichen-, Exponenten- und Mantissenfeld ergibt sich:
0 10000000001 0000000000000000000000000000000000000000000000000000
Die Subtraktion beider Werte ergibt
1,0000110011001100110011001100110011001100110011001101*22 -1,0000000000000000000000000000000000000000000000000000*22 --------------------------------------------------------------------- 0,0000110011001100110011001100110011001100110011001101*22
Die Bildung einer normalisierten Darstellung führt zu
1,1001100110011001100110011001100110011001100110100000*2-3
Bei der Bitverschiebung im Mantissenfeld werden von hinten Nullen aufgefüllt. Für Vorzeichen-, Exponenten- und Mantissenfeld ergibt sich:
0 01111111001 001100110011001100110011001100110011001100110100000
Der dezimale Bruch 0,2 ergibt im Dualsystem den unendlich langen Bruch 0,00110011, in normierter Schreibweise 1,10011*2-3. Für Vorzeichen-, Exponenten- und Mantissenfeld ergibt dies:
0 01111111001 001100110011001100110011001100110011001100110011010
Hier werden die letzten beiden Bits des Mantissenfeldes durch das Aufrunden beim Abschneiden der nachfolgenden Stellen von "01" zu "10".
Die Gleitkommazahldarstellungen von (4,2 - 4) und 0,2 unterscheiden sich voneinander, der Test auf Gleichheit schlägt fehl.
Die Ausgaberoutinen der printf()-Familie gehen davon aus, dass nur eine bestimmte Anzahl signifikanter Stellen in der Ausgabe sinnvoll ist, deshalb wird als Berechnungsergebnis zwar 0,2 angezeigt während andererseits der Test auf Gleichheit mit 0,2 fehlschlägt.
Im diesem Beispiel wird der Wert 2251799813685248.375 in einer double-Variablen gespeichert. Anschließend werden davon 2251799813685248.0 subtrahiert und das Ergebnis ausgegeben.
#include <stdio.h>
int main(void)
{
double a;
a = 2251799813685248.375;
a = a - 2251799813685248.0;
printf("a = %g\n", a);
return 0;
}
Überraschenderweise ergibt sich als Ergebnis 0,5.
Die Dezimalzahl 2251799813685248.375 ergibt im Dualsystem
1000000000000000000000000000000000000000000000000000,011.
Die normierte Darstellung ergibt
1,000000000000000000000000000000000000000000000000000011*251.
Die Zahl ist nichtnegativ, das Vorzeichenfeld hat somit den Wert
0.
Das 11 Bit breite Exponentenfeld ergibt sich, indem zum Exponenten
51 der BIAS 1023 addiert wird zu 10000110010.
Von der Mantisse wird die führende 1 nicht in das Mantissenfeld
übernommen. Der Nachkommateil ist zu lang für das Mantissenfeld,
die letzten beiden 1-Bits passen nicht hinein. Das letzte noch
hineinpassende Bit ist die 0 davor, hier wird wegen der
weggelassenen folgenden 1-Bits zu 1 aufgerundet. Gespeichert ist
also der Wert
1,0000000000000000000000000000000000000000000000000001*251.
Wird hiervon
1,0000000000000000000000000000000000000000000000000000*251
abgezogen, verbleibt
0,0000000000000000000000000000000000000000000000000001*251,
was in normierter Darstellung
1,0000000000000000000000000000000000000000000000000000*2-1
ergibt.
Im Beispiel wird die Variable a auf den Wert
9007199254740992.0 gesetzt.
Für Variable b wird nochmals 1 hinzuaddiert.
Wider Erwarten ist allerdings b nicht größer als
a.
#include <stdio.h>
int main(void)
{
double a, b;
a = 9007199254740992.0;
b = a + 1.0;
if (b > a) {
printf("b ist groesser als a\n");
}
else {
printf("b ist nicht groesser als a\n");
}
return 0;
}
Im Exponenten- und Mantissenfeld von a wird die normierte
Darstellung
1,00000000000000000000000000000000000000000000000000*253
gespeichert. Die zu addierende 1 wird normiert als
1,0000000000000000000000000000000000000000000000000000*20
dargestellt.
Vor der Addition werden beide Werte in eine Form mit gleichen
Exponenten (in diesem Fall 53) gebracht:
1,0000000000000000000000000000000000000000000000000000 * 253 + 0,00000000000000000000000000000000000000000000000000001 * 253 -------------------------------------------------------------------------- 1,0000000000000000000000000000000000000000000000000000 * 253
Bei der Addition werden nur die Teile der Mantisse berücksichtigt, die auch nach dem Gleichmachen der Exponenten noch in das Mantissenfeld passen. Aufgrund der Größe von a wandert das 1-Bit der 1.0 zu weit nach hinten, so dass es nicht mehr berücksichtigt wird.
Im Beispiel wird die Variable a auf den Wert ⅐ gesetzt.
Die Variable b wird auf 0 gesetzt, in jedem
Schleifendurchlauf wird a addiert, anschließend wird
b zur int-Variable i konvertiert.
Im 7. Schleifendurchlauf würden wir erwarten, dass b nicht mehr
kleiner als 1 ist und die Konvertierung zu int eine 1 ergibt.
Beim praktischen Test zeigt sich, dass zwar im 7. Durchlauf von
printf() für den Wert von b eine 1 angezeigt wird,
der Test allerdings noch immer sagt, dass b kleiner als 1
ist und die Konvertierung zu int eine 0 ergibt.
#include <stdio.h>
int main(void)
{
double a, b;
int i;
a = 1.0 / 7.0;
b = 0.0;
do {
b = b + a;
i = (int)b;
printf("Addition ergab b = %g\n", b);
printf("Konvertierung zu int ergibt i = %d\n", i);
if (b < 1.0) {
printf("b ist kleiner als 1.0\n\n");
}
} while (b < 1.0);
return 0;
}
Der dezimale Bruch ⅐ ergibt im Dualsystem den unendlichen Bruch
0,001001, in
normierter Darstellung 1,001*2-3. Bei der
Übernahme in das Mantissenfeld erfolgt das Abschneiden an einer
Stelle 0|01001, beim
Abrunden wird die 0 übernommen, der Rest wird abgeschnitten.
Der in a gespeicherte Wert ist somit geringfügig kleiner als
⅐.
Beim Aufsummieren wird nach 7 Durchläufen somit ein Wert erreicht,
der zwar sehr nahe an 1 liegt, aber noch immer geringfügig kleiner
ist.
Bei der Konvertierung zu int mittels Typecast wird der
Nachkommateil abgeschnitten, das Ergebnis ist also der betragsmäßig
nächstkleinere int-Wert, in diesem Fall 0.
Die Ausgabefunktionen der printf()-Familie geben nur eine
bestimmte Stellenanzahl aus und runden das Ergebnis vor der Ausgabe
so, dass es in diese Stellenanzahl passt. Daher die Ausgabe einer
1.
Die Beispiele veranschaulichen u.a. einen Teil der nachfolgenden Sachverhalte:
#define EPSILON 1.0e-5 if ( fabs(a-b) < (EPSILON*fabs(a+b)) ) { /* a und b sind naeherungsweise gleich */ } if ( (b-a) > (EPSILON*fabs(b)) ) { /* b ist groesser als a */ } if ( (a-b) > (EPSILON*fabs(b)) ) { /* b ist kleiner als a */ }
double rint(double x);aus math.h verwendet werden.
/* SO NICHT */ #define PI 3.141 #define SQRT2 1.414was natürlich zu ungenaueren Ergebnissen führt.
Quelle: CERT C Coding Standard 〈3〉
[BMB1996] | Bruce M. Bush: The Perils of Floating Point http://www.lahey.com/float.htm |
1 | http://www.netlib.org |
2 | https://software.intel.com/en-us/mkl |
3 | http://www.securecoding.cert.org |