Fehler bei Fließkomma-Berechnungen auf x86-Architektur

Der Eintrag "Fließender Ärger" des Hostbloggers machte mich schon etwas stutzig. Zur Erinnerung, das war:

#include "stdio.h"

void main()
{
    float f = 35.91;
    printf("%.2f => %i / %i\n", f, (int)(f*100), (int)(f*1000));
    // erwartet: 35.91 => 3591 / 35910
    // tatsächlich: 35.91 => 3590 / 35909
};

Die in den Kommentaren geäußerte Vermutung, das in diesem speziellen Fall ein Float-Rundungfehler vorliegt, ist hier zwar richtig, es kann aber auch eine ganz andere Ursache haben.

Ich habe das Beispiel etwas abgeändert, um das zu verdeutlichen:

#include "stdio.h"

int main()
{
    float f = 35.91;
    float g = 35.91 * 100.0;

    printf("%.2f => %i / %i / %i \n", f, (int)(f*100.0), (int)(f*1000.0), (int)(g));
    // Ausgabe: 35.91 => 3590 / 35909 / 3591
};

Nanu, es macht offenbar einen Unterschied, ob das Berechnungsergebnis der Multiplikation einer Variablen zugewiesen wird oder ein Zwischenergebnis Verwendung findet! Wie kann das sein? Mal weiter abstrahieren und statt float double benutzen:

#include "stdio.h"

int main()
{
    double a=1.0f, b=3.0f, c;
    c = a / b;
    if( c == a / b )
    {
	printf("ok\n");
    }
    else
    {
	printf("nok\n");
    }
    // Ausgabe: nok
}

OK, es gibt Rundungsfehler bei Fließkommaberechnungen und man sollte zwei double-Zahlen nie mit "==" auf Gleichheit überprüfen. Aber hier wurde in C exakt die gleiche Berechnung durchgeführt, es sollte also auch in der Binärdarstellung der Zahlen exakt das Gleiche rauskommen.

Der Grund für das Verhalten ist in der x86/x87 – Hardware zu finden. Die FPU rechnet nämlich im Gegensatz zu den meisten andern Prozessoren nicht mit double, sondern mit extended precision, zumindest in der Standardeinstellung.

Wird was zurück in den Speicher geschrieben, dann wird das Ergebnis nach double konvertiert. Das Tückische dabei ist: Steht zufällig der andere Vergleichswert schon im CPU-Register in extended precision, dann geht der Vergleich schief. Da man in einer Hochsprache wie C nicht zuletzt durch Optimierungen im Compiler im Wesentlichen keinen direkten Einfluß auf die Datenhaltung (Speicher/CPU-Register)hat, führt das leicht zu nichtdeterministischem Verhalten.

Das Ganze scheint ein spezielles Problem der x86-Architektur zu sein. Man kann das Problem umgehen, indem man die FPU anweist, tatsächlich mit double zu rechnen.

Beim Testen mit den verschiedenen Compiler-Optimierungs-Optionen muß man übrigens auch aufpassen: Wenn man obiges Beispiel beim gcc mit "-O1" übersetzt, scheint auf einmal alles i.O. zu sein ("ok" wird ausgegeben). Schaut man sich aber den Assemblercode an (gcc -S datei.c), stellt man fest, daß der String "nok" nicht mehr vorkommt – der Compiler hat den else-Zweig einfach wegoptimiert.

Das folgende Beispiel bringt die Abhilfe:

#include "stdio.h"
void set_fpu (unsigned int mode)
{
  asm ("fldcw %0" : : "m" (*&mode));
}

int main()
{
    set_fpu(0x27F);

    double a=1.0f, b=3.0f, c;
    c = a / b;
    if( c == a / b )
    {
        printf("ok\n");
    }
    else
    {
        printf("nok\n");
    }
}

Das ist aber dann doch ganz schön erschreckend, wenn in einer Hochsprache ein Vergleich ein und desselben Datums mit sich selbst schiefgehen kann:

#include <stdlib.h>
#include <stdio.h>

#define axiom_order(a,b)  !(a < b && b < a) 
#define axiom_eq(a)       a == a 
#define third             ((double)atoi("1")/atoi("3"))

int main(void) {
        if (axiom_order(third, third))
                puts("ok");
        else
                puts("error");

        if (axiom_eq(third))
                puts("ok");
        else
                puts("error");

        return 0;
}

Beim gcc ist dieses Verhalten seit 2000 auch als Bug erfasst und mittlerweile (runterscrollen bis ins Jahr 2009) teilweise gelöst. Von dort stammt auch das letzte Beispiel.

Zum Weiterlesen und für weitere Beispiele siehe auch: