Vergleich von Gleitkommazahlen – knifflig, aber machbar!

Gleitkomma-Mathematik ist mit recht diffizilen und subtilen Problemen belastet. Der Vergleich von Werten macht da keine Ausnahme. In diesem Artikel diskutieren wir häufige Fallstricke, untersuchen mögliche Lösungen und versuchen, „Boost“-Probleme zu überlisten.

Wenn Sie eine nicht ganzzahlige Zahl in einer der etablierten Programmiersprachen darstellen wollen, werden Sie letzten Endes wahrscheinlich Gleitkomma-Zahlen nach IEEE 754 verwenden. Nach ihrer Normung 1984 wurden sie einfach allgegenwärtig. Fast alle modernen CPUs – und viele Mikroprozessoren – haben für sie eine spezielle Hardware, die sogenannten Floating-Point Units (FPUs, Gleitkomma-Einheiten).

Jede Variable des Typs float besteht aus einem Vorzeichenbit, einigen Bits, die den Exponenten repräsentieren und Bits, die einen Bruchteil, auch Mantisse genannt, darstellen. In den meisten Fällen kann der Wert einer Gleitkomma-Zahl durch diese Gleichung dargestellt werden.

Hierbei stellt s unser Vorzeichenbit dar, m einen Bruchteil, der durch die Bits der Mantisse repräsentiert wird, e eine vorzeichenlose Ganzzahl, die durch die Bits des Exponenten dargestellt wird und c ist der halbe Maximalwert von e, das heißt 127 für eine 32-Bit Gleitkomma-Zahl und 1023 für eine 64-Bit Gleitkomma-Zahl.

Es gibt aber auch einige Spezialfälle. Sind zum Beispiel alle Bits des Exponenten Null, ändert sich die Formel entsprechend. Bemerkenswert ist der Entfall der „1“ vor der Mantisse. Dies erlaubt es uns, kleine Werte nahe Null zu speichern, sogenannte normalisierte oder denormalisierte Werte. Und wenn alle Bits im Exponenten Einsen sind, repräsentieren bestimme Werte der Mantisse +∞, - ∞ und keine Zahl (NaN / Not a Number), das Ergebnis einer undefinierten oder nicht darstellbaren Rechenoperation wie die Division durch Null.

Hier sind zwei Beobachtungen festzuhalten, die sich in Kürze als nützlich herausstellen.

1. Variablen vom Typ float können keine beliebigen realen Zahlen speichern, ganz zu schweigen von beliebigen rationalen Zahlen. Sie können ausschließlich Zahlen speichern, die durch die oben gezeigten Gleichungen dargestellt werden können. Wenn ich zum Beispiel die folgende Variable deklariere:

wird f als 0.100000001490116119384765625 dargestellt, was dem an 0,1 nächstliegenden 32-Bit Gleitkomma-Wert entspricht.

2. Da die Gleichungen exponentiell sind, nimmt der Abstand auf der Zahlengerade zwischen benachbarten Werten zu – und zwar exponentiell, je weiter man sich von der Null entfernt. Der Abstand zwischen 1,0 und dem nächsten möglichen Wert ist ungefähr 1,19 x 10^-7, aber der Abstand zwischen zwei benachbarten Gleitkomma-Zahlen in der Nähe von 6,022 x 10²³ ist ungefähr 3,6 x 10¹⁶. Das wird sich als eine unserer größten Herausforderung herausstellen: Wenn wir Gleitkomma-Zahlen vergleichen wollen, dann wollen wir Eingaben in der Nähe von Null genauso handhaben, wie Zahlen in der Nähe der Avogadro-Konstante.

Was ist Gleichheit?

Da das Ergebnis jeder Gleitkomma-Operation auf den nächstmöglichen Wert gerundet werden muss, verhält sich die Mathematik nicht so wie bei realen Zahlen. Abhängig von Ihrer Hardware, dem Compiler und den Compiler-Einstellungen, kann die Operation 0,1 x 10 ein anderes Ergebnis produzieren als ∑(n=1 bis 10) 0,1 (Auf meiner Maschine ergibt das Erstere 1,0 und das Letztere 1,000000119209290.).

