15. Gleitkommaarithmetik: Probleme und Einschränkungen
******************************************************

Gleitkommazahlen werden in der Computerhardware als Brüche zur Basis 2
(binär) dargestellt.  Beispielsweise hat der **dezimale** Bruch
"0.625" den Wert 6/10 + 2/100 + 5/1000, und ebenso hat der **binäre**
Bruch "0.101" den Wert 1/2 + 0/4 + 1/8. Diese beiden Brüche haben
identische Werte; der einzige wirkliche Unterschied besteht darin,
dass der erste in der Bruchschreibweise zur Basis 10 und der zweite in
der zur Basis 2 geschrieben ist.

Leider lassen sich die meisten Dezimalbrüche nicht exakt als
Binärbrüche darstellen. Dies hat zur Folge, dass die von Ihnen
eingegebenen Dezimal-Gleitkommazahlen im Allgemeinen nur annähernd
durch die tatsächlich im Rechner gespeicherten Binär-Gleitkommazahlen
wiedergegeben werden.

Das Problem lässt sich zunächst im Dezimalsystem leichter verstehen.
Betrachten wir den Bruch 1/3. Man kann ihn als Bruch im Dezimalsystem
näherungsweise darstellen:

   0.3

oder besser noch:

   0.33

oder besser noch:

   0.333

und so weiter. Ganz gleich, wie viele Stellen man notieren möchte, das
Ergebnis wird niemals genau 1/3 sein, sondern eine immer bessere
Annäherung an 1/3 darstellen.

Ebenso lässt sich der Dezimalwert 0,1 – ganz gleich, wie viele Ziffern
im Zweier-System man verwenden möchte – nicht exakt als Bruch im
Zweier-System darstellen. Im Zweier-System ist 1/10 der sich unendlich
wiederholende Bruch

   0.0001100110011001100110011001100110011001100110011...

Hält man bei einer beliebigen endlichen Anzahl von Bits an, erhält man
eine Näherung. Auf den meisten heutigen Rechnern werden
Gleitkommazahlen mithilfe eines binären Bruchs approximiert, wobei der
Zähler aus den ersten 53 Bits (beginnend mit dem höchstwertigen Bit)
besteht und der Nenner eine Zweierpotenz ist.  Im Fall von 1/10 lautet
der Binärbruch "3602879701896397 / 2 ** 55", was nahe am tatsächlichen
Wert von 1/10 liegt, diesem aber nicht exakt entspricht.

Vielen Benutzern ist diese Annäherung aufgrund der Art und Weise, wie
Werte angezeigt werden, nicht bewusst. Python gibt lediglich eine
dezimale Annäherung an den tatsächlichen Dezimalwert der vom Rechner
gespeicherten binären Annäherung aus. Würde Python auf den meisten
Rechnern den tatsächlichen Dezimalwert der für 0,1 gespeicherten
binären Annäherung ausgeben, müsste es Folgendes anzeigen:

   >>> 0.1
   0.1000000000000000055511151231257827021181583404541015625

Das sind mehr Stellen, als die meisten Menschen für sinnvoll halten.
Deshalb sorgt Python dafür, dass die Anzahl der Stellen überschaubar
bleibt, indem stattdessen ein gerundeter Wert angezeigt wird:

   >>> 1 / 10
   0.1

Denke daran: Auch wenn das ausgegebene Ergebnis genau wie der Wert
1/10 aussieht, ist der tatsächlich gespeicherte Wert der
nächstgelegene darstellbare Binärbruch.

Interessanterweise gibt es viele verschiedene Dezimalzahlen, die
denselben nächstliegenden binären Näherungsbruch haben. Beispielsweise
werden die Zahlen "0.1", "0.10000000000000001" und
"0.1000000000000000055511151231257827021181583404541015625" alle durch
"3602879701896397 / 2 ** 55" approximiert. Da alle diese Dezimalwerte
dieselbe Approximation haben, könnte jeder einzelne von ihnen
angezeigt werden, ohne dass die Invariante "eval(repr(x)) == x"
verloren geht.

