如何解决考虑到硬件限制和使用限制的符合 IEEE 754 的 sqrtf() 实现
IEEE 754 conformant sqrt() implementation for double type 的后续问题。
上下文：考虑到以下硬件限制和使用限制，需要实现符合 IEEE 754 的 sqrtf()：
-
提供特殊指令
qseed.f来求平方根倒数的近似值（结果精度不低于6.75位，因此始终在精度结果的±1%以内). -
单精度 FP：
一个。 Support by HW (SP FPU)：有支持；
B. SW（库）支持：有支持；
c.支持次正规数：不支持（
FLT_HAS_SUBNORM为 0）。 -
双精度 FP：
一个。硬件支持（DP FPU）：不支持；
B. SW（库）支持：有支持；
c.支持次正规数：不支持（
DBL_HAS_SUBNORM为 0）。
我找到了 John Harrison 的一个演示文稿并最终实现了这个实现（请注意，这里的 qseed.f 被替换为 rsqrtf()）：
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
// https://github.com/nickzman/hyperspace/blob/master/frsqrt.hh
#if 1
float rsqrtf ( float x )
{
const float xhalf = 0.5f * x;
int i = *(int*) & x;
i = 0x5f375a86 - ( i >> 1 );
x = *(float*) & i;
x = x * ( 1.5f - xhalf * x * x );
x = x * ( 1.5f - xhalf * x * x );
x = x * ( 1.5f - xhalf * x * x );
return x;
}
#else
float rsqrtf ( float x )
{
return 1.0f / sqrtf( x );
}
#endif
float sqrtfr_jh( float x,float r )
{
/*
* John Harrison,Formal Verification Methods 5: Floating Point Verification,* Intel Corporation,12 December 2002,document name: slides5.pdf,page 14,* slide "The square root algorithm".
* URL: https://www.cl.cam.ac.uk/~jrh13/slides/anu-09_12dec02/slides5.pdf
*/
double rd,b,z0,s0,d,k,h0,e,t0,s1,c,d1,h1,s;
static const double half = 0.5;
static const double one = 1.0;
static const double three = 3.0;
static const double two = 2.0;
rd = (double)r;
b = half * x;
z0 = rd * rd;
s0 = x * rd;
d = fma( -b,half );
k = fma( x,rd,-s0 );
h0 = half * rd;
e = fma( three / two,one );
t0 = fma( d,k );
s1 = fma( e,s0 );
c = fma( d,one );
d1 = fma( -s1,x );
h1 = c * h0;
s = fma( d1,s1 );
return (float)s;
}
float my_sqrtf( float x )
{
/* handle special cases */
if (x == 0) {
return x + x;
}
/* handle normal cases */
if ((x > 0) && (x < INFINITY)) {
return sqrtfr_jh( x,rsqrtf( x ) );
}
/* handle special cases */
return (x < 0) ? NAN : (x + x);
}
/*
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 */
float arg,ref,res;
uint64_t argi64;
uint32_t refi,resi;
uint64_t count = 0;
float spec[] = {0.0f,1.0f,INFINITY,NAN};
printf ("test a few special cases:\n");
for (int i = 0; i < sizeof (spec)/sizeof(spec[0]); i++) {
printf ("my_sqrt(%a) = %a\n",spec[i],my_sqrtf(spec[i]));
printf ("my_sqrt(%a) = %a\n",-spec[i],my_sqrtf(-spec[i]));
}
printf ("test %lu random cases:\n",N);
do {
argi64 = KISS64;
memcpy (&arg,&argi64,sizeof arg);
if ( fpclassify(arg) == FP_SUBNORMAL )
{
continue;
}
++count;
res = my_sqrtf (arg);
ref = sqrtf (arg);
memcpy (&resi,&res,sizeof resi);
memcpy (&refi,&ref,sizeof refi);
if ( ! ( isnan(res) && isnan(ref) ) )
if (resi != refi) {
printf ("\rerror @ arg=%a (%e)\n",arg,arg);
printf ("\rerror @ res=%a (%e)\n",res,res);
printf ("\rerror @ ref=%a (%e)\n",ref);
return EXIT_FAILURE;
}
if ((count & 0xfffff) == 0) printf ("\r[%lu]",count);
} while (count < N);
printf ("\r[%lu]",count);
printf ("\ntests PASSED\n");
return EXIT_SUCCESS;
}
它似乎工作正常（至少对于某些随机情况）：它报告：
[10000000000]
tests PASSED
现在的问题是：由于原始约翰哈里森的sqrtf()算法仅使用单精度计算（即类型float），因此在仅使用（转换除外）double时可以减少运算次数精度计算（即类型 double）并且仍然符合 IEEE 754？
附言由于用户 @njuffa 和 @chux - Reinstate Monica 在 FP 方面很强，我邀请他们参与。但是，欢迎所有有能力的 FP 用户。
解决方法
通过双精度代码计算单精度平方根将是低效的，尤其是在硬件不提供本机双精度运算的情况下。
以下假设硬件符合 IEEE-754 (2008)，除了不支持次正规并刷新为零。支持融合乘加 (FMA)。它还假定 ISO-C99 编译器将 float 映射到 IEEE-754 binary32，并将硬件的单精度 FMA 指令映射到标准数学函数 fmaf()。
从最大相对误差为 2-6.75 的倒数平方根的硬件起始近似值，可以通过两次 Newton-Raphson 迭代得到精确到 1 个单精度 ulp 的倒数平方根.将其与原始参数相乘可提供平方根的准确估计。从原始参数中减去这个近似值的平方来计算平方根的近似误差。然后使用该误差对平方根近似值进行修正，从而得到正确舍入的平方根。
然而，这种简单的算法对于由于中间计算中的下溢或溢出而非常小的参数会失效，特别是当底层算法在将次正规数刷新为零的闪存到零模式下运行时。对于这样的参数，我们可以构建一个慢路径代码，将输入扩展到统一，并在计算平方根后相应地缩减结果。此慢路径代码中还添加了用于处理特殊操作数（例如零、无穷大、NaN 和非零的负参数）的代码。
应对无效操作的慢路径代码生成的 NaN 进行调整，以匹配系统现有的操作。例如，对于基于 x86 的系统，这将是一个称为 INDEFINITE 的特殊 QNaN，其位模式为 0xffc00000，而对于运行 CUDA 的 GPU，它将是规范的单精度 NaN，其位模式为 {{ 1}}。
出于性能原因，内联快速路径代码同时使慢速路径代码成为调用的概述子例程可能很有用。应始终针对“黄金”参考实现对具有单个参数的单精度数学函数进行详尽的测试，这在现代硬件上只需几分钟。
0x7fffffff版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。