Implementación de sqrt() conforme a IEEE 754 para tipo doble
Estoy tratando de implementar la double __ieee754_sqrt(double x)función que usa instrucciones de hardware para obtener la primera aproximación:
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;
}
Sin embargo, la prueba paranoia.c ( link , link ) se queja:
Square root is neither chopped nor correctly rounded.
Observed errors run from -6.0493828e-01 to 5.0000000e-01 ulps.
Pregunta: ¿cómo implementar lógica adicional para chopping and correct rounding?
UPD. El hardware no es compatible de forma nativa sqrt(). El hardware solo admite la obtención del recíproco de la raíz cuadrada (precisión de 6,75 bits).
UPD2.
- Usé la solución de njuffa (¡muchas gracias!) con cambios menores: usar
qseeddf()en lugar deqseedf()=> usarfma()en lugar defmaf(). ¿Por qué? Porque omitedouble<=>floatconversiones y por lo tanto más rápido. - Sí, el hardware admite instrucciones fusionadas de multiplicación y suma (FMA).
- ¡Gracias a todos por participar en la discusión y por las respuestas detalladas!
- Para todos los interesados en el tema, aquí está la lista de
sqrt()implementaciones:- De las matemáticas de Cygwin. biblioteca (
libm):cygwin-snapshot-20200710-1/newlib/libm/math/e_sqrt.c: con derechos de autorCopyright (C) 1993 by Sun Microsystems. - De la biblioteca GNU C (
glibc):glibc-2.31/sysdeps/ieee754/dbl-64/e_sqrt.c: tituladoIBM Accurate Mathematical Library.glibc-2.31/sysdeps/powerpc/fpu/e_sqrt.c: usando__builtin_fma()funciones.
- De las matemáticas de Cygwin. biblioteca (
Respuestas
Antes de embarcarse en la construcción de una implementación propia, es recomendable buscar en Internet para verificar si se dispone de un código fuente abierto adecuado y bien probado.
Los algoritmos iterativos comunes usan iteraciones sin división para la raíz cuadrada recíproca con la precisión deseada, multiplican hacia atrás con el argumento para calcular la raíz cuadrada y finalmente redondean usando el modo de redondeo deseado. Las iteraciones para la raíz cuadrada recíproca pueden usar iteraciones de Newton-Raphson con convergencia cuadrática (aproximadamente duplicando el número de bits correctos) o iteraciones de Halley con convergencia cúbica (aproximadamente triplicando el número de bits correctos). Si bien existen iteraciones de orden superior, normalmente no se utilizan.
Para mantener el código simple, es recomendable reducir el argumento a un solo intervalo estrecho que comprenda dos binadas consecutivas en el caso de la aritmética binaria de coma flotante. Tenga en cuenta que esto generalmente no da como resultado la implementación de mayor rendimiento debido a la necesidad de manipulación de exponentes. Por razones de rendimiento, las iteraciones iniciales para una implementación de precisión doble a menudo se realizan en precisión simple.
En el ejemplo de implementación de ISO-C99 a continuación, muestro cómo se puede implementar una raíz cuadrada de doble precisión redondeada correctamente a lo largo de esas líneas. Supongo que doublese asigna a IEEE-754 binary64y que floatse asigna a IEEE-754 binary32. Estoy restringiendo a un modo sqrtimplementado con IEEE-754 redondo al más cercano o par.
Muy importante , asumo que el hardware del proceso proporciona instrucciones combinadas de suma y multiplicación y que estas se exponen correctamente a través de las funciones estándar de la biblioteca matemática fmafy fma. En los comentarios, pedí una aclaración de OP sobre la disponibilidad de FMA, pero decidí comenzar con el código antes de que los comentarios estuvieran disponibles. Las implementaciones sin FMA son posibles, pero mucho más desafiantes, y un tratamiento lo suficientemente completo probablemente exceda el espacio de una respuesta de Stackoverflow.
Dado que OP no especificó la arquitectura de destino ni proporcionó detalles de la aproximación inicial, estoy usando mi propia aproximación inicial a continuación basada en una aproximación minimax polinomial en el intervalo [0.25, 1] al que se reducen todos los argumentos no excepcionales. qseedf()los resultados tienen una precisión de aproximadamente 7 bits, por lo que es un poco mejor que el hardware integrado de OP. No puedo evaluar si esta diferencia es significativa.
El algoritmo, en particular la lógica de redondeo, se basa en las ideas de Peter Markstein, por lo que estoy razonablemente seguro de que el algoritmo es correcto por construcción. He implementado solo pruebas muy rudimentarias aquí. La mejor práctica de la industria es probar matemáticamente la corrección de dichos algoritmos; consulte las publicaciones de David Russinoff y John Harrison, por ejemplo. En un apuro, uno podría salirse con la suya con una prueba exhaustiva en dos binadas consecutivas (factible en estos días con un pequeño grupo funcionando durante unos días), junto con pruebas aleatorias y basadas en patrones que ejercitan todas las binadas.
#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 salida del programa anterior debería ser similar a esto:
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 */
...
Esto puede ser más rápido.
Y detener una iteración antes (creo).
Cuando te detengas, compara z*zy x. El z*zserá (creo) no más pequeño que x. Reste 1ulp de zy verifique z*zvs x. No es una comprobación perfecta del "redondeo correcto", pero puede ser "suficientemente bueno" para decidir entre zy z - 1ulp.
Dado que obtuvo una gama tan amplia de errores, me preocupa que el resto del 'hardware' de punto flotante sea descuidado en lo que respecta al redondeo, o incluso a la precisión.
Vaya, lo olvidé. Hubo una razón para darle una aproximación a 1/z-- Continúe aproximando 1/z; puede hacerlo con multiplicaciones en lugar de divisiones, siendo así (en la mayoría del hardware) significativamente más rápido y posiblemente con menos redondeo.
z = ( z + x * z) * 0.5; /* 1st Newton-Raphson iteration */
...
z = 1 / z;
Además, vea si hay una manera de disminuir el exponente en lugar de hacer una multiplicación por / 2.