Ein wenig Theorie zu den Fließkommazahlen
in der Softwaretechnik mit C und C++

Es hat sich meiner Erfahrung nach herumgesprochen, dass die Verwendung von Fließkommazahlen – auch Gleitkommazahlen genannt – in der Software­entwicklung nicht immer ganz problemfrei ist. Manch Programmierer scheut deshalb deren Verwendung wie der Teufel das Weihwasser. Andere Entwickler verwenden die ensprechenden Datentypen float und double etwas naiver und reagieren erst wenn die Software zu erkennen gibt, dass sie nicht so laufen möchte, wie es sich der Urheber eigentlich dachte. Doch was ist nun eigentlich das Problem mit diesen Zahlen?

#include <stdio.h>
int main()
{
      double d1 = 1.2;
      double d2 = 0.8;
      d2 += 0.4;
      puts( ( d1 == d2 ) ? "gleich" : "ungleich" );
      return 0;
}
Listing 1: Beispiel in C zur Ungenauigkeit eines
errechneten Gleitkommawertes.

Das Darstellungsproblem

Einfache Vergleiche von errechneten Fließkommazahlen unter Verwendung der Vergleichs­operatoren in C und C++ führen oft nicht zum erwarteten Ergebnis. Das C-Beispiel im nebenstehenden Kasten lässt sich mit jedem Standard C Compiler übersetzen. Zum Beispiel mit dem C Compiler unter UNIX/Linux:

cc bsp1.c -o bsp1
Der Aufruf ./bsp1 zeigt ungleich.

Das liegt daran, dass der Wert 1.2 wie die meisten Zahlen binär überhaupt nicht exakt darstellbar ist. Es wird ein Wert codiert, der knapp neben dem zu erwarteten Ergebnis liegt. Generell kann gesagt werden, dass keine Darstellungsform von Realen Zahlen in der Lage ist, alle Zahlen abzubilden. Auch nicht die dezimale Form, die wir für gewöhnlich im Alltag verwenden. Zwischen zwei darstellbaren Zahlen liegen immer unendlich viele nicht darstellbare. Dieses grundsätzliche Problem ist also keine Spezialität der Binärdarstellung in Computersystemen. Man kann die meisten Zahlen bestenfalls als Näherung in ausreichender Genauigkeit darstellen und nur wenige exakt. Darüber hinaus gibt es aber - theoretisch unendlich - viele Zahlen, die im Dezimalsystem exakt darstellbar sind, im Binärsystem jedoch unedlich viele Nachkomma­stellen haben. Die darstellbaren Zahlen des Dezimalsystems sind also nicht komplett auf das Binärsystem abbildbar und man muss mit Näherungswerten arbeiten. Weiteres ist im Abschnitt über die Ungenauigkeit der Fließkommadarstellung beschrieben.

Die technische Darstellung von Fließkommazahlen in C und C++

Die in C und C++ gebräuchlichen Gleitkommatypen float und double sind im Standard IEEE 754 definiert.

Darin ist geregelt, dass eine Zahl z nach dem Prinzip z = m × 2e dargestellt wird. Wobei m die Mantisse und e der Exponent ist. Der Datentyp float hat 32 Bit, die sich in eine Mantisse mit 23 Bit, einen Exponenten mit 8 Bit und ein Vorzeichenbit aufteilen.

Bei double sieht die Aufteilung genauso aus. Nur die Breiten ändern sich gegenüber float. Die Mantisse hat 52 und der Exponent 11 Bit. Zusammen mit dem Vorzeichenbit kommen wir auf eine Größe von 64 Bit.

Es gibt mit long double in den C und  C++ Standards noch einen breiteren Datentyp für die Darstellung von Gleitkommawerten. Allerdings ist dieser Datentyp nicht strikt definiert und kann auf unter­schied­lichen Plattformen oder bei der Verwendung unter­schied­licher Compiler in seiner Darstellung variieren. Übliche Darstellungs­formate für long double sind beispiels­weise das 80 Bit breite Extended Precision Format oder das 128 Bit breite Quadrupel Format aus dem Standard IEEE 754.
Da die technische Darstellung dieses Datentyps nicht einheitlich definiert ist, verzichte ich in diesem Artikel auf eine genauere Beschreibung. Gleichwohl gilt die Beschreibung des Problem­bereichs der Fließkomma­darstellung grundsätzlich auch für diesen breiteren Typ. Auch das Prinzip der Codierung ist technisch analog zu den beiden Datentypen float und double gelöst.

