考虑到硬件限制和使用限制的符合 IEEE 754 标准的 sqrtf() 实现

Question

的后续问题。

上下文：需要实现 IEEE 754 一致性 sqrtf()，同时考虑以下硬件限制和使用限制：

提供特殊指令qseed.f求平方根倒数的近似值（结果精度不低于6.75位，因此始终在±1%以内的准确结果）。
单精度浮点数：

一个。硬件支持（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, z0, half );
    k = fma( x, rd, -s0 );
    h0 = half * rd;
    e = fma( three / two, d, one );
    t0 = fma( d, s0, k );
    s1 = fma( e, t0, s0 );
    c = fma( d, e, one );
    d1 = fma( -s1, s1, x );
    h1 = c * h0;
    s = fma( d1, h1, 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, 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

现在的问题是：由于最初的 John Harrison 的 sqrtf() 算法仅使用单精度计算（即类型 float），因此可以减少仅使用时的运算次数（转换除外） ) 双精度计算（即类型 double）并且仍然符合 IEEE 754？

P.S。由于用户@njuffa和@chux - Reinstate Monica的FP很强，所以邀请他们参与。不过，有能力的FP用户还是欢迎的。

Answer 1

通过双精度代码计算单精度平方根的效率会很低，尤其是在硬件不提供本机双精度运算的情况下。

以下假设硬件符合 IEEE-754 (2008)，除了不支持次正规并刷新为零。支持融合乘加 (FMA)。它进一步假定一个 ISO-C99 编译器将 float 映射到 IEEE-754 binary32，并将硬件的单精度 FMA 指令映射到标准数学函数 fmaf().

从最大相对误差为 2^-6.75 的倒数平方根的硬件开始近似，可以得到精确到 1 单精度 ulp 的倒数平方根两次 Newton-Raphson 迭代。将其与原始参数相乘可提供平方根的准确估计。从原始参数中减去该近似值的平方，以计算平方根的近似误差。然后使用此错误对平方根近似值进行校正，从而得到正确舍入的平方根。

但是，这种简单的算法会因中间计算中的下溢或溢出而导致非常小的参数失效，特别是当底层算法以闪存归零模式运行时，会将次正规数刷新为零。对于这样的论点，我们可以构建一个慢速路径代码，将输入缩放到统一，并在计算出平方根后相应地缩减结果。用于处理特殊操作数（例如零、无穷大、NaN 和非零的负参数）的代码也添加到此慢路径代码中。

慢速路径代码为无效操作生成的 NaN 应进行调整以匹配系统的现有操作。例如，对于基于 x86 的系统，这将是一个名为 INDEFINITE 的特殊 QNaN，其位模式为 0xffc00000，而对于 GPU 运行 CUDA，它将是规范的单精度 NaN 0x7fffffff.

的模式

出于性能原因，在使慢速路径代码成为被调用的概述子例程的同时内联快速路径代码可能很有用。应始终针对“黄金”参考实现对具有单个参数的单精度数学函数进行详尽测试，这在现代硬件上只需几分钟。

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <math.h>

float uint32_as_float (uint32_t);
uint32_t float_as_uint32 (float);
float qseedf (float);
float sqrtf_slowpath (float);

/* Square root computation for IEEE-754 binary32 mapped to 'float' */
float my_sqrtf (float arg)
{
    const uint32_t upper = float_as_uint32 (0x1.fffffep+127f);
    const uint32_t lower = float_as_uint32 (0x1.000000p-102f);
    float rsq, sqt, err;

    /* use fastpath computation if argument in [0x1.0p-102, 0x1.0p+128) */
    if ((float_as_uint32 (arg) - lower) <= (upper - lower)) {
        /* generate low-accuracy approximation to rsqrt(arg) */
        rsq = qseedf (arg);
        
        /* apply two Newton-Raphson iterations with quadratic convergence */
        rsq = fmaf (fmaf (-0.5f * arg * rsq, rsq, 0.5f), rsq, rsq);
        rsq = fmaf (fmaf (-0.5f * arg * rsq, rsq, 0.5f), rsq, rsq);
        
        /* compute sqrt from rsqrt, round result to nearest or even */
        sqt = rsq * arg;
        err = fmaf (sqt, -sqt, arg);
        sqt = fmaf (0.5f * rsq, err, sqt);
    } else {
        sqt = sqrtf_slowpath (arg);
    }
    return sqt;
}

/* interprete bit pattern of 32-bit unsigned integer as IEEE-754 binary32 */
float uint32_as_float (uint32_t a)
{
    float r;
    memcpy (&r, &a, sizeof r);
    return r;
}

/* interprete bit pattern of  IEEE-754 binary32 as a 32-bit unsigned integer */
uint32_t float_as_uint32 (float a)
{
    uint32_t r;
    memcpy (&r, &a, sizeof r);
    return r;
}

/* simulate low-accuracy hardware approximation to 1/sqrt(a) */
float qseedf (float a)
{
    float r = 1.0f / sqrtf (a);
    r = uint32_as_float (float_as_uint32 (r) & ~0x1ffff);
    return r;
}

/* square root computation suitable for all IEEE-754 binary32 arguments */
float sqrtf_slowpath (float arg)
{
    const float FP32_INFINITY = uint32_as_float (0x7f800000);
    const float FP32_QNAN = uint32_as_float (0xffc00000); /* system specific */
    const float scale_in  = 0x1.0p+26f;
    const float scale_out = 0x1.0p-13f;
    float rsq, err, sqt;

    if (arg < 0.0f) {
        return FP32_QNAN;
    } else if ((arg == 0.0f) || !(fabsf (arg) < FP32_INFINITY)) { /* Inf, NaN */
        return arg + arg;
    } else {
        /* scale subnormal arguments towards unity */
        arg = arg * scale_in;
        
        /* generate low-accuracy approximation to rsqrt(arg) */
        rsq = qseedf (arg);
        
        /* apply two Newton-Raphson iterations with quadratic convergence */
        rsq = fmaf (fmaf (-0.5f * arg * rsq, rsq, 0.5f), rsq, rsq);
        rsq = fmaf (fmaf (-0.5f * arg * rsq, rsq, 0.5f), rsq, rsq);
        
        /* compute sqrt from rsqrt, round to nearest or even */
        sqt = rsq * arg;
        err = fmaf (sqt, -sqt, arg);
        sqt = fmaf (0.5f * rsq, err, sqt);

        /* compensate scaling of argument by counter-scaling the result */
        sqt = sqt * scale_out;
        
        return sqt;
    }
}

int main (void)
{
    uint32_t ai, resi, refi;
    float a, res, reff;
    double ref;

    ai = 0x00000000;
    do {
        a = uint32_as_float (ai);
        res = my_sqrtf (a);
        ref = sqrt ((double)a);
        reff = (float)ref;
        resi = float_as_uint32 (res);
        refi = float_as_uint32 (reff);
        if (resi != refi) {
            printf ("error @ %08x %15.8e   res=%08x %15.8e  ref=%08x %15.8e\n",
                    ai, a, resi, res, refi, reff);
            return EXIT_FAILURE;
        }
        

        ai++;
    } while (ai);
    return EXIT_SUCCESS;
}

考虑到硬件限制和使用限制的符合 IEEE 754 标准的 sqrtf() 实现

IEEE 754 conformant sqrtf() implementation taking into account hardware restrictions and usage limitations

c

floating-point

precision

ieee-754

sqrt