Rechenfehler mit Double-Gleitkommewerten

[MitjaStachowiak]

Hallo,
folgender Code

Code: Alles auswählen

var
  a, b : double;
  m : QWord;
  e : double; 

m := 400001;
e := 1e-9;
a := e * m;
m := 40000100000000003;
e := 1e-20;
b := e * m;

if (b < a) then
  WriteLn('Weird!');

Gibt "Weird" zurück. b müsste ja eigentlich etwas größer sein, als a. Durch begrenzte Nachkommastellen könnte ich noch verstehen, wenn es gleich groß ist. Aber wie ein b<a zustande kommt, erschließt sich mir nicht.

Wenn ich den Datentyp von e in extended ändere, dann ist b echt größer als a. Double hat im Ergebnis nicht zu wenig stellen... Wenn ich für e extendet verwende aber double(e)*m rechne, ist es wieder falsch. Also es muss irgendwie an der Multiplikation mit Double-Genauigkeit liegen.

Ist mir aufgefallen, weil meine selbst programmierte StrToFloat-Routine in einigen Fällen nicht das gleiche ergibt, wie .. naja die Konvertierung von im Code konstant vorgegebenen Werten oder dem StrToFloat von Lazarus. Mit Extended Precision für die Exponententabelle stimmen die Ergebnisse exakt überein.

[Jorg3000]

Hi!
Ohne dass ich auch nur den Hauch einer Ahnung habe, wie Fließkomma-Arithmetik in der CPU durchgeführt wird, trotzdem kurz mein Gedankengang:

Der Wert von m hat 17 Dezimalstellen vor dem Komma, und e:=1e-20 hat 20 Nachkommastellen. Würde man versuchen beide Werte temporär auf die gleiche Exponenten-Basis zu schieben, vielleicht weil das für eine Multiplikation intern erforderlich sein könnte (keine Ahnung), bräuchte man eine Genauigkeit von 37 Dezimalstellen. Dies ist jedoch weder bei Double (~16 Stellen) noch bei Extended (~19 Stellen) gegeben.
Vielleicht ist zum Multiplizieren aber gar keine gleiche Exponenten-Basis nötig.

Zweiter Gedanke, jetzt sehe ich es: Der Wert von m sprengt die Double-Genauigkeit, denn er hat 17 Stellen und Double ist nur auf 16 Stellen genau. Dadurch wird quasi aus der letzten Ziffer 3 eine Null (dezimal-bildlich gesprochen). Könnte dies schon für die Abweichung ausreichend sein?

Grüße, Jörg

[six1]

Das kommt daher, dass nicht jede beliebige Zahl des Dezimalsystems im Binärsystem darstellbar ist und du diesen Fehler mit sehr großen Zahlen auch noch multiplizierst.

[MitjaStachowiak]

Naja, die Gleitkommezahlen bestehen aus einer Basis, einem Exponenten als Integer mit Bias und einem Vorzeichen-Bit (deswegen ist +0 im bitweisen Vergleich nicht gleich -0).

Außer bei minimalem Expontenten berechnet sich der Wert als 1,mmmmm * 2^(e-bias) * [-1, wenn sign-bit gesetzt]. Die m sind die Bits der Mantisse.

Für eine Multiplikation werden die 1,mmm-Werte multipliziert und die Exponenten addiert. Ein Expontentenausgleich ist nur bei Addition erforderlich. Der Fehler kommt, wie six1 beschrieben hat, wohl daher, das 10^-20 nicht direkt in eine Basis 2^ex umrechenbar ist, sondern man unter Zuhilfenahme der Mantisse einen Wert nahe an 10^-20 verwenden muss, der hier wohl etwas darunter liegt und so der Fehler entsteht.

Ich bekomme die Werte von einer externen Anwendung nunmal als Strings geliefert und das sind sortierte Werte einer X-Achse. Durch den Konvertierungsfehler ändert sich aber plötzlich die Reihenfolge, sodass bestimmte Codes, die sich auf eine Sortierung verlassen, nicht mehr zuverlässig arbeiten.

Im konkreten Fall klappt es mit Extendet Precision - allerdings wirft das die Frage auf, ob Programme, die StrToFloat machen auf Plattformen, die kein Extended können, noch sicher funktionieren.? Oder ob bei anderen Zahlenkonstellationen vielleicht nicht auch mit Extended ein solches Sortierungsproblem auftreten kann? Und ob strToFloat(floatToStr(f)) = f garantiert immer true ergibt...

Jedenfalls speichere ich meine Ergebnisse wenn möglich lieber als binäre Buffer oder Base64 davon. Alleine weil dieses Umgerechne elend lahm ist.

[wp_xyz]

Also wenn man so knapp an der Genauigkeitgrenze der Gleitkommazahlen arbeitet, hat man sowieso ein Problem, das früher oder später an die Oberfläche kommt. Welche Werte stehen denn an der x-Achse?

[MitjaStachowiak]

Das sind einfach Zeitschritte einer SPICE-Simulation. Wenn irgendwo ein Transistor schaltet, dann macht die Simulation in der Nähe des Schaltvorgangs sehr viele, sehr kleine Rechenschritte...

[wp_xyz]

Und wie macht es dann SPICE? Die können doch unmöglich an der Genauigkeitgrenze der Gleitkommazahlen rechnen? Arbeiten die intern mit Nanosekunden oder Pikosekunden, so dass sie mit größeren Werten zu tun haben, und wandeln für die Ausgabe in Sekunden um (aber dann hätten sie letzendlich dasselbe Problem wie du)? Kannst du SPICE überreden, die Originalzahlen auszugeben?

[fliegermichl]

Gibt es da eigentlich eine Unit / Package, welches diese Beschränkungen wegen der binären Zahlendarstellung/Verwaltung nicht hat.

In PHP zum Beispiel kann man mit beliebig großen/kleinen Zahlen rechnen.

In meinem CAD Programm bin ich auch häufig schon über diese Probleme gestolpert und die Kundschaft sagt. "Mir doch egal wie dein blöder Computer rechnet. Wenn er falsche Ergebnisse liefert, ist er Mist"

und da muß ich denen leider Recht geben.

[PascalDragon]

Wenn ich den Datentyp von e in extended ändere, dann ist b echt größer als a. Double hat im Ergebnis nicht zu wenig stellen... Wenn ich für e extendet verwende aber double(e)*m rechne, ist es wieder falsch. Also es muss irgendwie an der Multiplikation mit Double-Genauigkeit liegen.

Du hast keinerlei Aussage darüber getroffen welche Plattform du verwendest. Auf x86_64-win64 und allen nicht-x86 Plattformen ist Extended nur ein Alias für Double.

Zum Beispiel unter x86_64-win64 bekomme ich die Ausgabe Weird unter x86_64-linux jedoch nicht.

Gibt es da eigentlich eine Unit / Package, welches diese Beschränkungen wegen der binären Zahlendarstellung/Verwaltung nicht hat.

Das Paket gmp zum Beispiel, welches Teil von FPC ist. Es benötigt jedoch eine externe Bibliothek.

[kupferstecher]

Ich bekomme die Werte von einer externen Anwendung nunmal als Strings geliefert und das sind sortierte Werte einer X-Achse. Durch den Konvertierungsfehler ändert sich aber plötzlich die Reihenfolge, sodass bestimmte Codes, die sich auf eine Sortierung verlassen, nicht mehr zuverlässig arbeiten.

Bei Taschenrechnern ist es wohl so, dass mit verdeckten Stellen gerechnet wird und die Anzeige dann daraus gerundet wird. Das hat dann den Effekt dass 1/3= 0,33333333 und anschließendes * 3 wieder 1 ergibt anstatt 0,99999999, was bei händisch eingegebenen 0,33333333 der Fall wäre. Vielleicht kannst du etwas ähnliches machen. Wenn man die letzten 2 Stellen wegrundet, sollte es zumindest nicht mehr zur Vertauschung der Reihenfolge kommen.

[MitjaStachowiak]

Kannst du SPICE überreden, die Originalzahlen auszugeben?