Immer wenn wir berechnete Werte miteinander vergleichen, sollten wir einen gewissen Spielraum berücksichtigen. Ein Vergleich der genauen Werte mit == bringt es nicht. Stattdessen sollten wir zwei unterschiedliche Zahlen a und b als gleich ansehen, wenn für ein ausreichend kleines ε gilt: │a – b│ ≤ ε. Doch wie es der Zufall will, enthält die C-Standardbibliothek eine Definition FLT_EPSILON. Probieren wir es aus:

Eigentlich würden wir jetzt darauf hoffen, fertig zu sein, aber damit würden wir falsch liegen. Ein Blick in den Standard der Sprache zeigt, dass FLT_EPSILON gleich dem Abstand zwischen 1,0 und dem nächstliegenden Wert ist.

Aber, wie wir vorher festgestellt haben, sind Gleitkomma-Zahlen nicht äquidistant! Für Werte kleiner 1, wird FLT_EPSILON schnell zu groß, um nützlich zu sein. Für Werte größer 2, ist FLT_EPSILON kleiner als der Abstand zweier benachbarter Werte, so dass fabs(a - b) <= FLT_EPSILON immer false ergibt.

Was passiert, wenn wir zur Lösung dieser Probleme ε proportional zu unseren Eingangswerten skalieren?

Das funktioniert besser als unsere ursprüngliche Version, aber es ist nicht sofort offensichtlich, welche Werte von maxRelativeDiff wir für verschiedene Fälle verwenden wollen. Die Tatsache, dass wir maxRelativeDiff mit beliebigen Eingangswerten skalieren, kann bedeuten, dass wir demselben Rundungsmechanismus zum Opfer fallen, über den wir uns anfangs Sorgen gemacht hatten.

Was ist mit Boost?

Boost, die populäre Sammlung von C++ Bibliotheken, stellt Funktionen für ähnliche Zwecke zur Verfügung (Siehe „Floating Point Comparison“ im Abschnitt über die Gleitkomma-Hilfsprogramme im mathematischen Toolkit von Boost.). Entfernt man die Templates und die Behandlung von Randfällen wie ±∞ und NaNs, sehen sie ungefähr so aus:

Unglücklicherweise scheinen diese Funktionen unser Problem nicht zu lösen. Da die Division in relative_difference oftmals sehr kleine Ergebnisse liefert – der Wert für relative_difference beispielsweise zwischen 42 und dem nächstmöglichen Gleitkomma-Wert ist ungefähr 9,08 x 10^-8 – stellt sich die Frage, was wohl ein guter Schwellenwert wäre. Durch die Division mit FLT_EPSILON versucht epsilon_difference einen einfacheren Wert zur Verfügung zu stellen. Aber wir haben die Gefahren von FLT_EPSILON gerade gesehen! Dieses Schema wird umso fragwürdiger, je weiter sich die Eingangswerte von Eins entfernen.

Was ist mit ULPs?

Es wäre schön, wenn Vergleiche mit etwas konkreterem als willkürlich gewählten Schwellenwerten durchgeführt werden könnten. Idealerweise würden wir gerne die Anzahl möglicher Gleitkomma-Werte – manchmal als „units of least precision (ULPs)“ bezeichnet – zwischen den Eingabewerten wissen.

Wenn ich irgendeinen Wert a habe und ein anderer Wert b ist nur zwei oder drei ULPs entfernt, dann können wir beide wahrscheinlich als gleich betrachten, in der Annahme, dass es sich um Rundungsfehler handelt. Und wichtiger: Dies gilt unabhängig vom Abstand zwischen a und b auf der Zahlengerade.

Boost stellt eine Funktion mit dem Namen float_distance zur Verfügung, um den Abstand zwischen Werten in ULPs zu bekommen, aber sie ist ungefähr eine Größenordnung langsamer als die Ansätze, die wir bisher diskutiert haben. Mit ein wenig Bitmanipulation können wir das besser machen.

Betrachten wir ein beliebige Gleitkomma-Zahl x, bei der jedes Bit in der Mantisse Eins ist. x + 1ULP muss dann den nächstgrößeren Exponenten verwenden und alle Bits in der Mantisse müssen dann Null sein. Schauen wir uns als Beispiel 1,99999988 und 2 an:

Die Eigenschaft gilt für denormalisierte Zahlen, obwohl sie eine unterschiedliche Wertgleichung besitzen. Schauen wir uns die größte denormalisierte und die kleinste normalisierte Zahl an:

