Moyen numérique stable de calculer SQRT ((b² * c²) / (1-c²)) pour c dans [-1, 1]

La partie la plus intéressante de cette stabilité est le dénominateur, sqrt (1 - c * c) . Pour cela, tout ce que vous avez à faire est de l'étendre sous le nom de sqrt (1 - c) * sqrt (1 + c) . Je ne pense pas que cela se qualifie vraiment de "truc intelligent", mais c'est tout ce qui est nécessaire.

pour un format de point flottant binaire typique (par exemple IEEE 754 Binary64, mais d'autres formats courants devraient se comporter également bien, avec les Exception possible de choses désagréables comme les Double-double Format), Si c est proche de 1 alors 1 - c sera calculé exactement, par lemme de sterbenz , tandis que 1 + c n'a aucun problème de stabilité. De même, si c est proche de -1 alors 1 + c sera calculé exactement, et 1 - C sera calculé avec précision. Les opérations de racine carrée et de multiplication n'introduiront pas une nouvelle erreur significative.

Voici une démonstration numérique, en utilisant Python sur une machine avec IEEE 754 Binary64 Floating Point et un SQRT correctement articulé opération.

Prenons un c proche de (mais plus petit que) 1 :

>>> better = sqrt(1 - c) * sqrt(1 + c)
>>> better
1.4142136208440158e-05
>>> float((Decimal(better) - good) / Decimal(ulp(float(good))))
-0.7170147200803595

Nous devons être un peu prudents ici: notez que la valeur décimale affichée, 0.999999999 , est une approximation à la valeur exacte de c . La valeur exacte est comme indiqué dans la construction de la chaîne hexadécimale, ou sous forme de fraction, 562949953365017/562949953421312 , et c'est cette valeur exacte pour laquelle nous nous soucions d'obtenir de bons résultats pour.

La valeur exacte de l'expression sqrt (1 - c * c) , arrondi à 100 décimales après le point, est:

>>> from math import ulp
>>> float((Decimal(naive) - good) / Decimal(ulp(float(good))))
208701.28298527992

i a calculé ceci en utilisant le décimal vérifié le résultat en utilisant pari / gp . Voici le calcul de Python:

>>> from math import sqrt
>>> naive = sqrt(1 - c*c)
>>> naive
1.4142136208793713e-05

Si nous calculons naïvement, nous obtenons ce résultat:

>>> from decimal import Decimal, getcontext
>>> getcontext().prec = 1000
>>> good = (1 - Decimal(c) * Decimal(c)).sqrt().quantize(Decimal("1e-100"))
>>> print(good)
0.0000141421362084401590649378320134409069878639187055610216016949959890888003204161068184484972504813

Nous pouvons facilement calculer le Nombre approximatif d'erreur ULPS (avec des excuses pour la quantité de conversion de type en cours - float et Decimal ne peut pas être mélangé directement dans les opérations arithmétiques):

0.0000141421362084401590649378320134409069878639187055610216016949959890888003204161068184484972504813

Ainsi, le résultat naïf est sorti de quelques centaines de mille ulps - à peu près, nous avons perdu environ 5 décimales de précision.

Essayons maintenant avec La version élargie:

>>> c = float.fromhex('0x1.ffffffff24190p-1')
>>> c
0.9999999999

Alors ici, nous sommes exacts à une erreur ULP meilleure. Pas parfaitement correctement arrondi, mais la meilleure chose suivante.

Avec un peu plus de travail, il devrait être possible d'indiquer et de prouver une limite supérieure absolue sur le nombre d'erreur ULPS dans l'expression sqrt ( 1 - C) * SQRT (1 + C) , sur le domaine -1 , en supposant le mode d'arrondi à point flottant binaire IEEE 754 , et les opérations correctement strictes partout. Je ne l'ai pas fait, mais je serais très surpris si cette limite supérieure s'est avérée être plus de 10 ulps.

Mark Dickinson fournit un bon Réponse Pour le cas général, j'ajouterai à cela avec une approche quelque peu plus spécialisée.

De nombreux environnements informatiques de nos jours fournissent une opération appelée multiplié par multiplié, ou FMA pour Short, qui a été spécialement conçu avec des situations comme celle-ci à l'esprit. Dans le calcul de fma (a, b, c) , le produit complet a * b (non tronqué et non lié) entre dans l'addition avec c , alors un seul arrondissement est appliqué à la fin.

expédient actuellement des GPU et des CPU, y compris ceux basés sur les architectures ARM64, X86-64, et Power, incluent généralement une implémentation matérielle rapide de FMA, qui est exposés dans les langages de programmation des familles C et C ++ ainsi que de nombreuses autres en tant que fonction mathématique standard fma () . Certains environnements logiciels - généralement plus anciens, utilisent une émulation de logiciels de FMA, et certaines de ces émulations se sont révélées défectueuses. De plus, de telles émulations ont tendance à être assez lentes.

où la FMA est disponible, l'expression peut être évaluée numériquement stable et sans risque de débordement prématuré et sous-flux comme Fabs (b * c) / sqrt (FMA (C, -C, 1.0)) , où Fabs () est l'opération de valeur absolue pour les opérandes à virgule flottante et sqrt () calcule le calcul de la manière racine carrée. Certains environnements offrent également une opération de racine carrée réciproque, souvent appelée rsqrt () , 1.0)) . L'utilisation de rsqrt () évite la division relativement coûteuse et est donc généralement plus rapide. Cependant, de nombreuses implémentations de rsqrt () ne sont pas correctement arrondies comme sqrt () , donc la précision peut être un peu pire.

