IEEE 754 konforman sqrt () implementasi untuk tipe ganda
Saya mencoba menerapkan double __ieee754_sqrt(double x)
fungsi yang menggunakan instruksi perangkat keras untuk mendapatkan perkiraan pertama:
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;
}
Namun, tes paranoia.c ( link , link ) mengeluh:
Square root is neither chopped nor correctly rounded.
Observed errors run from -6.0493828e-01 to 5.0000000e-01 ulps.
Pertanyaan: bagaimana cara menerapkan logika tambahan chopping and correct rounding
?
UPD. Perangkat keras tidak mendukung secara bawaan sqrt()
. Perangkat keras hanya mendukung perolehan kebalikan dari akar kuadrat (akurasi 6,75 bit).
UPD2.
- Menggunakan solusi njuffa (terima kasih banyak!) Dengan perubahan kecil: gunakan
qseeddf()
sebagai gantiqseedf()
=> gunakanfma()
sebagai gantifmaf()
. Mengapa? Karena itu menghilangkandouble<=>float
konversi dan karenanya lebih cepat. - Ya, petunjuk multiply-add fused (FMA) didukung oleh perangkat keras.
- Terima kasih untuk semua yang telah berpartisipasi dalam diskusi dan untuk jawaban yang mendetail!
- Untuk semua yang tertarik dengan topik ini, berikut adalah daftar
sqrt()
penerapannya:- Dari Cygwin math. library (
libm
)cygwin-snapshot-20200710-1/newlib/libm/math/e_sqrt.c
:: berhak ciptaCopyright (C) 1993 by Sun Microsystems
. - Dari perpustakaan GNU C (
glibc
):glibc-2.31/sysdeps/ieee754/dbl-64/e_sqrt.c
: berjudulIBM Accurate Mathematical Library
.glibc-2.31/sysdeps/powerpc/fpu/e_sqrt.c
: menggunakan__builtin_fma()
fungsi.
- Dari Cygwin math. library (
Jawaban
Sebelum memulai konstruksi implementasinya sendiri, disarankan untuk mencari di internet untuk memeriksa apakah tersedia kode sumber terbuka yang sesuai dan teruji dengan baik.
Algoritme iteratif umum menggunakan iterasi bebas pembagian untuk akar kuadrat timbal balik ke akurasi yang diinginkan, perkalian balik dengan argumen untuk menghitung akar kuadrat, dan akhirnya membulatkan menggunakan mode pembulatan yang diinginkan. Iterasi untuk akar kuadrat timbal balik dapat menggunakan iterasi Newton-Raphson dengan konvergensi kuadrat (secara kasar menggandakan jumlah bit yang benar) atau iterasi Halley dengan konvergensi kubik (kira-kira melipatgandakan jumlah bit yang benar). Meskipun ada iterasi tingkat tinggi, mereka biasanya tidak digunakan.
Untuk menjaga agar kode tetap sederhana, disarankan untuk mengurangi argumen menjadi satu interval sempit yang terdiri dari dua binade berurutan dalam kasus aritmatika titik mengambang biner. Perhatikan bahwa ini umumnya tidak menghasilkan implementasi kinerja tertinggi karena diperlukan manipulasi eksponen. Untuk alasan performa, iterasi awal untuk implementasi presisi ganda sering dilakukan dalam presisi tunggal.
Dalam contoh penerapan ISO-C99 di bawah ini, saya menunjukkan bagaimana akar kuadrat presisi ganda yang dibulatkan dengan benar dapat diterapkan di sepanjang garis tersebut. Saya berasumsi bahwa double
memetakan ke IEEE-754 binary64
dan float
memetakan ke IEEE-754 binary32
. Saya membatasi untuk sqrt
diimplementasikan dengan IEEE-754 round-to-nearest-or-even mode.
Sangat penting saya berasumsi bahwa perangkat keras proses menyediakan instruksi multiply-add yang menyatu dan ini diekspos dengan benar melalui fungsi perpustakaan matematika standar fmaf
dan fma
. Dalam komentar saya telah meminta klarifikasi dari OP mengenai ketersediaan FMA, tetapi memutuskan untuk memulai kode sebelum umpan balik tersedia. Penerapan tanpa FMA dimungkinkan tetapi jauh lebih menantang, dan penanganan yang cukup lengkap kemungkinan akan melebihi ruang jawaban Stackoverflow.
Karena OP tidak menentukan arsitektur target atau memberikan rincian perkiraan awal, saya menggunakan perkiraan awal saya sendiri di bawah ini berdasarkan perkiraan minimax polinomial pada interval [0,25, 1] di mana semua argumen non-pengecualian dikurangi. qseedf()
hasilnya akurat hingga sekitar 7 bit, jadi sedikit lebih baik daripada perangkat keras bawaan OP. Apakah perbedaan ini signifikan, saya tidak dapat menilai.
Algoritme, khususnya logika pembulatan, bergantung pada gagasan Peter Markstein, oleh karena itu saya cukup yakin bahwa algoritme tersebut benar berdasarkan konstruksi. Saya hanya menerapkan pengujian yang sangat mendasar di sini. Praktik industri terbaik adalah membuktikan secara matematis kebenaran algoritme tersebut, lihat publikasi oleh David Russinoff dan John Harrison misalnya. Dalam keadaan darurat, seseorang mungkin bisa lolos dengan pengujian menyeluruh di dua binade berturut-turut (mungkin saat ini dengan cluster kecil yang berjalan selama beberapa hari), ditambah dengan pengujian acak dan berbasis pola yang menjalankan semua binade.
#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;
}
Output dari program di atas akan terlihat seperti ini:
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 */
...
Ini mungkin lebih cepat.
Dan hentikan satu iterasi lebih cepat (saya pikir.)
Saat Anda berhenti, bandingkan z*z
dan x
. The z*z
akan (saya pikir) tidak lebih kecil dari x
. Subtrace 1ulp dari z
dan memeriksa z*z
vs x
. Ini bukan pemeriksaan sempurna untuk "pembulatan yang benar", tetapi mungkin "cukup baik" untuk memutuskan antara z
dan z - 1ulp
.
Karena Anda mendapatkan begitu banyak kesalahan, saya khawatir bahwa 'perangkat keras' floating point lainnya ceroboh dalam hal pembulatan, atau bahkan presisi.
Ups, saya lupa. Ada alasan untuk memberi Anda perkiraan untuk 1/z
- Lanjutkan ke perkiraan 1 / z; Anda dapat melakukannya dengan mengalikan alih-alih membagi, sehingga (di sebagian besar perangkat keras) secara signifikan lebih cepat dan mungkin dengan lebih sedikit pembulatan.
z = ( z + x * z) * 0.5; /* 1st Newton-Raphson iteration */
...
z = 1 / z;
Juga, lihat apakah ada cara untuk menurunkan eksponen daripada melakukan perkalian / 2
.