Wie werden nun die Gleikommazahlen codiert?

Die Normalisierung von Fließkommawerten

Wenn man Werte binär so darstellt, dass sie vor dem Komma genau eine Stelle mit der Eins haben, wird diese Darstellung normalisiert genannt. In der Binärdarstellung steht also immer die Eins an der ersten Stelle wenn es sich um den sogenannten normalisierten Zahlenbereich handelt und eine Null wenn es der denormalisierte Bereich ist. Der Wert 0.0 kann deshalb nicht normalisiert werden. Aber bei den meisten Zahlen, die in der Exponenten­darstellung darstellbar sind, ist eine Normalisierbarkeit möglich. Man verschiebt dabei das Komma mit Hilfe des Exponenten so, dass nur eine Stelle davor übrig bleibt. Die normalisierten Zahlen bewegen sich für float im Bereich FLT_MIN bis FLT_MAX im positiven, wie auch im negativen Bereich. Um die 0.0 herum ist der Bereich der nicht normalisierten oder eben denormalisierten Werte. Für double gibt es die Konstanten DBL_MIN und DBL_MAX.

Durch das Zusammenwirken von Mantisse und Exponent wird ein Wert dargestellt der durch das Vorzeichenbit nur noch mit einem Vorzeichen versehen wird. Der Exponent verschiebt einfach nur das Komma - Im binären Zahlensystem natürlich. Dabei gibt es noch eine Besonderheit bezüglich des Exponenten.

Die Verschiebedarstellung des Exponenten

Der Exponent wird in einer sogenannten Verschiebedarstellung repräsentiert (engl. Biased Representation). Dazu wird ein Verschiebewert − Bias oder Exzess − aufaddiert. Wenn man also den echten Exponenten eines Wertes erhalten möchte, muss man den entsprechenden Bias vom dargestellten Exponentenwert abziehen. Der Bias für den Exponenten errechnet sich aus 2n−1−1 für n = Bits des Exponenten in der binären Darstellung. Der Bias für float ist 27−1 = 127. Für double ist der Bias 210−1 = 1023.

Die Verschiebedarstellung führt dazu, dass sich Gleitkommazahlen bei größer/kleiner-Vergleichen wie ganze Zahlen vergleichen lassen. Der Exponent wird damit in eine Form gebracht, in der sein direkter Betrag verglichen werden kann. Bei der Mantisse ist das sowieso schon der Fall. Da im Bitmuster der Exponent die Mantisse nach vorne erweitert, können beide zusammen genommen direkt verglichen werden. Die Zahl, bei der am weitesten vorne eine 1 steht, während die andere an der gleichen Position eine 0 hat, ist die größere. Das ist technisch als Hardware­operation sehr einfach zu realisieren und deshalb auch sehr performant. Es entspricht einem Vergleich zwischen zwei ganz­zahligen Werten.

Sonderfälle

Die denormalisierten Zahlen

Wenn der Exponent nur aus Nullen besteht, spricht man von denormalisierten Zahlen (engl. subnormal). In diesem Fall wird keine 1 vor dem Komma angenommen. Außerdem wird mit dem Exponenten −126 für float und −1022 für double gerechnet.

Die Null

Die 0 ist eine denormalisierte Zahl und kann mit beiden Vorzeichen dargestellt werden.

+ 0 :    0 00000000000 0000000000000000000000000000000000000000000000000000
− 0 :    1 00000000000 0000000000000000000000000000000000000000000000000000


#include <stdio.h>
#include <float.h>
int main()
{
      if( NAN != NAN )
      {
          puts( "NAN != NAN" );
      }
      return 0;
}
Listing 2: Beispiel in C zum Vergleich mit NaN.

Not a Number - NaN

Wenn der Exponent nur aus Einsen besteht und die Mantisse mindestens eine Eins enthält so handelt es sich um eine Darstellung für NaN. Das Vorzeichenbit spielt für die Darstellung keine Rolle. Undefinierte Rechenoperationen können NaN als Resultat haben. Es gilt außerdem die Regel, dass jeder Vergleich mit NaN false ergibt − auch ein Vergleich auf sich selbst. Das nebenstehende C-Programm hat die Ausgabe NAN != NAN. In der C-Headerdatei float.h und dem C++-Pendant cfloat ist eine Konstante für NAN definiert.

