Alles über Gleitkommazahlen und Rundungsfehler

Einführung

Probieren Sie in GWBASIC.EXE folgende Berechnungen aus:

Ok
print 81!/10!;84!/10!
 8.100001  8.399999
Ok
print int(11! - (2.9 + 2.7 + 2.6) + .2)
 2 <= sollte doch 3 geben...
Ok
print int(2.9# + 2.8# + 3.1# + 0.2#)
 8 <= sollte doch 9 geben...
Ok
_

Benützer von QBASIC.EXE probieren als Rechnungen einmal

PRINT 82# / 10#
PRINT 83# / 10#

aus. In diesem Artikel erfahren Sie, wie überhaupt diese merkwürdigen Phänomene entstehen. Gleichzeitig zeige ich Ihnen anhand einer Rundungsroutine, welche ohne PRINT USING arbeitet, wie man solche Rundungsprobleme lösen kann.

Wie werden Gleitkommazahlen überhaupt gespeichert?

Microsoft BASIC speichert intern sämtliche Gleitkommazahlen (ganze Zahlen übrigens auch) immer als binäre Werte (Bitfolgen) ab, wie dies bei jeder Rechnerarchitektur und Programmiersprache üblich ist. Die Verarbeitung solcher Gleitkommawerte erfolgt über den mathematischen Co-Prozessor, welcher bei allen CPUs ab 80846DX und höher von Intel serienmässig integriert ist und früher noch als Option (der bekannte 8087/80287/80387-Chip) nachgerüstet werden musste. Die Interpreter-Versionen von Microsoft BASIC unterstützen allerdings diesen Hilfsprozessor noch nicht, so dass sie Gleitkommawerte über eine Software-Routine verarbeiten. Diese Tatsache können Sie übrigens beim Mandelbrotmengen-Berechungsprogramm sehr gut beobachten: Die BASIC-Version nimmt es wegen der fehlenden Co-Prozessorunterstützung auch auf einem modernen Pentium-System recht gemütlich, während die C++-Version wirklich die volle Leistung aus Ihrem Rechner herausholt.

Analyse im Detail

Mit Hilfe der Umwandlungsfunktionen MKS$() und MKD$() (C-Programmierer kennen dies als sog. Type Casting) können Sie mit einem sehr kleinen Programm Gleitkommazahlen analysieren. Im folgenden beziehe ich mich auf das in GWBASIC.EXE verwendete MBF-Format, doch mehr dazu noch später.

10 INPUT "Gleitkommazahl";F!
20 B$=""
30 FOR I%=1 TO 4
40 BY%=ASC(MID$(MKS$(F!),I%,1))
50 FOR J%=1 TO 8
60 B$=CHR$(48+BY% MOD 2)+B$
70 BY%=BY%\2
80 NEXT J%
90 NEXT I%
100 PRINT B$

Genereller Aufbau:

           EEEEEEEESMMMMMMMMMMMMMMMMMMMMMMM
           \__  __/|\_________  __________/
              \/   |          \/
           Exponent|       Mantisse
              Vorzeichen: 0=positiv, 1=negativ

Formel für einfache Genauigkeit:

Zahl! = (-1)^S * (1+MMMM...MMM/2^23) * 2^(EEEEEEEE-129)

wobei MMMM....MMM, S und EEEEEEEE den dezimalen Werten der binären Zahlen entsprechen, und zwar als ganze Zahlen interpretiert. Beispiel: 55:

55₁₀ = 110111₂ = +1.10111₂*10₂^101₂

Eine wichtige Anmerkung zur Mantisse: Die vorderste, signifikante und auch gleichzeitig die höchstwertige Stelle vor dem Komma beträgt immer 1, d.h. die Mantisse weist im Binärsystem immer die Form 1,xxx auf. Aus diesem Grund wird zur Vermeidung von Redundanz diese führende Eins weggelassen, so dass Sie anstelle von 11011100000000000000000 nur 10111000000000000000000 vorfinden.

Unendliche, periodische Binärbrüche

Und jetzt komme ich auf einen zentralen Punkt, welcher auch für Sie wichtig ist, sobald Sie beispielsweise in einer Finanzanwendung gebrochene Geldbeträge wie

10 GELDBETR1!=14.35
20 GELDBETR2!=171.29

verwenden. Dazu werfen wir einen Blick auf solche Dezimalbrüche:

Vom Dezimalsystem her sind Ihnen periodische Dezimalbrüche von Divisionen, welche nicht »aufgehen«, bestens bekannt:

1/3 = 0.333333333333... (Periode 3 wiederholt sich unendlich oft)

1/7 = 0.142857142857142857... (die Ziffernfolge 142857 wiederholt sich unendlich)

Falls früher Mathematik zu Ihren Lieblingsfächern an der Schule gehörte, mögen Sie sich vielleicht noch daran erinnern, dass eine Division nur dann endlich ausgeführt werden kann, wenn sich der Divisor aus den beiden Primfaktoren 2 und 5 zusammensetzt, was durch die Tatsache 2×5=10 beim Dezimalsystem (Primfaktoren-Zusammensetzung der Basis) begründet liegt.

Beim Binärsystem mit der Basis 2 steht dementsprechend nur der Primfaktor 2 zur Verfügung, so dass auch ein Faktor 5 genauso eine unendliche Periode bewirkt. Aus diesem Grund sind Ihr PC und natürlich auch sämtliche übrigen Rechnerarchitekturen, Betriebssysteme und Programmiersprachen nicht imstande, Dezimalbrüche exakt zu speichern. Sie arbeiten also immer mit Näherungswerten, was beispielsweise bei

PRINT 2.9 + 2.7 + 2.6

als Ergebnis 8.200001 anstelle von 8.2 zur Folge hat. Ziemlich fatale Folgen bewirkt dieses Beispiel:

PRINT INT(11 - (2.9 + 2.7 + 2.6) + .2)

Anstelle von 3 liefert Ihnen der PC nur 2, weil sich beim Addieren solche Rundungsfehler zusammengetragen haben.

Seien Sie daher in solchen Fällen immer vorsichtig mit den Ergebnissen speziell in Finanzanwendungen!

Wie wird eigentlich die Null dargestellt?

Bei der Analyse der Beispielwerte oben stellten wir fest, dass die redundante führende 1 vor dem Komma bei der Mantisse weggelassen wird. Da aber S und EEEEEEEE im Exponent der oben dargestellten Formel erscheinen, müsste es doch rein mathematisch gesehen unmöglich sein, den Wert Null darzustellen. Aber wie kann die Null dennoch dargestellt werden?

Ok
run
Gleitkommazahl? 0
00000000000000000000000000000000
Ok
_

Mit der oberen Formel erhalten wir strenggenommen

(-1)^0*(1+0/2^23)*2^(0-129) = 2^(-129)

Die Null ist jedoch als Spezialwert (engl. magic number) für den mathematischen Coprozessor bzw. Software-Routine definiert.

Doppeltgenaue Zahlen

Ein »schlauer« Kopf unter Ihnen würde als »Abhilfe« einfach hingehen und sämtliche Variablen und Konstanten durch doppeltgenaue Werte ersetzen. Jedoch löst dieses Vorgehen keinesfalls das Grundproblem, wie Sie sogleich sehen werden.

Analyse und Beispiele doppeltgenauer Werte

Dazu können Sie das obige Programm leicht modifizieren:

10 INPUT "Gleitkommazahl";F#
20 B$=""
30 FOR I%=1 TO 8
40 BY%=ASC(MID$(MKD$(F#),I%,1))
50 FOR J%=1 TO 8
60 B$=CHR$(48+BY% MOD 2)+B$
70 BY%=BY%\2
80 NEXT J%
90 NEXT I%
100 PRINT B$

Allgemeiner Aufbau: EEEEEEEESMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM

Der Aufbau ist grundsätzlich ganz ähnlich, ausserdem können Sie auch einen gleich langen Exponenten feststellen, was auch erklärt, warum 1.7D+38 die grösstmögliche Zahl auch bei den doppeltgenauen Werten darstellt. Bei fast allen übrigen Programmiersprachen und Mathematik-Coprozessoren werden jedoch 11 statt 8 Bits für den Exponenten verwendet, was dann auch die bekannte D+308/D-308-Limite begründet. Dies geht natürlich auf Kosten der Genauigkeit, weil in der Mantisse genau diese drei Bits fehlen. Doch wie sieht die Sache bei doppeltgenauen Dezimalbrüchen aus?

Analog den einfachgenauen Zahlen entstehen die Rundungsfehler genauso.

Umwandlung in die höhere Genauigkeit

Mit allem bisher Gesagten können Sie nun den folgenden Effekt analysieren und begründen:

Ok
print cdbl(.1)
 .1000000014901161
Ok
_

Wie kommt dieser Genauigkeitsfehler zustande? Dazu werfen wir zuerst einen Blick auf die Darstellung von 0.1 als einfachgenaue Zahl:

Die CDBL()-Funktion bekommt aber nur diesen binären Wert und kann die Genauigkeit nur durch Auffüllen von Nullen innerhalb der Mantisse erweitern:

CDBL(01111101010011001100110011001101) =
     0111110101001100110011001100110100000000000000000000000000000000
                                    ^\______________  ______________/
                                    |               \/
Aufrundung => daher 1.000000xxxx ___/    Auffüllung mit Nullen

Der genauere Wert für 0.1# wäre ansonsten

     0111110101001100110011001100110011001100110011001100110011001101

Nun sollte Ihnen wiederum klar sein, wie folgender Effekt zustande kommt:

Ok
print 2.9# + 2.8# + 3.1# + 0.2#
 9
Ok
print int(2.9# + 2.8# + 3.1# + 0.2#)
 8        \___________  __________/
Ok                    \/
_          derselbe Rechenausdruck

Auch hier sind wiederum verursachen nichtabbrechende Binärbrüche solche Rundungsfehler.

Nur GW-BASIC: Mathematische Funktionen mit doppelter Genauigkeit

Eine minimale Portion an Vorsicht ist bei der Verwendung mathematischer Funktionen angebracht:

Ok
print sqr(2#); sin(1#); log(3#)
 1.414214  .841471  1.098612
Ok
_

Hoppla! Nur einfache Genauigkeit? :-(

Verwenden Sie auf gar keinen Fall CDBL(SQR(2#)) wegen der vorhin beschriebenen Auffüllung mit Nullen!

Der Grund für dieses Verhalten liegt historisch begründet: Als GWBASIC.EXE im Jahre 1981 herauskam, war Arbeitsspeicher noch sündhaft teuer, so dass mit jedem einzelnen Byte sehr geizig umgegangen werden musste. Daher wurde GWBASIC.EXE der Schalter /D spendiert, mit welchem man die zugehörige Software-Routine optional laden kann.

Das dieses Relikt bis zur letzten Version geblieben ist, müssen wir GWBASIC.EXE mit angehängtem /D neustarten:

Ok
system

C:\BASICPRG>gwbasic /d

GW-BASIC 3.23
(C) Copyright Microsoft 1983,1984,1985,1986,1987,1988
60300 Byte frei
Ok
? sqr(2#);sin(1#);log(3#)
 1.414213562373095  .8414709848078965  1.09861228866811
Ok
_

Jetzt erst arbeiten alle mathematischen Funktionen auf Wunsch doppeltgenau.

Prüfabfrage für eigene Programme

Damit die Endbenutzer im Falle einer hochpräzisen CAD-Anwendung wie dem CAD-Zahnradgenerator nicht aus Versehen Resultate in der einfachen Genauigkeit erhalten, sollten Sie in solchen Fällen, eine Prüfabfrage zu Beginn einbauen:

10 ' Prüfabfrage
20 IF LEN(STR$(SIN(1#))) < 14 THEN PRINT "Bitte GW-BASIC mit angehängtem /D neustarten!": SYSTEM

Entwicklung einer Rundungsroutine

GWBASIC.EXE und auch QuickBASIC bieten normalerweise den PRINT USING-Befehl an, um nebst einer Rundung auch eine Formatierung zu bewirken:

PRINT USING "CHF ####.##"; Geldbetr!

Was Microsoft allerdings vergessen hat, ist eine FORMAT$-Funktion, wie sie nur MaxonBASIC auf dem Amiga kennt, welche funktionsmässig sprintf() aus ANSI C/C++ entspricht, also die formatierte Ausgabe als String zurückgibt.

Variante 1: Temporäre Datei

Eine mögliche Tricklösung wäre das Schreiben in eine temporäre Datei:

10 RANDOMIZE TIMER
20 INPUT "Wert"; W!
30 TP$ = ENVIRON$("TEMP")
40 IF TP$=;"" THEN PRINT "TEMP nicht gesetzt! Ergänzen Sie die AUTOEXEC.BAT!": SYSTEM
50 IF RIGHT$(TP$, 1) <> "\" THEN TP$ = TP$ + "\"
60 TP$ = TP$ + "~"
70 FOR I%=1 TO 3: TP$ = TP$ + CHR$(65 + CINT(INT(26! * RND))): NEXT I%
80 TP$ = TP$ + MID$(TIME$, 4, 2) + RIGHT$(TIME$, 2) + ".TMP"
90 OPEN TP$ FOR OUTPUT AS 1
100 ' Hier erfolgt die Formatierung
110 PRINT#1, USING "##########.##"; W!
120 CLOSE 1
130 OPEN TP$ FOR INPUT AS 1
140 LINE INPUT#1, F$
150 CLOSE 1
160 KILL TP$
170 ' weitere Verarbeitung von F$
180 PRINT "Ausgabe: >"; F$; "<"

Betreffend dem in Zeilen 30 bis 80 erzeugten temporären Dateinamen verweise ich Sie auf den Artikel über Umgebungsvariablen.

Als grossen Nachteil fällt Ihnen sicherlich auf, dass diese Methode nur sehr langsam arbeitet, da viele Festplattenzugriffe entstehen. Mit Hilfe eines virtuellen Laufwerks mit RAMDRIVE.SYS könnten Sie zwar Ihre Festplatte vor einer intensiveren Abnutzung verschonen, aber langsam bleibt diese Variante aufgrund der vielen Systemaufrufe weiterhin.

Variante 2: Mit Rundungsfunktion

Als ich noch auf dem Commodore 64 programmierte, war folgende Variante in vielen BASIC-Bücher sehr verbreitet und auch gleichzeitig die einzige Möglichkeit:

10 INPUT "Wert"; W!
20 F$ = STR$(INT(100! * W! + .5) / 100!)
30 IF MID$(F$, 2, 1) = "." THEN F$ = LEFT$(F$, 1) + "0" + MID$(F$, 2)
40 IF LEN(F$) > 13 THEN F$ = RIGHT$(F$, 13)
50 F$ = SPACE$(13 - LEN(F$)) + F$
60 ' weitere Verarbeitung von F$
70 PRINT "Ausgabe: >%quot;; F$; "<"

Der Algorithmus ist recht einfach: Die INT()-Funktion rundet immer auf die nächstkleinere ganze Zahl. Mit der Addition von 0.5 wird daher eine korrekte kaufmännische Rundung bewirkt. Und mit dem Faktor 100 mit anschliessender Rückdivision erfolgt die Rundung auf eine beliebig wählbare Granularität ähnlich der VRUNDEN()-Funktion aus Microsoft Excel. Doch genau bei dieser Rückdivision kann etwas Fatales geschehen, wenn Sie dieses Programm einmal mit 9.1 starten:

Ok
run
Wert? 9.1
Ausgabe: >     9.100001<
Ok
_

Aufgrund der bisherigen Ausführungen sollte es Ihnen nicht mehr schwerfallen, diesen Effekt begründen zu können! Zugleich versagt diese Lösung auch auf .x0 oder gar .00 endenden Werten, wobei dieses Problem mit geeigneten Kommasuch-Zeichenkettenoperationen beseitigt werden könnte.

Variante 3: Rückdivisionsfreie Rundung

Das Problem mit den nichtabbrechenden Binärbrüchen kann durch Vermeidung der Rückdivision beseitigt werden, in dem Sie das Komma nur rein optisch durch Zeichenkettenoperationen einfügen:

10 INPUT "Wert"; W!
20 F$ = STR$(INT(100! * W! + .5))
30 IF LEN(F$) < 4 THEN F$ = LEFT$(F$, 1) + STRING$(4 - LEN(F$), "0") + MID$(F$, 2)
40 IF LEN(F$) > 12 THEN F$ = RIGHT$(F$, 12)
50 F$ = SPACE$(12 - LEN(F$)) + F$
60 F$ = LEFT$(F$, 10) + "." + RIGHT$(F$, 2)
70 ' weitere Verarbeitung von F$
80 PRINT "Ausgabe: >%quot;; F$; "<"

Kurz zum Algorithmus: In Zeile 30 erfolgt die eigentliche Rundung, und zwar liegt das Ergebnis zunächst ganzzahlig vor. Bevor wir in Zeile 50 mit den nötigen führenden Leerzeichen auffüllen bzw. in Zeile 40 abschneiden sowie in Zeile 60 den Dezimalpunkt optisch einfügen können, müssen wir mit Zeile 30 dafür sorgen, dass bei Bedarf genügend führende Nullen bei Werten wie 0.01 vorhanden sind.

Die Darstellungsnorm IEEE 754

Für die interne Darstellung von Gleitkommawerten existiert seit 1985 auch der vom Institute of Electrical and Electronics Engineers (IEEE) international festgelegter Standard IEEE 754, an dem sich sowohl QuickBASIC (auch QBASIC.EXE) als auch AmigaBASIC halten. Auch Intel hat sich bei der Entwicklung des mathematischen Coprozessors an diese Norm gehalten. Vorhin verwendete praktisch jeder Hersteller sein eigenes Format. So hatte auch Microsoft sein proprietäres Gleitkommazahlenformat Microsoft Binary Format (MBF) entwickelt und 1980 in GWBASIC.EXE eingesetzt, wie ich es immer in den vorherigen Beispielen gezeigt habe. Als Überbleibsel aus dieser Zeit befinden sich in QuickBASIC die bekannten 4 Funktionen MKSMBF$(), MKDMBF$(), CVSMBF() und CVDMBF(), mit welcher Sie im Falle der Migration eines alten GW-BASIC-Programms nach QuickBASIC die Dateikompatibilität sicherstellen können.

Aufbau von IEEE-Zahlen

Für das MBF-Format zeigte ich Ihnen den Aufbau EEEEEEEESMMMMMMMMMMMMMMMMMMMMMMM sowie die Formel

Zahl! = (-1)^S * (1+MMMM...MMM/2^23) * 2^(EEEEEEEE-129)

Die Zahl 129 wird in den einschlägigen Spezifikationen als sog. Bias bezeichnet. IEEE-Zahlen sind dabei etwas anders aufgebaut: Das Vorzeichenbit steht bereits am Anfang, danach folgt der Exponent mit etwas anderem Bias-Wert und zum Schluss folgt noch die Mantisse.

Ansonsten ist bei IEEE 754 die Bedeutung des Vorzeichen-Bits mit MBF identisch (0=+, 1=-), auch wird die immer vorhandene Eins vor dem Komma weggelassen, und ebenso wird auch der Wert Null als 00000000000000000000000000000000 bzw. 0000000000000000000000000000000000000000000000000000000000000000 aus lauter Nullen dargestellt. IEEE 754 spezifiziert darüber hinaus noch einige weitere solche Spezialwerte wie unendlich, ungültig usw. Einfachgenaue Zahlen belegen also jeweils 4 Bytes Speicher, während doppeltgenaue jeweils 8 Bytes belegen, was Ihnen sicherlich im Zusammenhang mit DIM und Feldvariablen bekannt sein dürfte.

Die beiden Analysierprogramme von vorher können Sie auch in QuickBASIC zur Analyse verwenden. Bei AmigaBASIC ist noch zusätzlich zu beachten, dass aufgrund der umgekehrten Byte-Reihenfolge der Motolora-Prozessorfamilie auch bei Gleitkommawerten die Bytes umgekehrt im Speicher liegen, also das höchstwertige Byte zuerst, ansonsten hat sich auch Commodore an den IEEE-Standard gehalten.

Weiterführende Informationen

Hier möchte ich Sie auf den Gleitkommazahlen-Artikel Q42980 von Microsoft selber verweisen.

Wieder zurück zur Übersicht