Naja, das ist PLECS, um genau zu sein, keine klassische Spice-Simulation, aber so ähnlich. Man kann wohl irgendwie die Ausgabe auch als Matlab-File bekommen, wo dann Binärdaten drin stehen. Aber da bin ich noch nicht durchgestiegen.

Ich arbeite auf x86_64 Linux. Werde also in Annahme, dass ohne Extended Precision die Hin- und Rückkonvertierung von Double-Werten fehlerhaft sein kann, das in der Doku erwähnen.

[Soner]

Teste dein Programm mit 32Bit-Compiler. Ich hatte letztes Jahr ähnliches Problem. Ich bin bei 32Bit-Compiler geblieben.

Siehst du dieses Beispiel von hier:

Code: Alles auswählen

var
  f: double;
  n: integer = 1758;
  m: integer = 0;
begin
  f := n * 1.2E6 + (2*m+1) * 50E3;

Der Wert von f ist bei 32Bit-FPC-Compiler 2109650000 (=richtig) und bei 64Bit-FPC-Compiler 2109650048 (=falsch).

[Niesi]

Interessant ...

Code: Alles auswählen

procedure tMainForm.Test;

var
  f: extended;
  i: integer = 1758;
  m: integer = 0;

begin                                   
  f := i * 1200000 + (2 * m + 1) * 50000;    { result ok }
  //f := i * 1.2e6 + (2 * m + 1) * 50e3;     { wrong result }
  //f := i * 1.2e6 + (2 * m + 1) * 50000;    { wrong result }
  //f := i * 12e5 + (2 * m + 1) * 50000;     { wrong result }
  //f := i * 1200000 + (2 * m + 1) * 50e3;   { wrong result }
  //f := i * 1200000 + (2 * m + 1) * 5e4;    { wrong result }
  MyLabel.Caption := Float2Str(f, 18, 0);
end;

So stimmt das Ergebnis, alles an Exponentialdarstellung macht es "kaputt" ...

[wp_xyz]

Lt dem zitierten wiki-Artikel interpretiert der Compiler Float-Konstanten ggfs als Single. Warum das denn? Bei Integern ist er übervorsichtig und rechnet im Bedarfsfall mit Int64, aber bei floats nimmt er es nicht so genau.

Jedenfalls kann man durch einen expliziten Type-Cast nach double das Problem auf 64-bit-Windows beheben:

Code: Alles auswählen

program Project1;
var
  f: double;
  n: integer = 1758;
  m: integer = 0;
begin
  f := n * 1.2E6 + (2*m+1) * 50E3;
  WriteLn(f:0:0);     // <----- falsches Ergebnis

  f := n * double(1.2E6) + (2*m+1) * double(50E3);
  WriteLn(f:0:0);     // <----- richtiges Ergebnis

  ReadLn;
end. 

Und der in dem Artikel erwähnte Trick mit {$excessprecision on} funktioniert auch nicht...

[Niesi]

Witzig ...

Code: Alles auswählen

procedure tMainForm.Test;

var
  f: extended;
  i: integer = 1758;
  m: integer = 0;
  g: extended;
  h: extended;

begin                                   
  //f := i * 1200000 + (2 * m + 1) * 50000;    { result ok }

  i := StrToInt(Edit_i.Text);
  m := StrToInt(Edit_m.Text);
  g := Str2Float(Edit_g.Text);
  h := Str2Float(Edit_h.Text);

  f := i * g + (2 * m + 1) * h;

  //f := i * 1.2e6 + (2 * m + 1) * 50e3;     { wrong result }
  //f := i * 1.2e6 + (2 * m + 1) * 50000;    { wrong result }
  //f := i * 12e5 + (2 * m + 1) * 50000;     { wrong result }
  //f := i * 1200000 + (2 * m + 1) * 50e3;   { wrong result }
  //f := i * 1200000 + (2 * m + 1) * 5e4;    { wrong result }
  MyLabel.Caption := Float2Str(f, 18, 0);
end;

So bekomme ich die Ergebnisse auch korrekt - was ist denn da los????