IEEE 754-konforme sqrt()-Implementierung für Double-Typ
Ich versuche, die double __ieee754_sqrt(double x)Funktion zu implementieren, die Hardwareanweisungen verwendet, um die 1. Annäherung zu erhalten:
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;
}
Der Test von paranoia.c ( link , link ) beschwert sich jedoch:
Square root is neither chopped nor correctly rounded.
Observed errors run from -6.0493828e-01 to 5.0000000e-01 ulps.
Frage: Wie implementiert man zusätzliche Logik für chopping and correct rounding?
UPD. Die Hardware unterstützt nicht nativ sqrt(). Die Hardware unterstützt nur das Erhalten des Kehrwertes der Quadratwurzel (6,75 Bit Genauigkeit).
UPD2.
- Benutzte die Lösung von njuffa (vielen Dank!) mit geringfügigen Änderungen: benutze
qseeddf()stattqseedf()=> benutzefma()stattfmaf(). Wieso den? Weil esdouble<=>floatUmrechnungen weglässt und somit schneller ist. - Ja, Fused Multiply-Add Instruktionen (FMA) werden von der Hardware unterstützt.
- Danke an alle für die Teilnahme an der Diskussion und für die ausführlichen Antworten!
- Für alle, die sich für das Thema interessieren, hier ist die Liste der
sqrt()Implementierungen:- Von Cygwin math. Bibliothek (
libm):cygwin-snapshot-20200710-1/newlib/libm/math/e_sqrt.c: urheberrechtlich geschütztCopyright (C) 1993 by Sun Microsystems. - Aus der GNU C-Bibliothek (
glibc):glibc-2.31/sysdeps/ieee754/dbl-64/e_sqrt.c: berechtigtIBM Accurate Mathematical Library.glibc-2.31/sysdeps/powerpc/fpu/e_sqrt.c:__builtin_fma()Funktionen verwenden.
- Von Cygwin math. Bibliothek (
Antworten
Bevor man mit dem Bau einer eigenen Implementierung beginnt, empfiehlt es sich, im Internet zu recherchieren, ob geeigneter und gut getesteter Open-Source-Code vorhanden ist.
Übliche iterative Algorithmen verwenden teilungsfreie Iterationen für die reziproke Quadratwurzel mit der gewünschten Genauigkeit, multiplizieren mit dem Argument zur Berechnung der Quadratwurzel zurück und runden schließlich unter Verwendung des gewünschten Rundungsmodus. Iterationen für die reziproke Quadratwurzel können entweder Newton-Raphson-Iterationen mit quadratischer Konvergenz (ungefähr Verdoppelung der Anzahl korrekter Bits) oder Halley-Iterationen mit kubischer Konvergenz (ungefähr Verdreifachung der Anzahl korrekter Bits) verwenden. Während Iterationen höherer Ordnung existieren, werden sie typischerweise nicht verwendet.
Um den Code einfach zu halten, ist es ratsam, das Argument im Fall einer binären Gleitkommaarithmetik auf ein einziges schmales Intervall zu reduzieren, das aus zwei aufeinanderfolgenden Binaden besteht. Beachten Sie, dass dies aufgrund der Notwendigkeit einer Exponentenmanipulation im Allgemeinen nicht zu der Implementierung mit der höchsten Leistung führt. Aus Leistungsgründen werden die anfängliche(n) Iteration(en) für eine Implementierung mit doppelter Genauigkeit häufig mit einfacher Genauigkeit durchgeführt.
In der beispielhaften ISO-C99-Implementierung unten zeige ich, wie eine korrekt gerundete Quadratwurzel mit doppelter Genauigkeit in diesem Sinne implementiert werden kann. Ich gehe davon aus, dass das doubleauf IEEE-754 binary64und das floatauf IEEE-754 abgebildet wird binary32. Ich beschränke mich auf einen sqrtmit IEEE-754 implementierten Round-to-Nearest-or-even-Modus.
Sehr wichtig gehe ich davon aus, dass die Prozesshardware verschmolzene Multiplikations-Additions-Anweisungen bereitstellt und dass diese korrekt über die Standard-Mathematikbibliotheksfunktionen und exponiert fmafwerden fma. In Kommentaren hatte ich OP um Klärung der Verfügbarkeit von FMA gebeten, aber beschlossen, mit dem Code zu beginnen, bevor Feedback verfügbar war. Implementierungen ohne FMA sind möglich, aber viel schwieriger, und eine ausreichend vollständige Behandlung würde wahrscheinlich den Raum einer Stackoverflow-Antwort überschreiten.
Da OP die Zielarchitektur nicht spezifiziert oder Details zur Startnäherung angegeben hat, verwende ich unten meine eigene Startnäherung basierend auf einer polynomischen Minimax-Näherung für das Intervall [0,25, 1], auf das alle nicht außergewöhnlichen Argumente reduziert werden. qseedf()Die Ergebnisse sind auf etwa 7 Bit genau, also etwas besser als die integrierte Hardware von OP. Ob dieser Unterschied signifikant ist, kann ich nicht beurteilen.
Der Algorithmus, insbesondere die Rundungslogik, beruht auf den Ideen von Peter Markstein, daher bin ich ziemlich zuversichtlich, dass der Algorithmus konstruktionsbedingt korrekt ist. Ich habe hier nur sehr rudimentäre Tests implementiert. Die beste Industriepraxis besteht darin , die Korrektheit solcher Algorithmen mathematisch zu beweisen , siehe beispielsweise Veröffentlichungen von David Russinoff und John Harrison. Zur Not könnte man mit einem erschöpfenden Test über zwei aufeinanderfolgende Binaden (heutzutage mit einem kleinen Cluster, der einige Tage läuft) davonkommen, gekoppelt mit zufälligen und musterbasierten Tests, die alle Binaden testen.
#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;
}
Die Ausgabe des obigen Programms sollte ungefähr so aussehen:
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 */
...
Dies kann schneller sein.
Und stoppen Sie eine Iteration früher (glaube ich).
Wenn Sie aufhören, vergleichen Sie z*zund x. Der z*zwird (glaube ich) nicht kleiner als x. Subtrahiere 1ulp von zund überprüfe z*zvs x. Es ist keine perfekte Überprüfung der "korrekten Rundung", aber es kann "gut genug" sein, um zwischen zund zu entscheiden z - 1ulp.
Da Sie so viele Fehler haben, mache ich mir Sorgen, dass der Rest der Gleitkomma-'Hardware' beim Runden oder sogar bei der Genauigkeit schlampig ist.
Ups, ich vergaß. Es gab einen Grund, Ihnen eine Annäherung an 1/z-- Fahren Sie mit der Annäherung an 1/z fort; Sie können es mit Multiplikationen anstelle von Divisionen tun, wodurch Sie (in der meisten Hardware) erheblich schneller und möglicherweise mit weniger Rundungen sind.
z = ( z + x * z) * 0.5; /* 1st Newton-Raphson iteration */
...
z = 1 / z;
Sehen Sie auch nach, ob es eine Möglichkeit gibt, den Exponenten zu dekrementieren, anstatt eine Multiplikation für durchzuführen / 2.