NaN :    0 11111111111 0000000000000000000000000000000000000000000000000001


Die Darstellung von Unendlich

Divisionen durch Null haben das Ergebnis Unendlich. Je nachdem ob eine positive oder negative Zahl durch Null geteilt wurde ergibt sich ±Unendlich. In der beschriebenen Gleitkommadarstellung besteht der Exponent nur aus Einsen und die Mantisse nur aus Nullen.

 :    0 11111111111 0000000000000000000000000000000000000000000000000000
−  :    1 11111111111 0000000000000000000000000000000000000000000000000000

Binäre Darstellung von beliebigen Gleitkommazahlen...

Zahl:   Vz.   Exp.                         Mantisse
 
 

Die Ungenauigkeit in der Darstellung von Fließkommawerten

Wie im Abschnitt über das Darstellungsproblem von Fließkommawerten schon beschrieben wurde, kann nur eine endliche Anzahl verschiedener Werte überhaupt dargestellt werden - gegenüber einer unendlichen Anzahl nicht darstellbarer Werte. Das heißt auch dass Werte, die aus einer Berechnung entstehen möglicherweise nicht exakt darstellbar sind. Mit hoher Wahrscheinlichkeit sogar. Natürlich hängt das von den verwendeten Operanden einer Operation und von der Operation selbst ab. Die Mehrheit der Rechenergebnisse wird aber nicht exakt darstellbar sein. Dabei muss auch die Konvertierung aus einer Fließkomma­darstellung im Dezimalsystem zum binären System als Operation verstanden werden. Daraus folgt, dass sogar Konstanten in einem Quellcode, die ja in eine binäre Darstellung überführt werden müssen, nicht immer exakt darstellbar sind. Das Codebeispiel in Listing 1 ist ein gutes Beispiel dafür: die dezimal dargestellte Zahl 1.2 ist im binären Exponential­format nicht exakt darstellbar. Sie kann nur näherungsweise repräsentiert werden. Da hilft auch die doppelt Präzision von double nichts. Wie kann also diese Ungenauigkeit erfasst und verstanden werden?

Dazu muss man verstehen, dass die Distanz der darstellbaren Zahlen auf dem Zahlenstrahl zueinander nicht immer gleich ist. Diese Distanz nimmt mit der größe des Betrags zu. Die geringsten Abstände zwischen den darstellbaren Werten sind also direkt um die 0 herum. Nach außen vergrößern sie sich, wie auf der folgenden Abbildung schematisch dargestellt:

                0       +          
· · · · · · · · · · · · · · · · · · · · ·

Zwischen den darstellbaren Zahlen, befinden sich immer unendlich viele nicht-darstellbare. Wie können diese Lücken nun mathematisch gefasst werden und wie get man praktisch mit ihnen um wenn man Fließkommatypen in der Software benötigt?

Das Epsilon als das numerische Maß der Ungenauigkeit

Die Ungenauigkeit in der Darstellung von Fließkommawerten erfasst man mit einem Wert mit dem Namen Epsilon. Das Epsilon ist ein relatives Maß für die Darstellungs­ungenauigkeit. Das Standardepsilon beschreibt die Ungenauigkeit um die Zahl 1.0 herum. Im IEEE754 Standard ist das Epsilon der kleinste Wert für den 1.0 + eps != 1.0 gilt. Da man aber mit Gleitkomma­datentypen Zahlen völlig unterschiedlicher Größen­ordnungen mit der gleichen Anzahl an binären Stellen Anzeigen möchte, skaliert die Ungenauigkeit mit der Größe der dargestellten Zahlen. Die beiden Zahlen 0.0001 und 100000.0 haben beispielweise völlig unterschiedliche Epsilons. Um das Epsilon einer bestimmten Zahl zu ermitteln muss man das Standardepsilon eps mit der Zahl multiplizieren:
eps100000 = 100000 × eps
Für eine beliebige Zahl n gilt also:
epsn = n × eps
Das Darstellung des Ergebnisses n einer einzigen mathematischen Operation kann also in einer Distanz von epsn um das tatsächliche Ergebnis herum liegen. Die Absolute Abweichung der Darstellung vom exakten Ergebnis variiert also mit der Größenordnung von n und wird duch das skalierte Epsilon beschrieben.

