双精度类型的 IEEE 754 一致性 sqrt() 实现
IEEE 754 conformant sqrt() implementation for double type
我正在尝试实现 double __ieee754_sqrt(double x) 函数，该函数使用硬件指令获得第一个近似值：
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 (link, link) 测试抱怨：
Square root is neither chopped nor correctly rounded.
Observed errors run from -6.0493828e-01 to 5.0000000e-01 ulps.
问题：如何为chopping and correct rounding实现额外的逻辑？
更新。硬件本身不支持 sqrt()。硬件只支持求平方根的倒数（6.75位精度）
UPD2.
- 使用 njuffa 的解决方案（非常感谢！），稍作改动：使用
qseeddf() 而不是 qseedf() => 使用 fma() 而不是 fmaf()。为什么？因为它省略了 double<=>float 转换，因此速度更快。
- 是的，硬件支持融合 multiply-add 指令 (FMA)。
- 感谢大家参与讨论和详细解答！
- 对于所有对该主题感兴趣的人，这里是
sqrt() 实现的列表：
- 来自 Cygwin 数学。库 (
libm)：cygwin-snapshot-20200710-1/newlib/libm/math/e_sqrt.c：受版权保护 Copyright (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() 函数。
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和x。 z*z 将（我认为）不小于 x。从 z 减去 1ulp 并检查 z*z 与 x。它不是“正确舍入”的完美检查，但它可能“足够好”以在 z 和 z - 1ulp 之间做出决定。
既然你得到了如此大的误差范围，我担心其余的浮点数 'hardware' 在四舍五入甚至精度方面都很草率。
糟糕，我忘记了。给你一个近似值 1/z 是有原因的——继续近似 1/z；你可以用乘法而不是除法来做，从而（在大多数硬件中）明显更快并且可能更少舍入。
z = ( z + x * z) * 0.5; /* 1st Newton-Raphson iteration */
...
z = 1 / z;
此外，看看是否有一种方法可以减少指数而不是对 / 2 进行乘法运算。
在着手构建自己的实现之前，建议在互联网上搜索以检查是否合适并且 well-tested open-source 代码可用。
常用迭代算法使用 division-free 次迭代来计算平方根的倒数以达到所需的精度，back-multiply 使用参数计算平方根，最后使用所需的舍入模式进行舍入。平方根倒数的迭代可以使用具有二次收敛的 Newton-Raphson 迭代（大约使正确位数加倍）或具有三次收敛的哈雷迭代（大约使正确位数增加三倍）。虽然存在 higher-order 次迭代，但通常不使用它们。
为了保持代码简单，在二进制 floating-point 算术的情况下，建议将参数减少为包含两个连续二进制数的单个窄区间。请注意，由于需要指数操作，这通常不会产生最高性能的实现。出于性能原因，double-precision 实现的初始迭代通常以单精度执行。
在下面的示例性 ISO-C99 实施中，我展示了如何按照这些原则实施正确舍入的 double-precision 平方根。我假设 double 映射到 IEEE-754 binary64 并且 float 映射到 IEEE-754 binary32。我限制使用 IEEE-754 round-to-nearest-or-even 模式实现的 sqrt。
非常重要 我假设处理器硬件提供融合 multiply-add 指令，并且这些指令通过标准数学库函数 fmaf 和fma。在评论中，我曾要求 OP 澄清 FMA 的可用性，但决定在获得反馈之前开始编写代码。没有 FMA 的实现是可能的，但更具挑战性，并且足够完整的处理可能会超过 Whosebug 答案的 space。
由于 OP 没有指定目标架构或提供起始近似的详细信息，我在下面使用我自己的起始近似，基于区间 [0.25, 1] 上的多项式极小极大近似，所有 non-exceptional 争论减少了。 qseedf() 结果精确到大约 7 位，因此比 OP 的 built-in 硬件略好。这种差异是否显着，我无法评估。
该算法，尤其是舍入逻辑，依赖于 Peter Markstein 的思想，因此我有理由相信该算法在结构上是正确的。我在这里只实施了非常基本的测试。最佳行业实践是从数学上证明 此类算法的正确性，例如，请参阅 David Russinoff 和 John Harrison 的出版物。在紧要关头，一个人可能能够通过两个连续的 binades 进行详尽的测试（这些天用小集群 运行 几天是可行的），再加上随机和 pattern-based 测试锻炼所有二进制文件。
#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 <gmars...@gmail.com>
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
我正在尝试实现 double __ieee754_sqrt(double x) 函数，该函数使用硬件指令获得第一个近似值：
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 (link, link) 测试抱怨：
Square root is neither chopped nor correctly rounded.
Observed errors run from -6.0493828e-01 to 5.0000000e-01 ulps.
问题：如何为chopping and correct rounding实现额外的逻辑？
更新。硬件本身不支持 sqrt()。硬件只支持求平方根的倒数（6.75位精度）
UPD2.
- 使用 njuffa 的解决方案（非常感谢！），稍作改动：使用
qseeddf()而不是qseedf()=> 使用fma()而不是fmaf()。为什么？因为它省略了double<=>float转换，因此速度更快。 - 是的，硬件支持融合 multiply-add 指令 (FMA)。
- 感谢大家参与讨论和详细解答！
- 对于所有对该主题感兴趣的人，这里是
sqrt()实现的列表：- 来自 Cygwin 数学。库 (
libm)：cygwin-snapshot-20200710-1/newlib/libm/math/e_sqrt.c：受版权保护Copyright (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 数学。库 (
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和x。 z*z 将（我认为）不小于 x。从 z 减去 1ulp 并检查 z*z 与 x。它不是“正确舍入”的完美检查，但它可能“足够好”以在 z 和 z - 1ulp 之间做出决定。
既然你得到了如此大的误差范围，我担心其余的浮点数 'hardware' 在四舍五入甚至精度方面都很草率。
糟糕，我忘记了。给你一个近似值 1/z 是有原因的——继续近似 1/z；你可以用乘法而不是除法来做，从而（在大多数硬件中）明显更快并且可能更少舍入。
z = ( z + x * z) * 0.5; /* 1st Newton-Raphson iteration */
...
z = 1 / z;
此外，看看是否有一种方法可以减少指数而不是对 / 2 进行乘法运算。
在着手构建自己的实现之前，建议在互联网上搜索以检查是否合适并且 well-tested open-source 代码可用。
常用迭代算法使用 division-free 次迭代来计算平方根的倒数以达到所需的精度，back-multiply 使用参数计算平方根，最后使用所需的舍入模式进行舍入。平方根倒数的迭代可以使用具有二次收敛的 Newton-Raphson 迭代（大约使正确位数加倍）或具有三次收敛的哈雷迭代（大约使正确位数增加三倍）。虽然存在 higher-order 次迭代，但通常不使用它们。
为了保持代码简单，在二进制 floating-point 算术的情况下，建议将参数减少为包含两个连续二进制数的单个窄区间。请注意，由于需要指数操作，这通常不会产生最高性能的实现。出于性能原因，double-precision 实现的初始迭代通常以单精度执行。
在下面的示例性 ISO-C99 实施中，我展示了如何按照这些原则实施正确舍入的 double-precision 平方根。我假设 double 映射到 IEEE-754 binary64 并且 float 映射到 IEEE-754 binary32。我限制使用 IEEE-754 round-to-nearest-or-even 模式实现的 sqrt。
非常重要 我假设处理器硬件提供融合 multiply-add 指令，并且这些指令通过标准数学库函数 fmaf 和fma。在评论中，我曾要求 OP 澄清 FMA 的可用性，但决定在获得反馈之前开始编写代码。没有 FMA 的实现是可能的，但更具挑战性，并且足够完整的处理可能会超过 Whosebug 答案的 space。
由于 OP 没有指定目标架构或提供起始近似的详细信息，我在下面使用我自己的起始近似，基于区间 [0.25, 1] 上的多项式极小极大近似，所有 non-exceptional 争论减少了。 qseedf() 结果精确到大约 7 位，因此比 OP 的 built-in 硬件略好。这种差异是否显着，我无法评估。
该算法，尤其是舍入逻辑，依赖于 Peter Markstein 的思想，因此我有理由相信该算法在结构上是正确的。我在这里只实施了非常基本的测试。最佳行业实践是从数学上证明 此类算法的正确性，例如，请参阅 David Russinoff 和 John Harrison 的出版物。在紧要关头，一个人可能能够通过两个连续的 binades 进行详尽的测试（这些天用小集群 运行 几天是可行的），再加上随机和 pattern-based 测试锻炼所有二进制文件。
#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 <gmars...@gmail.com>
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