การใช้งาน sqrt () ตามมาตรฐาน IEEE 754 สำหรับประเภทคู่
ฉันกำลังพยายามใช้double __ieee754_sqrt(double x)ฟังก์ชันที่ใช้คำแนะนำฮาร์ดแวร์เพื่อรับการประมาณครั้งที่ 1:
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;
}
อย่างไรก็ตาม paranoia.c ( การเชื่อมโยง , การเชื่อมโยง ) การทดสอบบ่น:
Square root is neither chopped nor correctly rounded.
Observed errors run from -6.0493828e-01 to 5.0000000e-01 ulps.
คำถาม: จะใช้ตรรกะเพิ่มเติมสำหรับchopping and correct rounding?
UPD. sqrt()ฮาร์ดแวร์ไม่สนับสนุนไม่ได้กำเนิด ฮาร์ดแวร์รองรับเฉพาะการรับส่วนกลับของสแควร์รูทเท่านั้น (ความแม่นยำ 6.75 บิต)
UPD2
- ใช้โซลูชันของ njuffa (ขอบคุณมาก!) กับการเปลี่ยนแปลงเล็กน้อย: ใช้
qseeddf()แทนqseedf()=> ใช้fma()แทนfmaf(). ทำไม? เนื่องจากจะละเว้นdouble<=>floatการแปลงและจึงเร็วกว่า - ใช่คำแนะนำการเพิ่มทวีคูณแบบผสม (FMA) ได้รับการสนับสนุนโดยฮาร์ดแวร์
- ขอบคุณทุกคนที่เข้าร่วมการสนทนาและคำตอบโดยละเอียด!
- สำหรับทุกคนที่สนใจในหัวข้อนี่คือรายการ
sqrt()การใช้งาน:- จาก Cygwin math. ไลบรารี (
libm)cygwin-snapshot-20200710-1/newlib/libm/math/e_sqrt.c:: copyrightedCopyright (C) 1993 by Sun Microsystems. - จากไลบรารี GNU C (
glibc):glibc-2.31/sysdeps/ieee754/dbl-64/e_sqrt.c: ได้รับสิทธิIBM Accurate Mathematical Library.glibc-2.31/sysdeps/powerpc/fpu/e_sqrt.c: ใช้__builtin_fma()ฟังก์ชั่น
- จาก Cygwin math. ไลบรารี (
คำตอบ
ก่อนที่จะเริ่มสร้างการนำไปใช้งานขอแนะนำให้ค้นหาทางอินเทอร์เน็ตเพื่อตรวจสอบว่ามีรหัสโอเพนซอร์สที่เหมาะสมและได้รับการทดสอบอย่างดีหรือไม่
อัลกอริธึมการวนซ้ำทั่วไปใช้การวนซ้ำแบบไม่ต้องหารสำหรับสแควร์รูทซึ่งกันและกันตามความแม่นยำที่ต้องการกลับคูณด้วยอาร์กิวเมนต์เพื่อคำนวณรากที่สองและสุดท้ายปัดโดยใช้โหมดการปัดเศษที่ต้องการ การวนซ้ำสำหรับสแควร์รูทซึ่งกันและกันสามารถใช้การวนซ้ำของนิวตัน - ราฟสันที่มีการบรรจบกันกำลังสอง (โดยประมาณเพิ่มจำนวนบิตที่ถูกต้องเป็นสองเท่า) หรือการวนซ้ำของฮัลเลย์ด้วยการบรรจบกันของลูกบาศก์ (โดยประมาณสามเท่าของจำนวนบิตที่ถูกต้อง) แม้ว่าจะมีการทำซ้ำลำดับที่สูงกว่า แต่โดยทั่วไปจะไม่ใช้
เพื่อให้โค้ดเรียบง่ายขอแนะนำให้ลดอาร์กิวเมนต์เป็นช่วงเวลาแคบ ๆ เดียวซึ่งประกอบด้วยสองไบนารีต่อเนื่องกันในกรณีของเลขคณิตจุดลอยตัวไบนารี โปรดทราบว่าโดยทั่วไปสิ่งนี้ไม่ได้ส่งผลให้เกิดการนำไปใช้งานที่มีประสิทธิภาพสูงสุดเนื่องจากความจำเป็นในการจัดการเลขชี้กำลัง ด้วยเหตุผลด้านประสิทธิภาพการทำซ้ำเริ่มต้นสำหรับการใช้งานที่มีความแม่นยำสองเท่ามักจะดำเนินการด้วยความแม่นยำเดียว
ในการใช้งาน ISO-C99 ที่เป็นแบบอย่างด้านล่างนี้ฉันกำลังแสดงให้เห็นว่าสแควร์รูทที่มีความแม่นยำสองเท่าที่ปัดเศษอย่างถูกต้องสามารถนำไปใช้ตามเส้นเหล่านั้นได้อย่างไร ฉันกำลังสมมติว่าdoubleแผนที่มาตรฐาน IEEE-754 binary64และแผนที่มาตรฐาน IEEE-754float binary32ฉันกำลัง จำกัด การsqrtใช้งานด้วยโหมด IEEE-754 round-to-closest-or-even
ที่สำคัญมากฉันสมมติว่าฮาร์ดแวร์ของกระบวนการให้คำแนะนำการเพิ่มทวีคูณแบบผสมและสิ่งเหล่านี้ถูกเปิดเผยอย่างถูกต้องผ่านฟังก์ชันไลบรารีคณิตศาสตร์มาตรฐานfmafและfma. ในความคิดเห็นฉันได้ขอคำชี้แจงจาก OP เกี่ยวกับความพร้อมใช้งานของ FMA แต่ตัดสินใจที่จะเริ่มต้นใช้งานโค้ดก่อนที่จะมีข้อเสนอแนะ การใช้งานโดยไม่มี FMA เป็นไปได้ แต่มีความท้าทายกว่ามากและการรักษาที่สมบูรณ์เพียงพอน่าจะเกินช่องว่างของคำตอบ Stackoverflow
เนื่องจาก OP ไม่ได้ระบุสถาปัตยกรรมเป้าหมายหรือให้รายละเอียดของการประมาณค่าเริ่มต้นฉันจึงใช้การประมาณเริ่มต้นของตัวเองด้านล่างตามการประมาณค่าต่ำสุดของพหุนามในช่วงเวลา [0.25, 1] ซึ่งอาร์กิวเมนต์ที่ไม่ใช่ข้อยกเว้นทั้งหมดจะลดลงqseedf()ผลลัพธ์มีความแม่นยำประมาณ 7 บิตดังนั้นดีกว่าฮาร์ดแวร์ในตัวของ OP เล็กน้อย ความแตกต่างนี้มีนัยสำคัญหรือไม่ฉันไม่สามารถประเมินได้
อัลกอริทึมโดยเฉพาะตรรกะการปัดเศษขึ้นอยู่กับความคิดของ Peter Markstein ดังนั้นฉันจึงมั่นใจพอสมควรว่าอัลกอริทึมถูกต้องโดยการสร้าง ฉันได้ใช้การทดสอบขั้นพื้นฐานที่นี่เท่านั้น แนวทางปฏิบัติที่ดีที่สุดในอุตสาหกรรมคือการพิสูจน์ความถูกต้องทางคณิตศาสตร์ของอัลกอริทึมดังกล่าวดูสิ่งพิมพ์ของ David Russinoff และ John Harrison เป็นต้น ในการบีบอัดเราอาจสามารถหลีกเลี่ยงการทดสอบอย่างละเอียดถี่ถ้วนในสอง binades ติดต่อกัน (เป็นไปได้ในวันนี้ที่มีคลัสเตอร์ขนาดเล็กทำงานอยู่สองสามวัน) ควบคู่ไปกับการทดสอบแบบสุ่มและตามรูปแบบที่ใช้ 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;
}
ผลลัพธ์ของโปรแกรมข้างต้นควรมีลักษณะดังนี้:
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 */
...
อาจเร็วกว่านี้
และหยุดการทำซ้ำหนึ่งครั้งเร็วขึ้น (ฉันคิดว่า)
เมื่อคุณหยุดการเปรียบเทียบและz*z จะ (ฉันคิด) ไม่มีมีขนาดเล็กกว่า Subtrace 1ulp จากและตรวจสอบเทียบ มันไม่ได้มีการตรวจสอบที่สมบูรณ์แบบของ "การปัดเศษที่ถูกต้อง" แต่มันอาจจะเป็น "ดีพอ" ที่จะตัดสินใจระหว่างและxz*zxzz*zxzz - 1ulp
เนื่องจากคุณมีข้อผิดพลาดมากมายเช่นนี้ฉันจึงกังวลว่าส่วนที่เหลือของ 'ฮาร์ดแวร์' จุดลอยตัวจะเลอะเทอะเมื่อต้องปัดเศษหรือแม้กระทั่งความแม่นยำ
อ๊ะฉันลืมไปแล้ว มีเหตุผลในการให้ค่าประมาณแก่คุณ1/z- ไปที่ค่าประมาณ 1 / z; คุณสามารถทำได้ด้วยการคูณแทนที่จะหารด้วยเหตุนี้ (ในฮาร์ดแวร์ส่วนใหญ่) จะเร็วขึ้นอย่างมากและอาจมีการปัดเศษน้อยลง
z = ( z + x * z) * 0.5; /* 1st Newton-Raphson iteration */
...
z = 1 / z;
นอกจากนี้ดูว่ามีวิธีลดเลขชี้กำลังแทนการคูณสำหรับ/ 2หรือไม่