Une expérience rapide avec la Le code ci-dessous semble indiquer que l'erreur maximale de l'expression basée sur FMA est d'environ 3 ULPS, tant que b est un numéro de point flottant normal . Je souligne que cela ne prouve pas une erreur liée à une erreur. Le Herbie Tool , qui essaie de trouver des réécritures numériquement avantageuses d'une expression de point flottante donné . Cela semble cependant être un résultat faux, car je ne peux ni penser à aucun avantage particulier ni en trouver un expérimental.

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <math.h>

#define USE_ORIGINAL  (0)
#define USE_HERBIE    (1)

/* function under test */
float func (float b, float c)
{
#if USE_HERBIE
     return fabsf (b * c) * sqrtf (1.0f / fmaf (c, -c, 1.0f));
#else USE_HERBIE
     return fabsf (b * c) / sqrtf (fmaf (c, -c, 1.0f));
#endif // USE_HERBIE
}

/* reference */
double funcd (double b, double c)
{
#if USE_ORIGINAL
    double b2 = b * b;
    double c2 = c * c;
    return sqrt ((b2 * c2) / (1.0 - c2));
#else
    return fabs (b * c) / sqrt (fma (c, -c, 1.0));
#endif
}

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;
}

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

double floatUlpErr (float res, double ref)
{
    uint64_t i, j, err, refi;
    int expoRef;
    
    /* ulp error cannot be computed if either operand is NaN, infinity, zero */
    if (isnan (res) || isnan (ref) || isinf (res) || isinf (ref) ||
        (res == 0.0f) || (ref == 0.0f)) {
        return 0.0;
    }
    /* Convert the float result to an "extended float". This is like a float
       with 56 instead of 24 effective mantissa bits.
    */
    i = ((uint64_t)float_as_uint32(res)) << 32;
    /* Convert the double reference to an "extended float". If the reference is
       >= 2^129, we need to clamp to the maximum "extended float". If reference
       is < 2^-126, we need to denormalize because of the float types's limited
       exponent range.
    */
    refi = double_as_uint64(ref);
    expoRef = (int)(((refi >> 52) & 0x7ff) - 1023);
    if (expoRef >= 129) {
        j = 0x7fffffffffffffffULL;
    } else if (expoRef < -126) {
        j = ((refi << 11) | 0x8000000000000000ULL) >> 8;
        j = j >> (-(expoRef + 126));
    } else {
        j = ((refi << 11) & 0x7fffffffffffffffULL) >> 8;
        j = j | ((uint64_t)(expoRef + 127) << 55);
    }
    j = j | (refi & 0x8000000000000000ULL);
    err = (i < j) ? (j - i) : (i - j);
    return err / 4294967296.0;
}

// Fixes via: Greg Rose, KISS: A Bit Too Simple. http://eprint.iacr.org/2011/007
static unsigned int z=362436069,w=521288629,jsr=362436069,jcong=123456789;
#define znew (z=36969*(z&0xffff)+(z>>16))
#define wnew (w=18000*(w&0xffff)+(w>>16))
#define MWC  ((znew<<16)+wnew)
#define SHR3 (jsr^=(jsr<<13),jsr^=(jsr>>17),jsr^=(jsr<<5)) /* 2^32-1 */
#define CONG (jcong=69069*jcong+13579)                     /* 2^32 */
#define KISS ((MWC^CONG)+SHR3)

#define N  (20)

int main (void)
{
    float b, c, errloc_b, errloc_c, res;
    double ref, err, maxerr = 0;
    
    c = -1.0f;
    while (c <= 1.0f) {
        /* try N random values of `b` per every value of `c` */
        for (int i = 0; i < N; i++) {
            /* allow only normals */
            do {
                b = uint32_as_float (KISS);
            } while (!isnormal (b));
            res = func (b, c);
            ref = funcd ((double)b, (double)c);
            err = floatUlpErr (res, ref);
            if (err > maxerr) {
                maxerr = err;
                errloc_b = b;
                errloc_c = c;
            }
        }
        c = nextafterf (c, INFINITY);
    }
#if USE_HERBIE
    printf ("HERBIE max ulp err = %.5f @ (b=% 15.8e c=% 15.8e)\n", maxerr, errloc_b, errloc_c);
#else // USE_HERBIE
    printf ("SIMPLE max ulp err = %.5f @ (b=% 15.8e c=% 15.8e)\n", maxerr, errloc_b, errloc_c);
#endif // USE_HERBIE
    
    return EXIT_SUCCESS;
}