Beachten Sie interessante logische Folge: Benachbarte Gleitkomma-Zahlen (mit dem gleichen Vorzeichen) besitzen benachbarte ganzzahlige Werte, wenn man sie als solche neu interpretiert. Diese Neuinterpretation wird manchmal als „type punning“ (Umgehung der Datentypen) bezeichnet und wir können es benutzen, um den Abstand zwischen Werten in ULPs zu berechnen.

Was ist mit ULPs?

Traditionellerweise benutzte man bisher in C und C++ den union-Trick:

Dies funktioniert in C immer noch, kann in C++ aber mit den strikten Aliasing-Regeln kollidieren (Siehe z.B. http://stackoverflow.com/q/11639947 und http://stackoverflow.com/q/17789928). Eine bessere Methode ist es, memcpy zu verwenden. Geht man vom normalen Gebrauch dieser Funktion aus, könnte man annehmen, dass dies weniger effizient ist, aber sehen wir uns den Code an:

Dieses Beispiel compiliert zu einer einzigen Instruktion, die den Wert von einem Gleitzahl-Register in ein Ganzzahl-Register kopiert. Das ist genau das, was wir wollen.

Nachdem wir dieses Problem gelöst haben, wird die Berechnung der ULPs zwischen Werten ziemlich unkompliziert. Wir haben nun einen recht portierbaren Code gefunden; er setzt nur voraus, dass die Plattform 32-Bit Integer unterstützt und dass die Gleitkomma-Zahlen gemäß IEEE 754 gespeichert werden. (Anmerkung: Der Typ float ist gemäß dem C++ Standard implementationsabhängig und beachtet nicht notwendigerweise IEEE 754. Siehe auch http://stackoverflow.com/a/24157568). Dies ist eventuell der Grund für Boost’s derzeitige Implementation von float_distance.)

Wir vermeiden den Vergleich von Werten mit unterschiedlichen Vorzeichen aus einer Reihe von Gründen:

• 1. ULPs sind, wie wir noch sehen werden, nicht das richtige Werkzeug um Werte in der Nähe von oder jenseits von Null zu vergleichen.

• 2. Nahezu alle modernen CPU benutzen Zweierkomplement-Arithmetik, während Gleitkommazahlen Vorzeichen und Größe benutzen. Die Umwandlung von einem Format in das andere benötigt zusätzliche Arbeit, um sinnvoll Zahlen mit unterschiedlichen Vorzeichen addieren oder subtrahieren zu können. Aus demselben Grund muss das Vorzeichen unseres Ergebnisses nicht das sein, das wir erwarten. Daher nehmen wir den Absolutwert. Uns interessieren nur die Abstände zwischen unseren beiden Eingabewerten.

• 3. Wenn die Subtraktion einen Über- oder Unterlauf erzeugt, bekommen wir ein undefiniertes Verhalten bei vorzeichenbehafteten Integerzahlen und bei modularer Arithmetik mit vorzeichenlosen Zahlen. Nichts davon ist hier erstrebenswert.

Wir berechnen den Absolutwert aus zwei Gründen selbst, statt std::abs zu verwenden. Erstens akzeptieren die Ganzzahl-Versionen von std::abs nur Datentypen, deren Größe plattform-spezifisch ist, wie int, long oder long long. Wir wollen Annahmen über die implizite Umwandlung zwischen diesen Typen und int32_t vermeiden (auch wenn das zugegebenermaßen etwas an Paranoia grenzt: int32_t ist einer der Datentypen, die es auf fast jeder relevanten Plattform gibt.)

Der zweite Grund ist ein sonderbarer Fallstrick, der mit der Anordnung von std::abs Überladungen in der C++ Standardbibliothek zu tun hat. Wenn Sie <cmath> einbinden, aber nicht <cstdlib>, stehen nur die Gleitkomma-Versionen von std::abs zur Verfügung. Mehrere Toolchains, die ich getestet habe, wandeln einen int32_t Wert in einen double um, auch wenn die Zielplattform nur über eine 32-Bit FPU verfügt und daher double in Integer-Registern emulieren muss (wie man vermuten kann, mit fürchterlichen Auswirkungen auf die Laufzeit). Die Verwendung von Compiler-Warnungen wie -Wconversion kann uns darauf aufmerksam machen, oder wir können alle diese Fallen generell vermeiden, indem wir den Absolutwert direkt berechnen. Jedenfalls ist das ein triviales Detail.

Es gibt keine Wunderwaffe

Relative Epsilons – einschließlich derer, die auf ULPs basieren – machen um Null herum keinen Sinn. Die exponentielle Natur von Gleitkomma-Zahlen bedeutet, dass sich dort viel mehr Werte auf der Zahlengerade versammeln, also irgendwo sonst. Obwohl es im Kontext vieler Berechnungen nur ein verhältnismäßig kleiner Wert ist, ist 0,1 über eine Milliarde ULPs von der Null entfernt!

Konsequenterweise sind also feste Epsilons möglicherweise die beste Wahl, wenn sie kleine Ergebnisse erwarten. Welches ε sie verwenden, hängt vollständig von den durchgeführten Berechnungen ab.

Mit diesem Wissen bewaffnet, besteht nun möglicherweise die Versuchung, eine Funktion wie die folgende zu schreiben, die all dem ein Ende setzt:

Aber ein aussagekräftiger Einsatz ist ohne die gerade diskutierte Theorie schwierig.

Ein kurzer Exkurs: Weitere ULP-basierte Funktionen

Dieselbe Technik können wir zum Schreiben weiterer nützlicher Funktionen benutzen. Zum Beispiel für eine Funktion, die eine Gleitkomma-Zahl um eine beliebige Zahl von ULPs inkrementiert. Boost bietet eine ähnlich Familie von Funktionen an (float_next, float_advance, usw.), aber wie bei der float_distance Funktion geht das auf Kosten der Leistung um Type Punning zu vermeiden.

Man könnte jetzt darauf hoffen, dass wir einfach unsere ULPs berechnen, unsere Addition durchführen und dann das Ergebnis mittels memcpy zurückgeben, wie beispielsweise diese Funktion:

Diese etwas naive Lösung funktioniert für positive Werte, aber auf der meisten Hardware bedeutet das Inkrementieren einer negativen Gleitkomma-Zahl um eine positive Anzahl von ULPs, dass wir uns von der Null wegbewegen. Das ist wahrscheinlich nicht das, was wir wollen.

Wir haben vorher erwähnt, dass Gleitkomma-Zahlen dem Schema „Vorzeichen – Exponent“ folgen, wogegen die meisten CPUs das Zweierkomplement benutzen.

Um also mit negativen Zahlen arbeiten zu können, müssen wir die Werte in das systemeigene Integerformat der CPU* umwandeln.

(*Auf esoterischer Hardware kann das systemeigene Format ebenfalls dem Vorzeichen – Exponent Schema folgen. In diesen Fällen vertrauen wir darauf, dass der Compiler unsere nutzlose Arbeit weglässt.)

Nach der Operation mit den ULPs müssen wir wieder zurückwandeln:

Da diese Funktionen nunmehr definiert sind, können wir uns wieder unserem Ziel zuwenden:

Was haben wir gelernt?

Wenn Gleitkomma-Zahlen verglichen werden, dürfen wir nicht vergessen:

• FLT_EPSILON ist nicht das Gleitkomma-Epsilon, das wir benötigen, außer in den Bereichen [-2,-1] und [1,2]. Der Abstand zwischen angrenzenden Werten hängt von den eigentlichen Werten ab.

• Beim Vergleich mit einigen bekannten Werten – speziell Null oder Werten in der Nähe von Null – verwenden Sie ein festes ε, das für Ihre Berechnungen Sinn ergibt.

• Beim Vergleich von Werten, die nicht Null sind, ist ein auf ULPs basierender Vergleich wahrscheinlich die beste Wahl.

• Könnten sich die Zahlen irgendwo auf der Zahlengeraden befinden, ist eine hybride Form der beiden notwendig. Wählen Sie Ihr Epsilon sorgfältig auf der Basis der erwarteten Ergebnisse.

Danksagungen

Vieles hier wurde von der fantastischen Abhandlung des Themas durch Bruce Dawson in seinem Blog Random ASCII übernommen. Vielen Dank auch an meine Kollegen Evan Thomson und Matt Drees für ihre Beiträge.

* Der Autor: Matt Kline ist Software Engineer für Fluke Networks und Betreiber des Software-Entwicklungs-Blogs bitbashing.io. Der Beitrag "Comparing Floating-Point Numbers Is Tricky" wurde lizensiert nach einer Creative Commons Attribution-ShareAlike 4.0 International License. (CC BY-SA 4.0). Übersetzung Richard Oed.

(ID:44970621)