Früher wählten die Python-Eingabeaufforderung und die integrierte
Funktion  "repr()"  die Zahl mit 17 signifikanten Stellen aus:
"0.10000000000000001". Ab Python 3.1 ist Python (auf den meisten
Systemen) nun in der Lage, die kürzeste dieser Zahlen auszuwählen und
einfach "0.1" anzuzeigen.

Beachte, dass dies in der Natur der binären Gleitkommazahlen liegt: Es
handelt sich weder um einen Fehler in Python noch um einen Fehler in
Ihrem Code.  Das gleiche Phänomen tritt in allen Sprachen auf, die die
Gleitkommaarithmetik Ihrer Hardware unterstützen (auch wenn manche
Sprachen den Unterschied standardmäßig oder in nicht allen Ausgabemodi
möglicherweise nicht *anzeigen*).

Um die Ausgabe übersichtlicher zu gestalten, kannst du die
Zeichenfolgenformatierung verwenden, um eine begrenzte Anzahl
signifikanter Stellen auszugeben:

   >>> format(math.pi, '.12g')  # gibt 12 signifikante Stellen an
   '3.14159265359'

   >>> format(math.pi, '.2f')   # gibt 2 Stellen nach dem Komma an
   '3.14'

   >>> repr(math.pi)
   '3.141592653589793'

Man muss sich bewusst machen, dass es sich hierbei im wahrsten Sinne
des Wortes um eine Illusion handelt: Man rundet lediglich die
*Anzeige* des tatsächlichen Maschinenwerts.

Eine Täuschung kann eine weitere nach sich ziehen. Da beispielsweise
0.1 nicht genau 1/10 ist, ergibt die Summe von drei Werten von 0.1
möglicherweise auch nicht genau 0.3:

   >>> 0.1 + 0.1 + 0.1 == 0.3
   Falsch

Da sich die Zahl 0.1 dem exakten Wert von 1/10 nicht weiter annähern
kann und 0.3 sich dem exakten Wert von 3/10 nicht weiter annähern
kann, hilft auch eine Vorrundung mit der Funktion  "round()"  nicht
weiter:

   >>> round(0.1, 1) + round(0.1, 1) + round(0.1, 1) == round(0.3, 1)
   False

Auch wenn die Zahlen nicht näher an ihre beabsichtigten exakten Werte
herangeführt werden können, kann die Funktion  "math.isclose()"  beim
Vergleich ungenauer Werte nützlich sein:

   >>> math.isclose(0.1 + 0.1 + 0.1, 0.3)
   True

Alternativ kann die Funktion  "round()"  verwendet werden, um grobe
Näherungswerte zu vergleichen:

   >>> round(math.pi, ndigits=2) == round(22 / 7, ndigits=2)
   True

Die binäre Gleitkommaarithmetik birgt viele solcher Überraschungen.
Das Problem mit "0,1" wird weiter unten im Abschnitt
"Darstellungsfehler" ausführlich erläutert. Unter Beispiele für
Gleitkomma-Probleme findest du eine anschauliche Zusammenfassung der
Funktionsweise der binären Gleitkommaarithmetik und der Arten von
Problemen, die in der Praxis häufig auftreten.  Siehe auch Die Tücken
der Gleitkommazahlen für eine umfassendere Darstellung weiterer
häufiger Überraschungen.

Wie es gegen Ende heißt: "Es gibt keine einfachen Antworten." Dennoch
solltest du gegenüber Gleitkommazahlen nicht übermäßig misstrauisch
sein! Die Fehler bei Gleitkommaoperationen in Python stammen von der
Gleitkomma-Hardware und liegen auf den meisten Rechnern bei höchstens
1 Teil pro 2**53 pro Operation.  Das ist für die meisten Aufgaben mehr
als ausreichend, aber man muss bedenken, dass es sich nicht um
Dezimalrechnen handelt und dass bei jeder Gleitkommaoperation ein
neuer Rundungsfehler auftreten kann.

Es gibt zwar Ausnahmefälle, doch bei den meisten alltäglichen
Anwendungen der Gleitkommaarithmetik erhältst du letztendlich das
erwartete Ergebnis, wenn du die Anzeige deiner Endergebnisse einfach
auf die gewünschte Anzahl von Dezimalstellen rundest.  "str()"  reicht
in der Regel aus; für eine feinere Steuerung findest du die
Formatbezeichner der Methode  "str.format()"  unter Syntax für
Formatzeichenketten.