Möchte man errechnete Fließkommawerte miteinander vergleichen, muss man den möglichen entstandenen Fehler in Form des Epsilons berücksichtigen. Man errechnet sich also das Epsilon für die entsprechenden Zahlen und macht Bereichsüberprüfungen. Für float, double und long double gibt es in der Standard­bibliothek vordefinierte Standardepsilons.

Die Abweichung denormalisierter Zahlen

Für die denormalisierten Zahlen kann kein Epsilon angegeben werden. Ist das Ergebnis einer Kalkulation denormalisiert, kann keine reguläre Abweichnung in Betracht gezogen werden. Das heißt, dass diese Zahlen zwar darstellbar sind, als Berechnungsergebnisse aber nicht taugen, da man nicht überprüfen kann wie stark sie vom wirklichen mathematischen Ergebnis abweichen. Der denormalisierte Bereich mit Ausnahme der Null kann also nicht als vollgültiger Darstellungsbereich der Fließkommazahlen angesehen werden.

Denormalisierter Bereich
außer der 0
Normalisierter Bereich,
effektiver Darstellungsbereich
float ± { 2−149 bis (1 − 2−23) × 2−126} ± { 2−126 bis (2 − 2−23) × 2127}
double ± { 2−1074 bis (1 − 2−52) × 2−1022} ± { 2−1022 bis (2 − 2−52) × 21023}

Das ULP als das technische Maß der Ungenauigkeit

ULP ist eine Abkürzung von „Unit in the Last Place“ und bezeichnet das am niedrigsten signifikante Bit für die Zahlendarstellung in der Mantisse. Also einfach das letzte Bit. Wenn ein Ergebnis einer Operation nicht darstellbar ist entscheidet es sich in diesem letzten Bit, ob die dargestellte Zahl gewissermaßen rechts oder links des wirklichen Ergebnisses liegt. Ob also das dargestellte Ergebnis größer oder kleiner ist. Natürlich hängt es von der Implemen­tierung der Operation ab, ob der kleinere oder größere Wert als Ergebnis genommen wird. Außerdem kann sich natürlich mehr an dem Bitmuster der Zahl ändern wenn man die nächste darstellbare Zahl erreicht.

Das ULP beschreibt also die Distanz zwischen zwei benachbarten darstellbaren Zahlen. Um das nächste ULP zu berechnen wird ein ein ganzzahliges Inkrement oder Dekrement auf die Mantisse ausgeführt. Ein gefülltes oder geleertes Bitmuster führt natürlich zu einem Kippen und einer Änderung des Exponenten. Die Entfernung zwischen den beiden Zahlen ist aber immer noch ein ULP. Das ULP ist also die kleinstmögliche Differenz zwischen zwei darstellbaren Zahlen. Betrachtet man das Bitmuster der Darstellung von Fließkommawerten ist das ULP ein technisch einfach verständlicher Begriff. Er führt zu genau der relativen Charakter­istik der Darstellung wie das Epsilon im vorange­gangenen Abschnitt. Welche mathematische Distanz ein ULP bedeutet entscheidet sich natürlich am Exponenten. Damit ist die Beschreibung der Ungenauigkeit über das ULP nur ein anderer Weg, um die gleiche Sache, die Abweichung der Genauigkeit in der Fließkokmma­darstellung, zu beschreiben. Das ULP, angewandt auf eine spezifische Zahl, entspricht dem skalierten Epsilon.

Der Zusammenhang von Epsilon und ULP

Das Epsilon ist genau so definiert, dass es bei der Addition auf 1.0 zum nächsten ULP von 1.0 führt. Skaliert man das Epsilon durch die Multiplikation mit einer Zahl im normalisierten Bereich, wird durch die Addition des skalierten Epsilons auf die Zahl das nächste ULP erreicht.
Es gilt also:
n + n × eps = nächstes ULP von n aus.
Die gleiche Regel gilt auch für die Substraktion skalierter Epsilons.

Konsequenzen der Ungenauigkeit von Fließkommawerten in der Codierung

