Gibt es eine Möglichkeit, mit der i387 fsqrt-Anweisung eine korrekte Rundung zu erhalten?

Question

Gibt es eine Möglichkeit, mit der i387 fsqrt-Anweisung eine korrekte Rundung zu erhalten?...

...abgesehen von der Änderung des Präzisionsmodus im x87-Steuerwort - ich weiß, dass das möglich ist, aber es ist keine vernünftige Lösung, da es unangenehme Probleme mit dem Wiedereintrittstyp gibt, bei denen der Präzisionsmodus falsch ist, wenn der sqrt-Vorgang unterbrochen wird.

Das Problem, mit dem ich mich befasse, ist wie folgt:der x87 fsqrt opcode führt eine korrekt gerundete (gemäß IEEE 754) Quadratwurzeloperation in der Genauigkeit der FPU-Register durch, von der ich annehme, dass sie eine erweiterte Genauigkeit (80 Bit) aufweist.Ich möchte es jedoch verwenden, um effiziente Quadratwurzelfunktionen mit einfacher und doppelter Genauigkeit zu implementieren, wobei die Ergebnisse korrekt gerundet sind (gemäß dem aktuellen Rundungsmodus).Da das Ergebnis eine übermäßige Genauigkeit aufweist, rundet der zweite Schritt der Konvertierung des Ergebnisses in einfache oder doppelte Genauigkeit erneut ab, wobei möglicherweise ein nicht korrekt gerundetes Ergebnis verbleibt.

Bei einigen Operationen ist es möglich, dies mit Verzerrungen zu umgehen.Zum Beispiel kann ich übermäßige Genauigkeit bei den Additionsergebnissen vermeiden, indem ich eine Vorspannung in Form einer Zweierpotenz hinzufüge, die die 52 signifikanten Bits eines Werts mit doppelter Genauigkeit in die letzten 52 Bits der 63-Bit-Mantisse mit erweiterter Genauigkeit zwingt.Aber ich sehe keinen offensichtlichen Weg, einen solchen Trick mit der Quadratwurzel zu machen.

Irgendwelche cleveren Ideen?

(Auch mit C gekennzeichnet, da die beabsichtigte Anwendung die Implementierung von C ist.) sqrt und sqrtf Funktion.)

Answer 1

Lassen Sie uns zuerst das Offensichtliche aus dem Weg räumen:sie sollten SSE anstelle von x87 verwenden.Die SSW sqrtss und sqrtsd anweisungen machen genau das, was Sie wollen, werden auf allen modernen x86-Systemen unterstützt und sind auch deutlich schneller.

Wenn Sie nun darauf bestehen, x87 zu verwenden, beginne ich mit den guten Nachrichten:sie müssen nichts für Float tun.Du brauchst 2p + 2 bits, um eine korrekt gerundete Quadratwurzel in einem p-Bit-Gleitkommaformat zu berechnen.Da 80 > 2*24 + 2, wird die zusätzliche Rundung auf einfache Genauigkeit immer korrekt gerundet, und Sie haben eine korrekt gerundete Quadratwurzel.

Jetzt die schlechten Nachrichten: 80 < 2*53 + 2, also kein Glück für doppelte Genauigkeit.Ich kann verschiedene Problemumgehungen vorschlagen;hier ist eine schöne einfache von meinem Kopf.

1. lassen y = round_to_double(x87_square_root(x));
2. verwenden Sie ein Dekker-Produkt (Kopf-Schwanz) zur Berechnung a und b so dass y*y = a + b genau.
3. berechnen Sie das Residuum r = x - a - b.
4. if (r == 0) return y
5. if (r > 0), lassen y1 = y + 1 ulp, und berechnen a1, b1 s.t. y1*y1 = a1 + b1.Vergleichen r1 = x - a1 - b1 zu r, und entweder zurückgeben y oder y1, abhängig davon, welches das kleinere Residuum hat (oder das mit null niederwertigem Bit, wenn die Residuen gleich groß sind).
6. if (r < 0), machen Sie dasselbe für y1 = y - 1 ulp.

Diese Vorgehensweise behandelt nur den Standardrundungsmodus;in den gerichteten Rundungsmodi ist es jedoch richtig, einfach auf das Zielformat zu runden.

Answer 2

OK, ich denke, ich habe eine bessere Lösung:

1. Berechnen y=sqrt(x) in erweiterter Präzision (fsqrt).
2. Wenn die letzten 11 Bits nicht sind 0x400, einfach in doppelte Genauigkeit konvertieren und zurückgeben.
3. Hinzufügen 0x100-(fpu_status_word&0x200) auf das niedrige Wort der erweiterten Präzisionsdarstellung.
4. In doppelte Genauigkeit konvertieren und zurückgeben.

Schritt 3 basiert auf der Tatsache, dass das C1-Bit (0x200) des Statusworts genau dann 1 ist, wenn fsqrtdas Ergebnis wurde aufgerundet.Dies ist gültig, weil aufgrund des Tests in Schritt 2, x war kein perfektes Quadrat;wenn es ein perfektes Quadrat wäre, y hätte keine Bits mit doppelter Genauigkeit.

Es kann schneller sein, Schritt 3 mit einem bedingten Gleitkommabetrieb auszuführen, anstatt an der Bitdarstellung zu arbeiten und neu zu laden.

Hier ist der Code (scheint in allen Fällen zu funktionieren):

sqrt:
    fldl 4(%esp)
    fsqrt
    fstsw %ax
    sub $12,%esp
    fld %st(0)
    fstpt (%esp)
    mov (%esp),%ecx
    and $0x7ff,%ecx
    cmp $0x400,%ecx
    jnz 1f
    and $0x200,%eax
    sub $0x100,%eax
    sub %eax,(%esp)
    fstp %st(0)
    fldt (%esp)
1:  add $12,%esp
    fstpl 4(%esp)
    fldl 4(%esp)
    ret

Answer 3

Es ist möglicherweise nicht das, was Sie wollen, da es die Vorteile des 387 nicht nutzt fsqrt Anleitung, aber es gibt eine überraschend effiziente sqrtf(float) In glibc implementiert mit 32-Bit-Integer-Arithmetik.Es verarbeitet sogar NaNs, Infs und Subnormals korrekt – es könnte möglich sein, einige dieser Prüfungen mit echten x87-Anweisungen/FP-Steuerwort-Flags zu eliminieren.sehen: glibc-2.14/sysdeps/ieee754/flt-32/e_sqrtf.c

Der dbl-64/e_sqrt.c Code ist nicht so freundlich.Es ist schwer, auf den ersten Blick zu erkennen, welche Annahmen getroffen werden.Kurioserweise ist der i386 der Bibliothek sqrt[f|l] Implementierungen rufen Sie einfach an fsqrt, aber laden Sie den Wert anders. flds für SP, fldl für DP.