Für Anwendungsfälle, die eine exakte Dezimaldarstellung erfordern,
solltest du das Modul  "decimal"  verwenden, das eine
Dezimalarithmetik implementiert, die für Buchhaltungsanwendungen und
Anwendungen mit hoher Genauigkeit geeignet ist.

Eine weitere Form der exakten Arithmetik wird vom Modul  "fractions"
unterstützt, das eine auf rationalen Zahlen basierende Arithmetik
implementiert (sodass Zahlen wie 1/3 exakt dargestellt werden können).

Wenn du häufig mit Gleitkommaoperationen arbeitest, solltest du das
NumPy-Paket sowie viele andere Pakete für mathematische und
statistische Operationen ansehen, die vom SciPy-Projekt bereitgestellt
werden. Siehe <https://scipy.org>.

Python bietet Werkzeuge, die in den seltenen Fällen hilfreich sein
können, in denen man den genauen Wert einer Gleitkommazahl wirklich
*wissen* möchte. Die Methode  "float.as_integer_ratio()"  gibt den
Wert einer Gleitkommazahl als Bruch an:

   >>> x = 3.14159
   >>> x.as_integer_ratio()
   (3537115888337719, 1125899906842624)

Da das Verhältnis exakt ist, lässt sich damit der ursprüngliche Wert
verlustfrei wiederherstellen:

   >>> x == 3537115888337719 / 1125899906842624
   True

Die Methode  "float.hex()"  gibt einen Float-Wert im Hexadezimalformat
(Basis 16) aus und liefert wiederum den exakten Wert, der auf Ihrem
Computer gespeichert ist:

   >>> x.hex()
   '0x1.921f9f01b866ep+1'

Anhand dieser exakten Hexadezimal-Darstellung lässt sich der Float-
Wert exakt rekonstruieren:

   >>> x == float.fromhex('0x1.921f9f01b866ep+1')
   True

Da die Darstellung exakt ist, eignet sie sich gut für die zuverlässige
Übertragung von Werten zwischen verschiedenen Python-Versionen
(Plattformunabhängigkeit) und den Datenaustausch mit anderen Sprachen,
die dasselbe Format unterstützen (wie beispielsweise Java und C99).

Ein weiteres nützliches Werkzeug ist die Funktion  "sum()" , die dazu
beiträgt, Präzisionsverluste bei der Summierung zu minimieren. Sie
nutzt erweiterte Präzision für die Zwischenrunden, während Werte zur
laufenden Summe hinzugefügt werden. Dies kann sich auf die
Gesamtgenauigkeit auswirken, sodass sich die Fehler nicht so weit
summieren, dass sie das Endergebnis beeinflussen:

   >>> 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 == 1.0
   False
   >>> sum([0.1] * 10) == 1.0
   True

Die  "math.fsum()"  geht noch einen Schritt weiter und verfolgt alle
"verlorenen Ziffern", während Werte zu einer laufenden Summe addiert
werden, sodass das Ergebnis nur einmal gerundet wird. Dies ist
langsamer als  "sum()" , liefert jedoch in seltenen Fällen, in denen
sich Eingaben großer Größenordnungen weitgehend gegenseitig aufheben
und eine Endsumme nahe Null ergibt, genauere Ergebnisse:

   >>> arr = [-0.10430216751806065, -266310978.67179024, 143401161448607.16,
   ...        -143401161400469.7, 266262841.31058735, -0.003244936839808227]
   >>> float(sum(map(Fraction, arr)))   # Exakte Summierung mit einmaliger Rundung
   8.042173697819788e-13
   >>> math.fsum(arr)                   # Einmalige Rundung
   8.042173697819788e-13
   >>> sum(arr)                         # Mehrfache Rundungen bei erweiterter Genauigkeit
   8.0421778034628478e-13
   >>> total = 0.0
   >>> for x in arr:
   ...     total += x                   # Mehrfache Rundungen in Standardgenauigkeit
   ...
   >>> total                            # Die einfache Addition liefert keine korrekten Ziffern!
   -0.0051575902860057365


15.1. Darstellungsfehler
========================

