Pourquoi la fonction hypot() est-elle si lente ?

Question

J'ai fait quelques tests avec C++ hypot() y Java Math.hypot . Ils semblent tous deux nettement plus lents que sqrt(a*a + b*b) . Est-ce en raison d'une meilleure précision ? Quelle méthode pour calculer une hypoténuse hypot fonction utilise ? Étonnamment, je n'ai pas trouvé d'indication de mauvaises performances dans la documentation.

Answer 1

Ce n'est pas une simple fonction sqrt. Vous devriez consulter ce lien pour l'implémentation de l'algorithme : http://www.koders.com/c/fid7D3C8841ADC384A5F8DE0D081C88331E3909BF3A.aspx

Il a une boucle while pour vérifier la convergence

/* Slower but safer algorithm due to Moler and Morrison.  Never
         produces any intermediate result greater than roughly the
         larger of X and Y.  Should converge to machine-precision
         accuracy in 3 iterations.  */

      double r = ratio*ratio, t, s, p = abig, q = asmall;

      do {
        t = 4. + r;
        if (t == 4.)
          break;
        s = r / t;
        p += 2. * s * p;
        q *= s;
        r = (q / p) * (q / p);
      } while (1);

EDIT (Mise à jour de J.M) :

Ici est l'article original de Moler-Morrison, et aquí est une belle suite due à Dubrulle.

Answer 2

