Existe-t-il un moyen d'obtenir un arrondi correct avec l'instruction i387 fsqrt ?

Question

Existe-t-il un moyen d'obtenir un arrondi correct avec l'instruction i387 fsqrt ?...

...à part changer le mode de précision dans le mot de contrôle x87 - je sais que c'est possible, mais ce n'est pas une solution raisonnable car il présente de vilains problèmes de type réentrance où le mode de précision sera erroné si l'opération sqrt est interrompue.

Le problème auquel je suis confronté est le suivant :le x87 fsqrt L'opcode effectue une opération de racine carrée correctement arrondie (conformément à la norme IEEE 754) dans la précision des registres fpu, ce qui, je suppose, est une précision étendue (80 bits).Cependant, je souhaite l'utiliser pour implémenter des fonctions de racine carrée efficaces en simple et double précision avec les résultats correctement arrondis (selon le mode d'arrondi actuel).Étant donné que le résultat présente une précision excessive, la deuxième étape de conversion du résultat en simple ou double précision est à nouveau arrondie, laissant éventuellement un résultat mal arrondi.

Avec certaines opérations, il est possible de contourner ce problème avec des biais.Par exemple, je peux éviter un excès de précision dans les résultats de l'addition en ajoutant un biais sous la forme d'une puissance de deux qui force les 52 bits significatifs d'une valeur double précision dans les 52 derniers bits de la mantisse à précision étendue de 63 bits. .Mais je ne vois aucun moyen évident de faire une telle astuce avec une racine carrée.

Des idées astucieuses ?

(Également étiqueté C car l'application prévue est l'implémentation du C sqrt et sqrtf les fonctions.)

Answer 1

Tout d’abord, éliminons l’évidence :vous devriez utiliser SSE au lieu de x87.L'ESS sqrtss et sqrtsd les instructions font exactement ce que vous voulez, sont prises en charge sur tous les systèmes x86 modernes et sont également beaucoup plus rapides.

Maintenant, si vous insistez pour utiliser x87, je vais commencer par la bonne nouvelle :vous n'avez rien à faire pour float.Vous avez besoin 2p + 2 bits pour calculer une racine carrée correctement arrondie dans un format à virgule flottante p-bits.Parce que 80 > 2*24 + 2, l'arrondi supplémentaire en simple précision arrondira toujours correctement et vous obtenez une racine carrée correctement arrondie.

Maintenant la mauvaise nouvelle : 80 < 2*53 + 2, donc pas de chance pour la double précision.Je peux suggérer plusieurs solutions de contournement ;en voici une belle et facile qui me vient à l'esprit.

1. laisser y = round_to_double(x87_square_root(x));
2. utiliser un produit Dekker (tête-queue) pour calculer a et b tel que y*y = a + b exactement.
3. calculer le résidu r = x - a - b.
4. if (r == 0) return y
5. if (r > 0), laisser y1 = y + 1 ulp, et calculer a1, b1 St. y1*y1 = a1 + b1.Comparer r1 = x - a1 - b1 à r, et revenez soit y ou y1, selon celui qui a le plus petit résidu (ou celui avec zéro bit de poids faible, si les résidus sont égaux en ampleur).
6. if (r < 0), fais la même chose pour y1 = y - 1 ulp.

Cette procédure ne gère que le mode d'arrondi par défaut ;cependant, dans les modes d'arrondi dirigé, le simple fait d'arrondir au format de destination fait la bonne chose.

Answer 2

OK, je pense que j'ai une meilleure solution :

1. Calculer y=sqrt(x) en précision étendue (fsqrt).
2. Si les 11 derniers bits ne sont pas 0x400, convertissez simplement en double précision et revenez.
3. Ajouter 0x100-(fpu_status_word&0x200) au mot bas de la représentation à précision étendue.
4. Convertissez en double précision et revenez.

L'étape 3 repose sur le fait que le bit C1 (0x200) du mot d'état est à 1 si et seulement si fsqrtLe résultat a été arrondi.Ceci est valable car, grâce au test de l'étape 2, x ce n'était pas un carré parfait ;si c'était un carré parfait, y n'aurait pas de bits au-delà de la double précision.

Il peut être plus rapide d'effectuer l'étape 3 avec une virgule flottante conditionnelle plutôt que de travailler sur la représentation des bits et le rechargement.

Voici le code (semble fonctionner dans tous les cas) :

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

Ce n'est peut-être pas ce que vous souhaitez, car il ne profite pas du 387 fsqrt instruction, mais il y a un étonnamment efficace sqrtf(float) dans glibc implémenté avec l'arithmétique entière de 32 bits.Il gère même correctement les NaN, les Infs et les anormalités - il pourrait être possible d'éliminer certaines de ces vérifications avec de véritables instructions x87 / indicateurs de mots de contrôle FP.voir: glibc-2.14/sysdeps/ieee754/flt-32/e_sqrtf.c

Le dbl-64/e_sqrt.c le code n'est pas si convivial.Il est difficile de déterminer d’un seul coup d’œil quelles sont les hypothèses qui sont formulées.Curieusement, le i386 de la bibliothèque sqrt[f|l] les implémentations appellent simplement fsqrt, mais chargez la valeur différemment. flds pour SP, fldl pour DP.