In diesem Abschnitt wird das Beispiel "0.1" ausführlich erläutert und
gezeigt, wie du Fälle wie diesen selbst genau analysieren kannst.
Grundkenntnisse über die binäre Gleitkommadarstellung werden
vorausgesetzt.

*Darstellungsfehler* bezieht sich auf die Tatsache, dass einige
(eigentlich die meisten) Dezimalbrüche nicht exakt als Binärbrüche
(Basis 2) dargestellt werden können. Dies ist der Hauptgrund dafür,
dass Python (oder Perl, C, C++, Java, Fortran und viele andere) oft
nicht genau die Dezimalzahl anzeigen, die man erwartet.

Warum ist das so?  1/10 lässt sich nicht exakt als binärer Bruch
darstellen.  Seit mindestens dem Jahr 2000 verwenden fast alle Rechner
die binäre IEEE-754-Gleitkommaarithmetik, und fast alle Plattformen
ordnen Python-Float-Werte den IEEE-754-binary64-Werten mit "doppelter
Genauigkeit" zu.  IEEE-754-binary64-Werte enthalten 53 Bits an
Genauigkeit, daher versucht der Computer bei der Eingabe, 0,1 in den
nächstgelegenen Bruch der Form *J*/2***N* umzuwandeln, wobei *J* eine
Ganzzahl mit genau 53 Bits ist. Umformulierung

   1 / 10 ~= J / (2**N)

als

   J ~= 2**N / 10

und unter Berücksichtigung, dass J genau 53 Bits hat (d. h. >= 2**52,
aber < 2**53), ist der beste Wert für N gleich 56:

   >>> 2**52 <=  2**56 // 10  < 2**53
   True

Das heißt, 56 ist der einzige Wert für *N*, bei dem *J* genau 53 Bits
hat. Der bestmögliche Wert für *J* ist dann dieser Quotient, gerundet:

   >>> q, r = divmod(2**56, 10)
   >>> r
   6

Da der Rest mehr als die Hälfte von 10 beträgt, erhält man die beste
Annäherung durch Aufrunden:

   >>> q+1
   7205759403792794

Daher lautet die bestmögliche Annäherung an 1/10 in IEEE
754-Doppelnauigkeit:

   7205759403792794 / 2 ** 56

Durch Division von Zähler und Nenner durch zwei lässt sich der Bruch
wie folgt vereinfachen:

   3602879701896397 / 2 ** 55

Beachte, dass dieser Wert, da wir aufgerundet haben, tatsächlich etwas
größer als 1/10 ist; hätten wir nicht aufgerundet, wäre der Quotient
etwas kleiner als 1/10 gewesen. Aber auf keinen Fall kann er *genau*
1/10 betragen!

Der Computer "sieht" also niemals 1/10: Was er sieht, ist genau der
oben angegebene Bruch, die bestmögliche IEEE-754-Doppelgenauigkeits-
Annäherung, die er erzielen kann:

   >>> 0.1 * 2 ** 55
   3602879701896397.0

Wenn wir diesen Bruch mit 10**55 multiplizieren, erhalten wir den Wert
mit 55 Dezimalstellen:

   >>> 3602879701896397 * 10 ** 55 // 2 ** 55
   1000000000000000055511151231257827021181583404541015625

Das bedeutet, dass die im Computer gespeicherte exakte Zahl dem
Dezimalwert 0,1000000000000000055511151231257827021181583404541015625
entspricht. Anstatt den vollständigen Dezimalwert anzuzeigen, runden
viele Sprachen (einschließlich älterer Python-Versionen) das Ergebnis
auf 17 signifikante Stellen:

   >>> format(0.1, '.17f')
   '0.10000000000000001'

Die Module  "fractions"  und  "decimal"  erleichtern diese
Berechnungen:

   >>> from decimal import Decimal
   >>> from fractions import Fraction

   >>> Fraction.from_float(0.1)
   Fraction(3602879701896397, 36028797018963968)

   >>> (0.1).as_integer_ratio()
   (3602879701896397, 36028797018963968)

   >>> Decimal.from_float(0.1)
   Decimal('0,1000000000000000055511151231257827021181583404541015625')

   >>> format(Decimal.from_float(0.1), '.17')
   '0.10000000000000001'