Vergleiche zwischen Fließkommawerten müssen die Ungenauigkeit der Darstellung berücksichtigen. Ein direkter Vergleich zweier Werte ohne die Abweichung mit einzubeziehen würde in sehr vielen Fällen zu falschen Resultaten führen. Listing 1 in diesem Artikel zeigt einen solchen fehlschlagenden Vergleich. Werden zwei berechnete Fließkommawerte verglichen, muss das skalierte Epsilon mit in den Vergleich einfließen, da davon ausgegangen werden muss, dass das Ergebnis nicht exakt dargestellt werden kann und damit um ein Epsilon abweicht.

#include <stdio.h>
#include <float.h> // fuer DBL_EPSILON
#include <math.h>  // fuer fabs()
int main()
{
    double d1 = 1.2;
    double d2 = 0.8;
    d2 += 0.4;
    double e = DBL_EPSILON * fabs(d1);
    puts( ( fabs(d1-d2) <= 2.0*e ) ? "gleich" : "ungleich" );
    return 0;
}
Listing 3: Beispiel in C zum korrekten Vergleich
zweier Gleitkommawerte.

Wenn zwei errechnete Werte verglichen werden, muss die doppelte skalierte Epsilon­distanz angenommen werden, da beide Werte in entgegen­gesetzter Richtung abweichen können. Wenn man die Differenz bildet, kann diese Differenz gegen den doppelten skalierten Epsilonwert verglichen werden. Das Standardepsilon ist als Konstante DBL_EPSILON in der Headerdatei float.h definiert und kann durch eine einfache Multiplikation mit einem der beiden zu vergleichenden Werte skaliert werden. Wenn diese Werte sich wesentlich unterscheiden, fällt das Epsilon nicht ins Gewicht. Wenn sie in direkter Nachbarschaft sind, ist die jeweils andere durch das skalierte Epsilon erreichbar. Wenn beide Werte aus einer Berechnung hervorgegangen sind kann man das Epsilon auch mit dem darstellungsmäßig größeren Wert multiplizieren, um bei doppelter Distanz ein zu kleines Epsilon zu vermeiden. Dann würde die Skalierung so aussehen:

double e = DBL_EPSILON * fmax( fabs(d1), fabs(d2) );

Das Codebeispiel in Listing 3 macht den Vergleich korrekt und prüft ob die beiden verglichenen Werte nahe genug beieinander stehen. Nämlich in einer Distanz von maximal zwei skalierten Epsilons oder eben von zwei ULPs.

Vergleiche von Gleitkommawerten

Nicht nur der Test auf Gleichheit muss das Epsilon berücksichtigen. Auch der Größer- und Kleinervergleich muss abgesichert werden. Anstatt zwei Werte einfach mit dem Vergleichsoperator zu vergleichen, muss vorher das Epsilon in richtiger Weise addiert oder Subtrahiert werden. Falsch wäre if( d1 < d2 ).. , da die beiden Werte so nahe beieinander liegen können, dass sie als gleich betrachtet werden müssten. Der Operator jedoch prüft nur die tatsächlich dargestellten Zahlen ohne die Ungenauigkeit der Darstellung des Ergebnisses des vorangegangenen Rechenwegs in Betracht zu ziehen. Das Epsilon muss also in den Vergleich mit einbezogen werden: if( ( d1 + e ) < d2 ).. oder eben if( ( d1 + e ) < ( d2 - e ) ).. , wenn beide Werte aus Berechnungen hervor­gegangen sind und damit die doppelte Epsilon­distanz berücksichtigt werden muss.

Vorsicht Falle! Die Operatoren „<=“ und „>=

#include <stdio.h>
#include <float.h> // fuer DBL_EPSILON
#include <math.h>  // fuer fabs()
int main()
{
    double d1 = 1.2;
    double d2 = 0.8;
    d2 += 0.4;
    double e = DBL_EPSILON * fabs(d1);
    puts( ( fabs(d1-d2) <= 2.0*e ) ? "gleich" : "ungleich" );
    puts( ( (d1+e) < (d2-e) ) ? "kleiner" : "nicht kleiner" );
    // Fehler:
    puts( ( (d1+e) <= (d2-e) ) ? "kleiner-gleich"
                               : "nicht kleiner-gleich" );
    return 0;
}
Listing 4: Beispiel in C mit fehlerhafter Anwendung des Operators „<=“.

Der nebenstehende Code enthält einen Fehler. Er produziert die Ausgabe:

gleich
nicht kleiner
nicht kleiner-gleich


