Ist meine fma () kaputt?

Es ist Cygwins Schuld. Oder genauer gesagt, die newlib C-Bibliothek, die sie verwendet. Es sagt ausdrücklich es versucht nicht einmal, fma() emulation richtig zu bekommen.

Die GNU C-Bibliothek hat seit 2015 eine korrekte Emulation für fast alle fma-Varianten. Details und die Patches, die für die Implementierung verwendet werden, finden Sie im Sourceware-Fehler 13304 .

Wenn die Effizienz kein Problem ist, würde ich einfach z. B.

Ich persönlich benutze Windows überhaupt nicht, aber wenn jemand das tut (Windows benutzt und die fma-Emulation benötigt), würde ich vorschlagen, dass sie versuchen, einen Patch upstream mit einem Link zu GNU C-Bibliothek Diskussion über korrekte fma-Emulation .

Was mich interessiert, ist, ob es möglich wäre, nur die niedrigen M Bits des Ergebnisses zu untersuchen (in der Rundung verworfen), um den korrekten Wert des ULP im Ergebnis zu bestimmen, und das Ergebnis unter Verwendung der einfachen a x b nextafter() ; anstatt Multiplecision-Arithmetik zu verwenden, um die gesamte Operation zu implementieren.

Edit: Nein, weil die Addition überlaufen kann und ein zusätzliches Bit als MSB des verworfenen Teils verloren geht. Aus diesem Grund müssen wir die gesamte Operation durchführen. Ein anderer Grund besteht darin, dass, wenn a × b und c verschiedene Vorzeichen haben, dann anstelle der Addition eine kleinere Größe von einer größeren Größe ( Ergebnis unter Verwendung des Größer-Vorzeichens), was dazu führen kann, dass mehrere hohe Bits in dem größeren gelöscht werden, und das betrifft, welche Bits des gesamten Ergebnisses in der Rundung fallengelassen werden.

Aber für die IEEE-754-Binär64% -Co_de% auf x86- und x86-64-Architekturen glaube ich, dass die fma-Emulation unter Verwendung von 64-Bit (Ganzzahl) -Registern und 128-Bit-Produkt immer noch ziemlich durchführbar ist. Ich werde mit einer Darstellung experimentieren, bei der niedrige 2 Bits in einem 64-Bit-Register für die Rundungsentscheidungsbits verwendet werden (LSB ist ein logisches OR aller abgelegten Bits), 53 Bits für die Mantisse und ein Übertragsbit, wobei 8 übrig bleibt unbenutzte und ignorierte hohe Bits. Die Rundung wird durchgeführt, wenn die vorzeichenlose ganzzahlige Mantisse in ein (64-Bit-) Double konvertiert wird. Wenn diese Experimente etwas Nützliches ergeben, werde ich sie hier beschreiben.

Erste Ergebnisse: double Emulation auf einem 32-Bit-System ist langsam. Das 80-Bit-Zeug auf der 387-FPU ist hier im Grunde nutzlos, und die Implementierung der 53 × 53-Bit-Multiplikation (und Bitverschiebung) auf einem 32-Bit-System ist nur ... nicht die Mühe wert. Der glibc fma() Emulationscode, der mit oben verlinkt ist, ist meiner Meinung nach gut genug.

Weitere Erkenntnisse: Die Behandlung nicht endlicher Werte ist böse . (Subnormale sind nur etwas störend, erfordern eine spezielle Behandlung (da das implizite MSB in der Mantisse dann Null ist).) Wenn eines der drei Argumente nicht endlich ist (unendlich oder irgendeine Form von NaN), dann wird fma() (nicht fusioniert) zurückgegeben ) ist die einzige vernünftige Option. Die Behandlung dieser Fälle erfordert zusätzliche Verzweigungen, die die Emulation verlangsamen.

Endgültige Entscheidung: Die Anzahl der Fälle, die in optimierter Weise behandelt werden sollen (anstatt den Multiplecision- "Gliedmaßen" -Ansatz zu verwenden, wie er in der Glibc-Emulation verwendet wird) ist groß genug, um diesen Ansatz den Aufwand nicht wert zu machen. Wenn jedes Glied 64-Bit ist, ist jedes von a , b und c über höchstens zwei Gliedmaßen verteilt, und a × b über drei Gliedmaßen. (Bei 32-Bit-Gliedmaßen sind dies nur 3 bzw. 5 Gliedmaßen.) Abhängig davon, ob a × b und c dasselbe haben oder verschiedene Zeichen, es gibt nur zwei grundsätzlich unterschiedliche Fälle zu behandeln - im Fall der verschiedenen Zeichen wird die Addition in Subtraktion (kleiner von größer, Ergebnis erhalten das gleiche Zeichen wie der größere Wert).

Kurz gesagt, der Multipräzisionsansatz ist besser. Die tatsächlich benötigte Genauigkeit ist sehr gut begrenzt und benötigt keine dynamische Zuordnung. Wenn das Produkt der Mantissen von a und b effizient berechnet werden kann, ist der Multipräzisionsteil darauf beschränkt, das Produkt zu halten und die Addition / Subtraktion zu handhaben. Die abschließende Rundung kann durchgeführt werden, indem das Ergebnis in eine 53-Bit-Mantisse, einen Exponenten und zwei extra niedrige Bits umgewandelt wird (wobei das höhere das signifikanteste Bit ist, das beim Runden verloren geht und das niedrigere ein OR des Rests der verlorenen Bits ist) die Rundung). Im Wesentlichen können die Schlüsseloperationen unter Verwendung von Ganzzahlen (oder SSE / AVX-Registern) durchgeführt werden, und die abschließende Umwandlung von einer 55-Bit-Mantisse in das Doppel behandelt die Rundung gemäß den aktuellen Regeln.