Çift tip için IEEE 754 uyumlu sqrt () uygulaması
double __ieee754_sqrt(double x)
1. yaklaşımı elde etmek için donanım talimatını kullanan işlevi uygulamaya çalışıyorum :
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;
}
Bununla birlikte, paranoia.c ( bağlantı , bağlantı ) testi şikayet ediyor:
Square root is neither chopped nor correctly rounded.
Observed errors run from -6.0493828e-01 to 5.0000000e-01 ulps.
Soru: Ek mantık nasıl uygulanır chopping and correct rounding
?
UPD. Donanım yerel olarak desteklemiyor sqrt()
. Donanım, yalnızca karekökün tersinin elde edilmesini destekler (6.75 bit doğruluk).
UPD2.
- Küçük değişikliklerle njuffa'nın çözümü (çok teşekkürler!) Kullanıldı:
qseeddf()
yerineqseedf()
=> kullanınfma()
yerine kullanınfmaf()
. Neden? Çünküdouble<=>float
dönüşümleri atlar ve dolayısıyla daha hızlıdır. - Evet, sigortalı çoklu ekleme talimatları (FMA) donanım tarafından desteklenir.
- Tartışmaya katıldığınız ve ayrıntılı cevaplar için herkese teşekkürler!
- Konuyla ilgilenen herkes için, işte
sqrt()
uygulamaların listesi :- Cygwin Math'dan. kütüphane (
libm
)cygwin-snapshot-20200710-1/newlib/libm/math/e_sqrt.c
:: telif hakkıyla korunmaktadırCopyright (C) 1993 by Sun Microsystems
. - GNU C kitaplığından (
glibc
):glibc-2.31/sysdeps/ieee754/dbl-64/e_sqrt.c
: başlıklıIBM Accurate Mathematical Library
.glibc-2.31/sysdeps/powerpc/fpu/e_sqrt.c
:__builtin_fma()
işlevleri kullanma .
- Cygwin Math'dan. kütüphane (
Yanıtlar
Kendi uygulamasının inşasına başlamadan önce, uygun ve iyi test edilmiş açık kaynak kodunun mevcut olup olmadığını kontrol etmek için internette arama yapmanız önerilir.
Yaygın yinelemeli algoritmalar, karşılıklı karekök için istenen doğrulukta bölmesiz yinelemeler kullanır, karekökü hesaplamak için bağımsız değişkenle geri çarpma ve son olarak istenen yuvarlama modunu kullanarak yuvarlama. Karşılıklı karekök için yinelemeler, ya kuadratik yakınsamalı Newton-Raphson yinelemelerini (doğru bit sayısını kabaca ikiye katlayarak) ya da kübik yakınsamalı Halley yinelemelerini (doğru bit sayısını kabaca üçe katlayarak) kullanabilir. Daha yüksek dereceli yinelemeler var olsa da, genellikle kullanılmazlar.
Kodu basit tutmak için, ikili kayan nokta aritmetiği durumunda argümanı iki ardışık ikili içeren tek bir dar aralığa indirgemek tavsiye edilir. Bunun, üs değiştirme ihtiyacından dolayı genellikle en yüksek performans uygulamasına neden olmadığına dikkat edin. Performans nedenleriyle, çift duyarlıklı bir uygulama için ilk yinelemeler genellikle tek duyarlıkta gerçekleştirilir.
Aşağıdaki örnek ISO-C99 uygulamasında, bu çizgiler boyunca doğru şekilde yuvarlatılmış çift kesinlikli bir karekökün nasıl uygulanabileceğini gösteriyorum. Bunu varsayıyorum double
IEEE-754 haritalar binary64
ve float
IEEE-754 için haritalar binary32
. sqrt
IEEE-754 yuvarlaktan en yakına veya eşit modda uygulanmış bir modla kısıtlıyorum .
Çok önemli bir şekilde, işlem donanımının kaynaşmış çoklu ekleme komutları sağladığını ve bunların standart matematik kütüphanesi işlevleri fmaf
ve fma
. Yorumlarda, OP'den FMA'nın kullanılabilirliği konusunda açıklama istemiştim, ancak geri bildirim sağlanmadan önce kodu kullanmaya başlamaya karar verdim. FMA'sız uygulamalar mümkündür ancak çok daha zordur ve yeterince eksiksiz bir işlem muhtemelen bir Stackoverflow yanıtının alanını aşacaktır.
OP, hedef mimariyi belirtmediğinden veya başlangıç yaklaşımının ayrıntılarını sağlamadığından, istisnai olmayan tüm argümanların indirgendiği [0.25, 1] aralığına ilişkin bir polinom minimaks yaklaşımına dayanan aşağıdaki kendi başlangıç yaklaşımımı kullanıyorum. qseedf()
sonuçlar yaklaşık 7 bit kadar doğrudur, bu nedenle OP'nin yerleşik donanımından biraz daha iyidir. Bu farkın önemli olup olmadığını değerlendiremiyorum.
Algoritma, özellikle yuvarlama mantığı, Peter Markstein'ın fikirlerine dayanır, bu nedenle algoritmanın yapısal olarak doğru olduğundan makul ölçüde eminim. Burada sadece çok temel testler uyguladım. En iyi endüstri uygulaması, bu tür algoritmaların doğruluğunu matematiksel olarak kanıtlamaktır , örneğin David Russinoff ve John Harrison'ın yayınlarına bakınız. Bir tutamda, iki ardışık ikili boyunca kapsamlı bir testle (bu günlerde birkaç gün çalışan küçük bir küme ile uygulanabilir), tüm binadları deneyen rastgele ve desen tabanlı testlerle birleştirilmiş bir şekilde kurtulmak mümkün olabilir.
#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;
}
Yukarıdaki programın çıktısı şuna benzer görünmelidir:
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 */
...
Bu daha hızlı olabilir.
Ve bir yinelemeyi daha erken durdurun (sanırım)
Durduğunuzda, karşılaştırın z*z
ve x
. z*z
(Sanırım) daha küçük olacaktır x
. Dan 1ulp Subtrace z
ve check z*z
vs x
. "Doğru yuvarlamanın" mükemmel bir kontrolü değildir, ancak z
ve arasında karar vermek için "yeterince iyi" olabilir z - 1ulp
.
Bu kadar geniş bir hata aralığına sahip olduğunuz için, kayan nokta 'donanımının' geri kalanının yuvarlama, hatta hassasiyet söz konusu olduğunda özensiz olduğundan endişeleniyorum.
Oops, unuttum. Size bir tahmin vermenin bir nedeni vardı 1/z
- 1 / z'ye yaklaşmaya devam edin; bunu bölmeler yerine çarpmalarla yapabilirsiniz, böylece (çoğu donanımda) önemli ölçüde daha hızlı ve muhtemelen daha az yuvarlama ile yapabilirsiniz.
z = ( z + x * z) * 0.5; /* 1st Newton-Raphson iteration */
...
z = 1 / z;
Ayrıca, için çarpma yapmak yerine üssü azaltmanın bir yolu olup olmadığına bakın / 2
.