Calculer RD (sqrt (x)) avec une FPU en mode RU

Question

Les intervalles de bornes à virgule flottante peuvent être utilisés pour surestimer des ensembles de réels, à condition que la limite supérieure de tout intervalle de résultat soit calculée en arrondi vers le haut et la borne inférieure en arrondi vers le bas.

Une astuce recommandée consiste à calculer réellement la négation de la borne inférieure. Ceci permet de garder le FPU en arrondi vers le haut à tout moment (par exemple, "Handbook of Floating-Point Arithmetic", 2.9.2).

Cela fonctionne bien pour l'addition et pour la multiplication. D'autre part, l'opération de la racine carrée n'est pas symétrique dans les façons dont l'addition et la multiplication sont.

Il me semble que, pour calculer sqrt _RD, pour la limite inférieure, l'idiome suivant, en dépit de ses complications, pourrait être plus rapide sur une plate-forme ordinaire avec IEEE 754 double précision et FLT_EVAL_METHOD définie à 0 que changer le mode d'arrondi deux fois:

#include <fenv.h> 
#include <math.h> 
#pragma STDC FENV_ACCESS ON 
… 
/* assumes round-upwards */ 
double sqrt_rd(double l) { 
    feclearexcept(FE_INEXACT); 
    double candidate = sqrt(l); 
    if (fetestexcept(FE_INEXACT)) 
    return nextafter(candidate, 0); 
    return candidate; 
} 

Je me demande si c'est mieux, et si c'est le plus rapide. Comme une alternative possible, mais pas nécessairement la plus rapide, il me semble que FMA RU (candidat, candidat, -l) pourrait ne pas être toujours exacte (en raison de l'arrondi dirigé) mais pourrait être assez précise autour de 0 pour les suivantes au travail:

/* assumes round-upwards */ 
double sqrt_rd(double l) { 
    double candidate = sqrt(l); 
    if (fma(candidate, candidate, -l) != 0.0) 
    return nextafter(candidate, 0); 
    return candidate; 
} 

Quels autres moyens peu coûteux sont là pour détecter sqrt était inexacte? Quelle combinaison d'opérations en virgule flottante conduit au calcul le plus rapide de sqrt_rd sur un ensemble FPU moderne pour arrondir à la hausse?

Answer 1

Je pense que vous devriez être en mesure d'utiliser:

/* assumes round-upwards */ 
double sqrt_rd(double l) { 
    double u = sqrt(l); 
    double w = u*u; 
    if (w != l) 
    return nextafter(u, 0); 
    return u; 
} 

La justification ici étant que si u est inexacte, il sera strictement supérieur à √ l, ce qui implique que w> = u >l (puisque w est également calculé en mode RU). Et si u est exact, il en est de même pour w (puisque nous savons qu'il doit être représentable comme un double).

Answer 2

fma calcule le résultat avec une précision infinie, puis applique le mode d'arrondi.

Si votre candidat est trop grand, le résultat infiniment précis est supérieur à 0, et puisque vous êtes arrondi, il sera arrondi. Même si c'est seulement un tout petit peu plus grand que zéro. Pour vérifier cela, essayez d'abord l = 1 + 2eps, où (1 + eps) = sqrt (1 + 2eps + eps^2) est juste un tout petit peu trop grand; puis redescendre d'une puissance négative de 4 pour que l'eps^2 soit bien au-delà de la résolution des nombres dénormalisés, et vérifier cela aussi.