Implémentation sqrt() conforme IEEE 754 pour le type double
J'essaie d'implémenter la double __ieee754_sqrt(double x)fonction qui utilise l'instruction matérielle pour obtenir la 1ère approximation :
double __ieee754_sqrt(double x) {
double z;
/* get reciprocal of the square root (6.75 bits accuracy) */
__asm(" QSEED.DF %0,%1 \n": "=e" (z):"e" (x):);
z = 1 / z;
z = ( z + x / z) / 2; /* 1st Newton-Raphson iteration */
z = ( z + x / z) / 2; /* 2nd Newton-Raphson iteration */
z = ( z + x / z) / 2; /* 3rd Newton-Raphson iteration */
z = ( z + x / z) / 2; /* 4th Newton-Raphson iteration */
return z;
}
Cependant, le test paranoia.c ( lien , lien ) se plaint :
Square root is neither chopped nor correctly rounded.
Observed errors run from -6.0493828e-01 to 5.0000000e-01 ulps.
Question : comment implémenter une logique supplémentaire pourchopping and correct rounding ?
UPD. Le matériel ne prend pas en charge nativement sqrt(). Le matériel ne prend en charge que l'obtention de l'inverse de la racine carrée (précision de 6,75 bits).
UPD2.
- J'ai utilisé la solution de njuffa (merci beaucoup !) avec des modifications mineures : utilisez
qseeddf()au lieu deqseedf()=> utilisezfma()au lieu defmaf(). Pourquoi? Parce qu'il ometdouble<=>floatles conversions et donc plus rapide. - Oui, les instructions fusionnées de multiplication-addition (FMA) sont prises en charge par le matériel.
- Merci à tous d'avoir participé à la discussion et pour les réponses détaillées !
- A tous ceux que le sujet intéresse, voici la liste des
sqrt()implémentations :- D'après les mathématiques de Cygwin. bibliothèque (
libm):cygwin-snapshot-20200710-1/newlib/libm/math/e_sqrt.c: protégé par copyrightCopyright (C) 1993 by Sun Microsystems. - Depuis la bibliothèque GNU C (
glibc):glibc-2.31/sysdeps/ieee754/dbl-64/e_sqrt.c: intituléIBM Accurate Mathematical Library.glibc-2.31/sysdeps/powerpc/fpu/e_sqrt.c: utilisation__builtin_fma()des fonctions.
- D'après les mathématiques de Cygwin. bibliothèque (
Réponses
Avant de se lancer dans la construction de sa propre implémentation, il est conseillé de rechercher sur Internet pour vérifier si un code open-source adapté et bien testé est disponible.
Les algorithmes itératifs courants utilisent des itérations sans division pour la racine carrée réciproque avec la précision souhaitée, rétromultiplient avec l'argument pour calculer la racine carrée, et enfin arrondissent en utilisant le mode d'arrondi souhaité. Les itérations pour la racine carrée réciproque peuvent utiliser soit des itérations de Newton-Raphson avec convergence quadratique (doublant approximativement le nombre de bits corrects), soit des itérations de Halley avec convergence cubique (triplant approximativement le nombre de bits corrects). Bien qu'il existe des itérations d'ordre supérieur, elles ne sont généralement pas utilisées.
Pour garder le code simple, il est conseillé de réduire l'argument à un seul intervalle étroit comprenant deux binades consécutives dans le cas de l'arithmétique binaire à virgule flottante. Notez que cela n'entraîne généralement pas l'implémentation la plus performante en raison de la nécessité d'une manipulation des exposants. Pour des raisons de performances, la ou les itérations initiales d'une implémentation en double précision sont souvent réalisées en simple précision.
Dans l'exemple d'implémentation ISO-C99 ci-dessous, je montre comment une racine carrée double précision correctement arrondie peut être implémentée dans ce sens. Je suppose que cela doublecorrespond à IEEE-754 binary64et que cela floatcorrespond à IEEE-754 binary32. Je me limite à un sqrtmode implémenté avec IEEE-754 round-to-nearest-or-even.
Très important , je suppose que le matériel de processus fournit des instructions de multiplication-addition fusionnées et que celles-ci sont correctement exposées via les fonctions de bibliothèque mathématique standard fmafet fma. Dans les commentaires, j'avais demandé des éclaircissements à OP quant à la disponibilité de FMA, mais j'ai décidé de commencer le code avant que les commentaires ne soient disponibles. Les implémentations sans FMA sont possibles mais beaucoup plus difficiles, et un traitement suffisamment complet dépasserait probablement l'espace d'une réponse Stackoverflow.
Étant donné que OP n'a pas spécifié l'architecture cible ni fourni de détails sur l'approximation de départ, j'utilise ma propre approximation de départ ci-dessous basée sur une approximation polynomiale minimax sur l'intervalle [0,25, 1] auquel tous les arguments non exceptionnels sont réduits. qseedf()les résultats sont précis à environ 7 bits, donc légèrement meilleurs que le matériel intégré d'OP. Je ne peux pas évaluer si cette différence est significative.
L'algorithme, en particulier la logique d'arrondi, repose sur les idées de Peter Markstein, donc je suis raisonnablement confiant que l'algorithme est correct par construction. Je n'ai mis en œuvre que des tests très rudimentaires ici. Les meilleures pratiques de l'industrie consistent à prouver mathématiquement l'exactitude de ces algorithmes, voir les publications de David Russinoff et John Harrison par exemple. À la rigueur, on pourrait peut-être s'en tirer avec un test exhaustif sur deux binades consécutives (réalisable de nos jours avec un petit cluster fonctionnant pendant quelques jours), couplé à des tests aléatoires et basés sur des modèles exerçant toutes les binades.
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <math.h>
/* Approximate 1/sqrt(a) on [0.25, 1] with an accuracy of about 7 bits */
float qseedf (float a)
{
float r;
r = -2.43845296f;
r = fmaf (r, a, 6.22994471f);
r = fmaf (r, a, -5.91090727f);
r = fmaf (r, a, 3.11237526f);
return r;
}
double my_sqrt (double a)
{
const double QNAN_INDEFINITE = 0.0 / 0.0;
const double half = 0.5;
const double three_eighth = 0.375;
double refined_rsqrt_approx, sqrt_approx, sqrt_residual, result, b;
double rsqrt_approx, rsqrt_approx_err, rsqrt_approx_squared, reduced_arg;
float argf, approxf, approxf_err;
int e, t, f;
/* handle normal cases */
if ((a >= 0) && (a < INFINITY)) {
/* compute exponent adjustments */
b = frexp (a, &e);
t = e - 2*512;
f = t / 2;
t = t - 2 * f;
f = f + 512;
/* map argument into the primary approximation interval [0.25,1) */
reduced_arg = ldexp (b, t);
/* Compute initial low-precision approximation */
argf = (float)reduced_arg;
approxf = qseedf (argf);
/* Apply two Newton-Raphson iterations with quadratic convergence */
approxf_err = fmaf (-argf, approxf * approxf, 1.0f);
approxf = fmaf (0.5f * approxf, approxf_err, approxf);
approxf_err = fmaf (-argf, approxf * approxf, 1.0f);
approxf = fmaf (0.5f * approxf, approxf_err, approxf);
/* rsqrt approximation is now accurate to 1 single-precision ulp */
rsqrt_approx = (double)approxf;
/* Perform a Halley iteration wih cubic convergence. Based on the work
of Peter Markstein. See: Peter Markstein, "IA-64 and Elementary
Functions", Prentice Hall 2000
*/
rsqrt_approx_squared = rsqrt_approx * rsqrt_approx;
rsqrt_approx_err = fma (-reduced_arg, rsqrt_approx_squared, 1.0);
refined_rsqrt_approx = fma (fma (rsqrt_approx_err, three_eighth, half),
rsqrt_approx * rsqrt_approx_err, rsqrt_approx);
sqrt_approx = reduced_arg * refined_rsqrt_approx;
sqrt_residual = fma (-sqrt_approx, sqrt_approx, reduced_arg);
result = fma (sqrt_residual, half * refined_rsqrt_approx, sqrt_approx);
/* map back from primary approximation interval by jamming exponent */
result = ldexp (result, f);
} else {
/* handle special cases */
result = (a < 0) ? QNAN_INDEFINITE : (a + a);
}
return result;
}
/*
https://groups.google.com/forum/#!original/comp.lang.c/qFv18ql_WlU/IK8KGZZFJx4J
From: geo <[email protected]>
Newsgroups: sci.math,comp.lang.c,comp.lang.fortran
Subject: 64-bit KISS RNGs
Date: Sat, 28 Feb 2009 04:30:48 -0800 (PST)
This 64-bit KISS RNG has three components, each nearly
good enough to serve alone. The components are:
Multiply-With-Carry (MWC), period (2^121+2^63-1)
Xorshift (XSH), period 2^64-1
Congruential (CNG), period 2^64
*/
static uint64_t kiss64_x = 1234567890987654321ULL;
static uint64_t kiss64_c = 123456123456123456ULL;
static uint64_t kiss64_y = 362436362436362436ULL;
static uint64_t kiss64_z = 1066149217761810ULL;
static uint64_t kiss64_t;
#define MWC64 (kiss64_t = (kiss64_x << 58) + kiss64_c, \
kiss64_c = (kiss64_x >> 6), kiss64_x += kiss64_t, \
kiss64_c += (kiss64_x < kiss64_t), kiss64_x)
#define XSH64 (kiss64_y ^= (kiss64_y << 13), kiss64_y ^= (kiss64_y >> 17), \
kiss64_y ^= (kiss64_y << 43))
#define CNG64 (kiss64_z = 6906969069ULL * kiss64_z + 1234567ULL)
#define KISS64 (MWC64 + XSH64 + CNG64)
int main (void)
{
const uint64_t N = 10000000000ULL; /* desired number of test cases */
double arg, ref, res;
uint64_t argi, refi, resi, count = 0;
double spec[] = {0, 1, INFINITY, NAN};
printf ("test a few special cases:\n");
for (int i = 0; i < sizeof (spec)/sizeof(spec[0]); i++) {
printf ("my_sqrt(%22.13a) = %22.13a\n", spec[i], my_sqrt(spec[i]));
printf ("my_sqrt(%22.13a) = %22.13a\n", -spec[i], my_sqrt(-spec[i]));
}
printf ("test %llu random cases:\n", N);
do {
count++;
argi = KISS64;
memcpy (&arg, &argi, sizeof arg);
res = my_sqrt (arg);
ref = sqrt (arg);
memcpy (&resi, &res, sizeof resi);
memcpy (&refi, &ref, sizeof refi);
if (resi != refi) {
printf ("\rerror @ arg=%22.13a res=%22.13a ref=%22.13a\n",
arg, res, ref);
return EXIT_FAILURE;
}
if ((count & 0xfffff) == 0) printf ("\r[%llu]", count);
} while (count < N);
printf ("\r[%llu]", count);
printf ("\ntests PASSED\n");
return EXIT_SUCCESS;
}
La sortie du programme ci-dessus devrait ressembler à ceci :
test a few special cases:
my_sqrt( 0x0.0000000000000p+0) = 0x0.0000000000000p+0
my_sqrt( -0x0.0000000000000p+0) = -0x0.0000000000000p+0
my_sqrt( 0x1.0000000000000p+0) = 0x1.0000000000000p+0
my_sqrt( -0x1.0000000000000p+0) = -0x1.#IND000000000p+0
my_sqrt( 0x1.#INF000000000p+0) = 0x1.#INF000000000p+0
my_sqrt( -0x1.#INF000000000p+0) = -0x1.#IND000000000p+0
my_sqrt( 0x1.#QNAN00000000p+0) = 0x1.#QNAN00000000p+0
my_sqrt( -0x1.#QNAN00000000p+0) = -0x1.#QNAN00000000p+0
test 10000000000 random cases:
[10000000000]
tests PASSED
z = 1 / z;
z = ( z + x / z) / 2; /* 1st Newton-Raphson iteration */
...
-->
z = 1 / z;
z += ( x / z - z) * 0.5; /* 1st Newton-Raphson iteration */
...
Cela peut être plus rapide.
Et arrêter une itération plus tôt (je pense.)
Lorsque vous vous arrêtez, comparez z*zet x. Le z*zsera (je pense) pas plus petit que x. Soustracez 1 ulp de zet vérifiez z*zvs x. Ce n'est pas une vérification parfaite de "l'arrondi correct", mais cela peut être "assez bon" pour décider entre zet z - 1ulp.
Étant donné que vous avez une telle gamme d'erreurs, je crains que le reste du "matériel" à virgule flottante ne soit bâclé en matière d'arrondi, voire de précision.
Oups, j'ai oublié. Il y avait une raison pour vous donner une approximation de 1/z-- Continuez à approximer 1/z ; vous pouvez le faire avec des multiplications au lieu de divisions, ce qui est (dans la plupart des matériels) beaucoup plus rapide et éventuellement avec moins d'arrondi.
z = ( z + x * z) * 0.5; /* 1st Newton-Raphson iteration */
...
z = 1 / z;
Aussi, voyez s'il existe un moyen de décrémenter l'exposant au lieu de faire une multiplication pour / 2.