Voici une implémentation plus rapide, dont les résultats sont également plus proches de java.lang.Math.hypot : (NB : pour l'implémentation de Delorie, il faut ajouter la gestion des entrées NaN et +-Infini)

private static final double TWO_POW_450 = Double.longBitsToDouble(0x5C10000000000000L);
private static final double TWO_POW_N450 = Double.longBitsToDouble(0x23D0000000000000L);
private static final double TWO_POW_750 = Double.longBitsToDouble(0x6ED0000000000000L);
private static final double TWO_POW_N750 = Double.longBitsToDouble(0x1110000000000000L);
public static double hypot(double x, double y) {
    x = Math.abs(x);
    y = Math.abs(y);
    if (y < x) {
        double a = x;
        x = y;
        y = a;
    } else if (!(y >= x)) { // Testing if we have some NaN.
        if ((x == Double.POSITIVE_INFINITY) || (y == Double.POSITIVE_INFINITY)) {
            return Double.POSITIVE_INFINITY;
        } else {
            return Double.NaN;
        }
    }
    if (y-x == y) { // x too small to substract from y
        return y;
    } else {
        double factor;
        if (x > TWO_POW_450) { // 2^450 < x < y
            x *= TWO_POW_N750;
            y *= TWO_POW_N750;
            factor = TWO_POW_750;
        } else if (y < TWO_POW_N450) { // x < y < 2^-450
            x *= TWO_POW_750;
            y *= TWO_POW_750;
            factor = TWO_POW_N750;
        } else {
            factor = 1.0;
        }
        return factor * Math.sqrt(x*x+y*y);
    }
}

Answer 3

Je trouvais que Math.hypot() était d'une lenteur abyssale. J'ai découvert que je pouvais coder une version java rapide utilisant le même algorithme et produisant des résultats identiques. Cela peut être utilisé pour le remplacer.

/**
 * <b>hypot</b>
 * @param x
 * @param y
 * @return sqrt(x*x +y*y) without intermediate overflow or underflow. 
 * @Note {@link Math#hypot} is unnecessarily slow.  This returns the identical result to 
 * Math.hypot with reasonable run times (~40 nsec vs. 800 nsec). 
 * <p>The logic for computing z is copied from "Freely Distributable Math Library" 
 * fdlibm's e_hypot.c. This minimizes rounding error to provide 1 ulb accuracy.
 */
public static double hypot(double x, double y) {

    if (Double.isInfinite(x) || Double.isInfinite(y)) return Double.POSITIVE_INFINITY;
    if (Double.isNaN(x) || Double.isNaN(y)) return Double.NaN;

    x = Math.abs(x);
    y = Math.abs(y);

    if (x < y) {
        double d = x;
        x = y;
        y = d;
    }

    int xi = Math.getExponent(x);
    int yi = Math.getExponent(y);

    if (xi > yi + 27) return x;

    int bias = 0;
    if (xi > 510 || xi < -511) {
        bias = xi;
        x = Math.scalb(x, -bias);
        y = Math.scalb(y, -bias);           
    }

    // translated from "Freely Distributable Math Library" e_hypot.c to minimize rounding errors
    double z = 0; 
    if (x > 2*y) { 
        double x1 = Double.longBitsToDouble(Double.doubleToLongBits(x) & 0xffffffff00000000L);
        double x2 = x - x1;
        z = Math.sqrt(x1*x1 + (y*y + x2*(x+x1)));
    } else {
        double t = 2 * x;
        double t1 = Double.longBitsToDouble(Double.doubleToLongBits(t) & 0xffffffff00000000L);
        double t2 = t - t1;
        double y1 = Double.longBitsToDouble(Double.doubleToLongBits(y) & 0xffffffff00000000L);
        double y2 = y - y1;
        double x_y = x - y;
        z = Math.sqrt(t1*y1 + (x_y*x_y + (t1*y2 + t2*y))); // Note: 2*x*y <= x*x + y*y
    }

    if (bias == 0) {
        return z; 
    } else {
        return Math.scalb(z, bias);
    }
}

Answer 4

hypot() entraîne une surcharge pour éviter les débordements et les sous-exécutions dans les calculs intermédiaires, par rapport à l'implémentation naïve. sqrt(a*a+b*b) . Cela implique des opérations de mise à l'échelle. Dans les anciennes implémentations, la mise à l'échelle peut utiliser la division, qui est elle-même une opération plutôt lente. Même lorsque la mise à l'échelle est accomplie par la manipulation directe de l'exposant, elle peut être plutôt lente sur les anciennes architectures de processeur, car le transfert de données entre les UAL entières et les unités à virgule flottante était plutôt lent (impliquant par exemple un aller-retour dans la mémoire). Les branchements basés sur les données sont également communs à diverses approches de mise à l'échelle ; ils sont difficiles à prévoir.

Un objectif fréquent des concepteurs de bibliothèques mathématiques est d'obtenir un arrondi fidèle pour des fonctions mathématiques simples telles que hypot() c'est-à-dire une erreur maximale inférieure à 1 %. ulp . Cette amélioration de la précision par rapport à la solution naïve signifie que les calculs intermédiaires doivent être effectués avec une précision supérieure à la précision native. Une méthode classique consiste à diviser les opérandes en parties "haute" et "basse" et à simuler la précision étendue. Cela ajoute une surcharge pour le fractionnement et augmente le nombre d'opérations en virgule flottante. Enfin, la spécification ISO C99 pour hypot() (adoptée plus tard dans la norme C++) prescrit la gestion des NaN et des infinis qui ne découlent pas naturellement d'un calcul simple.

Un exemple représentatif des anciennes meilleures pratiques est __ieee754_hypot() en FDLIBM . Bien qu'il revendique une erreur maximale de 1 ulp, un test rapide contre une bibliothèque de précision arbitraire montre qu'il n'atteint pas réellement cet objectif, puisque des erreurs allant jusqu'à 1,15 ulp peuvent être observées.

Depuis que cette question a été posée, deux avancées dans les architectures de processeurs ont permis de rendre plus efficace hypot() mises en œuvre. L'une d'entre elles est l'introduction de la fonction de multiplication-addition fusionnée ( FMA ), qui utilise en interne le produit complet pour l'addition et n'applique qu'un seul arrondi à la fin. Cela évite souvent d'avoir à simuler une précision légèrement supérieure dans les calculs intermédiaires et, en combinant deux opérations en une seule instruction, améliore également les performances. L'autre avancée architecturale consiste à permettre aux opérations en virgule flottante et en nombre entier d'opérer sur le même ensemble de registres, ce qui rend la manipulation des bits des opérandes en virgule flottante peu coûteuse.

En conséquence, dans les bibliothèques mathématiques modernes à haute performance hypot() a généralement un débit deux fois moindre que celui de sqrt() comme on peut le voir dans le données de performance Intel publie pour MKL, par exemple.

Pour ses propres expériences, il est préférable de commencer par une mise en œuvre en simple précision de la fonction hypot() car cela permet de tester une plus grande partie de la combinaison totale de toutes les valeurs d'arguments possibles. Cela permet également d'explorer les compromis entre la vitesse et la précision lors de l'utilisation d'approximations de racine carrée réciproque de faible précision qui sont fournies par de nombreuses architectures de processeurs modernes. Un exemple de code pour une telle exploration est présenté ci-dessous.

Notez que les exigences relatives à la gestion des NaN imposées par la norme ISO C99 (et par extension, C++) sur fmin() fmax() ne correspondent pas toujours à la fonctionnalité des instructions machine minimum/maximum à virgule flottante, ce qui rend ces fonctions plus lentes qu'elles ne le pourraient. Puisque le code ci-dessous traite les NaNs séparément, il devrait être possible d'utiliser directement ces instructions machine. Le code devrait être compilé avec une optimisation complète mais aussi une stricte conformité IEEE-754.

Sur la base de tests effectués avec 50 milliards de paires d'arguments aléatoires, l'erreur maximale observée avec les variantes de code ci-dessous est la suivante : ( USE_BEEBE == 1 ) : 1,51 ulp ; ( USE_BEEBE == 0 && USE_SQRT == 1 ) : 1,02 ulp ; ( USE_BEEBE == 0 && USE_SQRT == 0 && USE_2ND_ITERATION == 0 ) : 2,77 ulp ; ( USE_BEEBE == 0 && USE_SQRT == 0 && USE_2ND_ITERATION == 1 ) : 1,02 ulp.

#include <stdint.h>
#include <string.h>
#include <float.h>
#include <math.h>
#include "immintrin.h"

uint32_t __float_as_uint32 (float a);
float __uint32_as_float (uint32_t a);
float sse_rsqrtf (float a);

float my_hypotf (float a, float b)
{
    float fa, fb, mn, mx, res, r, t;

    fa = fabsf (a);
    fb = fabsf (b);
    mx = fmaxf (fa, fb);
    mn = fminf (fa, fb);

#if USE_BEEBE
    /* Nelson H. F. Beebe, "The Mathematical Function Computation Handbook."
       Springer 2017
    */
    float s, c;

    r = mn / mx;
    t = fmaf (r, r, 1.0f);
    s = sqrtf (t);
    c = fmaf (-s, s, t) / (s + s);
    res = fmaf (mx, s, mx * c);
    if (mx == 0) res = mx; // fixup hypot(0,0)

#else // USE_BEEBE

    float scale_in, scale_out, s, v, w;
    uint32_t expo;

    /* compute scale factors */
    expo = __float_as_uint32 (mx) & 0xfc000000;
    scale_in = __uint32_as_float (0x7e000000 - expo);
    scale_out = __uint32_as_float (expo + 0x01000000);

    /* scale mx towards unity */
    mn = mn * scale_in;
    mx = mx * scale_in;

    /* mx in [2**-23, 2**6) */
    s = fmaf (mx, mx, mn * mn); // 0.75 ulp
#if USE_SQRT
    w = sqrtf (s);
#else // USE_SQRT
    r = sse_rsqrtf (s);
    /* use A. Schoenhage's coupled iteration for the square root */
    v = r * 0.5f;
    w = r * s;
    w = fmaf (fmaf (w, -w, s), v, w);
#if USE_2ND_ITERATION
    v = fmaf (fmaf (r, -w, 1), v, v);
    w = fmaf (fmaf (w, -w, s), v, w);
#endif // USE_2ND_ITERATION
    if (mx == 0) w = mx; // fixup hypot(0,0)
#endif // USE_SQRT

    /* reverse previous scaling */
    res = w * scale_out;
#endif // USE_BEEBE

    /* check for special cases: infinity and NaN */
    t = a + b;
    if (!(fabsf (t) <= INFINITY)) res = t; // isnan(t)
    if (mx == INFINITY) res = mx; // isinf(mx)
    return res;
}

uint32_t __float_as_uint32 (float a)
{
    uint32_t r;
    memcpy (&r, &a, sizeof r);
    return r;
}

float __uint32_as_float (uint32_t a)
{
    float r;
    memcpy (&r, &a, sizeof r);
    return r;
}

float sse_rsqrtf (float a)
{
    __m128 b, t;
    float res;
    int old_mxcsr;
    old_mxcsr = _mm_getcsr();
    _MM_SET_FLUSH_ZERO_MODE (_MM_FLUSH_ZERO_OFF); 
    _MM_SET_ROUNDING_MODE (_MM_ROUND_NEAREST);
    b = _mm_set_ss (a);
    t = _mm_rsqrt_ss (b);
    _mm_store_ss (&res, t);
    _mm_setcsr (old_mxcsr);
    return res;

}

Une implémentation simple en double-précision basée sur le FMA de hypot() ressemble à ça :

uint64_t __double_as_uint64 (double a)
{
    uint64_t r;
    memcpy (&r, &a, sizeof r);
    return r;
}

double __uint64_as_double (uint64_t a)
{
    double r;
    memcpy (&r, &a, sizeof r);
    return r;
}

double my_hypot (double a, double b)
{
    double fa, fb, mn, mx, scale_in, scale_out, r, s, t;
    uint64_t expo;

    fa = fabs (a);
    fb = fabs (b);
    mx = fmax (fa, fb);
    mn = fmin (fa, fb);

    /* compute scale factors */
    expo = __double_as_uint64 (mx) & 0xff80000000000000ULL;
    scale_in = __uint64_as_double (0x7fc0000000000000ULL - expo);
    scale_out = __uint64_as_double (expo + 0x0020000000000000ULL);

    /* scale mx towards unity */
    mn = mn * scale_in;
    mx = mx * scale_in;

    /* mx in [2**-52, 2**6) */
    s = fma (mx, mx, mn * mn); // 0.75 ulp
    r = sqrt (s);

    /* reverse previous scaling */
    r = r * scale_out;

    /* check for special cases: infinity and NaN */
    t = a + b;
    if (!(fabs (t) <= INFINITY)) r = t; // isnan(t)
    if (mx == INFINITY) r = mx; // isinf(mx)

    return r;
}

Answer 5

http://steve.hollasch.net/cgindex/math/pythag-Root.txt

suggère que l'implémentation plus rapide utilisant sqrt() est quadratique contre cubique pour Moler & Morrison, avec à peu près les mêmes caractéristiques de débordement.