Das ist natürlich nicht konsistent, denn wenn Werte gleich sind, dann sind sie formal auch kleiner-gleich. Die erwartete Ausgabe wäre also:

gleich
nicht kleiner
kleiner-gleich


Der Kleiner-Gleich-Operator „<=“ kann nicht naiv angewendet werden, als wäre er ein Kleiner-Operator „<“. Bei einem Kleiner-Vergleich wird das Epsilon so angewendet, dass der getestete angenommern kleinere Wert vergrößert wird. Er könnte aber bereits in seiner Darstellung größer sein. Die Addierung des Epsilong kann den Wert also aus dem Bereich bringen, in dem Gleichheit angenommen werden kann. Wenn also ein logischer Kleiner-Gleich-Vergleich gebraucht wird muss er aus einem logischen Kleiner-Vergleich und der logischen Überprüfung auf Gleichheit zusammen­gesetzt werden. Eine mögliche korrekte Lösung wäre also die folgende:

Download des korrigierten
Codes aus Listing 4.

puts( ( (d1+e) < (d2-e) || fabs(d1-d2) <= 2.0*e )
      ? "kleiner-gleich" : "nicht kleiner-gleich" );

Die Bedingung des Vergleichs besteht aus zwei veroderten Teilbedingungen. Also entweder der Wert d1 ist kleiner als d2 oder die Werte sind gleich. Hier zeigt sich ganz besonders, dass logische Vergleiche nicht mit den gleich­lautenden Zahlenwert­vergleichen verwechselt werden dürfen.

Dieser Fehler in der Anwendung der Operatoren „<=“ und „>=“ auf Gleitkommawerte wird meiner Erfahrung nach oft begangen. Selbst dann, wenn skalierte Epsilons in Vergleichen Anwendung finden. Das liegt auch sicherlich daran, dass die Materie nicht unmittelbar intuitiv ist und dass diese beiden Operatoren eine gesonderte Aufmerksamkeit erfordern. Sie dürfen nicht verwendet werden, um einen logischen Vergleich zwischen Fließkommawerten zu formulieren. Gleichwohl können sie angewendet werden, um die Zahlen­darstellung zu vergleichen, wie es beispielsweise im logischen Vergleich auf Gleichheit geschieht.

Die Anwendung von Gleitkommawerten und Hilfestellung aus der Bibliothek

Wer das volle Potential der Fließkommawerte in digitalen Systemen nutzen möchte tut gut daran, ein Studium der höheren Mathematik mit Vertiefungs­richtung Numerik absolviert zu haben. Dann dürfte es vielleicht noch angesagt sein, den eigenen Code gegenüber jeder Änderung von fremder Hand standhaft zu verteidigen, wenn diese Hand nicht im Auftrag einer mindestens ebenso spezifisch akademisch gestählten und verfeinerten Persönlichkeit handelt.
Da ich mich selbst nicht zu diesem Personenkreis zählen kann, besteht meine vorrangige Strategie der korrekten Verwendung von Gleitkomma­arithmetik in deren Vermeidung. Ich mache es nur wenn ich gezwungen werde!
Solange man etwas in eine ganzzahlige Darstellung bringen kann, sollte man es tun. Nur wenn die Rahmen­bedingungen die Nutzung von Gleitkomma­darstellung notwendig macht, wendet man sie an. Und dann mit einem möglichst tiefen Verständnis. Kommen wir nun zu einem sehr üblichen, meiner Erfahrung nach dem häufigsten Anwendungs­fall der Fließkomma­arithmetik:

Die Überprüfung von Messwerten anhand von Schwellenwerten

Was passiert, wenn Vergleiche anhand von Schwellenwerten durchgeführt werden? Im Allgemeinen kann gesagt werden, dass Werte auf eine bestimmte Maßeinheit umgerechnet werden, um sie schließlich vergleichen zu können. Um Werte auf eine Einheit zu normalisieren müssen sie gestreckt oder gestaucht und verschoben werden. Die Streckung oder Stauchung erfolgt durch eine Multiplikation oder Division, die Verschiebung durch eine Addition oder Substraktion. In den meisten Fällen der Messwert­bearbeitung kann also von der Erhebung des Wertes bis zum Test gegenüber dem Schwellenwert von zwei Rechen­operationen ausgegangen werden: von einer Punkt- und von einer Strichoperation. Glücklicherweise und höchst­wahr­scheinlich nicht ganz zufällig kennt der IEEE 754 eine solche Doppel­operation, bei der nur eine einfache Epsilon­abweichung an Ungenauigkeit anfällt. In der C-Bibliothek wird ab C99 diese Doppel­operation in Form von Funktionen für die drei Fließkomma­datentypen float, double und long double angeboten:

float  fmaf( float x, float y, float z );
double  fma( double x, double y, double z );
long double  fmal( long double x, long double y, long double z );

Diese Funktionen berechnen (x × y) + z mit nur einer einzigen maximalen möglichen ULP Abweichung im Ergebnis. Sie garantieren, dass keine Fehler­fort­schreibung zwischen den beiden mathematischen Operationen stattfindet. Damit steht eine Lösung für ein weitverbreitetes Standard­problem zur Verfügung. Die Normierung von Messwerten kann in einem Arbeitsschritt erfolgen der aus zwei arithmetischen Operationen besteht aber nur eine Standard­ungenauigkeit verursacht. Die oben beschriebenen Verfahren der Gleitkomma­vergleiche können also direkt zur Anwendung gebracht werden.

Weitere Bibliotheksfunktionen

Im Abschnitt über die ULPs wurde die Distanz zwischen den darstellbaren Zahlen beschrieben. Mathematisch werden diese Distanzen mit skalieren Epsilons ausgedrückt, technisch sind es die ULPs, die Units in the Last Place. Die Distanz zwischen zwei darstellbaren Zahlen ist also ein ULP oder eben ein skaliertes Epsilon. Das bedeutet, dass es für den Umgang mit der Darstellungs­ungenauigkeit von Gleitkomma­werten neben der mathematischen Lösung über die Berechnung des Epsilons auch eine technische Lösung existieren kann, die mit den ULPs operiert.
Dazu gibt es ab dem C-Standard C99 eine kleine Sammlung von Funktionen in der Headerdatei <math.h>. Für C++ steht diese Funktions­sammlung ab dem Standard C++11 in der Headerdatei <math> zur Verfügung. Beginnen wir mit der folgenden Funktion, die stell­vertretend für die ganze Funktions­gruppe gesehen werden kann:

double  nextafter( double f /* from */, double t /* to */ );

Diese Funktion liefert die nächste darstellbare Zahl in der Nachbarschaft von f in die Richtung des Wertes von t.

#include <math.h>
#include <stdio.h>

int main()
{
    double a = nextafter( 0.0,  1.0 );
    double b = nextafter( 0.0, -1.0 );
    printf( "Die kleinste positive Zahl des Typs double: %E\n", a );
    printf( "Die kleinste negative Zahl des Typs double: %E\n", b );
    return 0;
}
Listing 5: C-Beispiel zur Anwendung der Funktion nextafter().

Das Beispiel in Listing 5 zeigt die An­wendung der Funk­tion. Sie liefert aus­gehend vom Wert des ersten Para­meters je nachdem, ob der zweite Para­meter kleiner oder größer als der erste Para­meter ist, den nächsten dar­stell­baren kleineren oder eben den nächsten dar­stell­baren größeren Wert.

Die anderen Funktionen in der kleinen Funktions­familie machen eigenlich genau das gleiche, sie arbeiten nur mit anderen Typen. Also auch mit float und long double. Es gibt die Funktionen mit den Namen nextafter..() und solche mit den Namen nexttoward..(). An der Stelle der beiden Punkte kann entweder gar nichts, ein „f“ für float oder ein „l“ für long double stehen.
Die Funktionen mit den Namen nexttoward..() nehmen als zweiten Parameter immer ein long double. Die Funktionen sind unter anderem in der C++ Referenz hinter diesem Link beschrieben. Sie sind auch anderstwo recht gut dokumentiert und wenn man einen der Namen kennt, findet man leicht ausreichende Dokumentation.

Setzt man diese Funktionen nun ein, um Gleitkomma­werte korrekt miteinander zu vergleichen, muss man natürlich die mögliche Anzahl der voraus­gegangenen Berechnungs­schritte mit in Betracht ziehen, genau so, wie man es auch mit einer auf das Eposilon basierten Bereichsüberprüfung machen würde. Wenn beide Werte errechnet sind, muss die doppelte ULP Distanz herangezogen werden.


Links zum Thema Fließkommazahlen


Zuletzt geändert am 13.01.2024