IEEE 754平方根与Newton-Raphson

#ifndef DECOMPOSE_H
#define DECOMPOSE_H 1

#ifdef __cplusplus
extern "C" {
#endif

#include <float.h>

#define MANTISSA_SIZE 52
#define EXPONENT_SIZE 11
#define EXPONENT_BIAS 0x3ff
#define TYPE double
#define TYPE_MIN DBL_MIN
#define TYPE_MAX DBL_MAX
#define NAME_SFX(name) name

typedef unsigned long long mantissa;
typedef unsigned short exponent;
typedef unsigned char sign;

#include "decompose_generic.h"

#ifdef __cplusplus
}
#endif

#endif

#ifndef DECOMPOSE_GENERIC_H
#define DECOMPOSE_GENERIC_H 1

#ifdef __cplusplus
extern "C" {
#endif

#define SIGN_SIZE 1

union decompose {
    TYPE f;
    struct decompose_s {
        mantissa m:MANTISSA_SIZE;
        exponent e:EXPONENT_SIZE;
        sign s:SIGN_SIZE;
    } __attribute__((packed)) o;
};

static inline union decompose decompose(TYPE f)
{
    union decompose result;
    result.f = f;
    return result;
}

#define is_nan(x) NAME_SFX(isnan)(x)

#ifdef __cplusplus
}
#endif

#endif

#include <math.h>
#include <errno.h>

#include "decompose.h"

TYPE NAME_SFX(sqrt)(TYPE a)
{
    TYPE result_old2 = 0;
    unsigned short c = 0;
    union decompose r = decompose(a);

    if(a < -(TYPE)0)
    {
        errno = EDOM;
        return NAME_SFX(nan)("");
    }
    if(a == (TYPE)0 || is_nan(a) || a > TYPE_MAX)
    {
        return a;
    }

    /* Divide the exponent by 2 */
    r.o.e -= EXPONENT_BIAS;
    r.o.e = (r.o.e & 1) | (r.o.e >> 1);
    r.o.e += EXPONENT_BIAS;

    /* Newton-Raphson */
    while(result_old2 != r.f)
    {
        if(c % 2 == 0) result_old2 = r.f;
        r.f -= (r.f*r.f-a)/(2*r.f);
        c++;
    }

    return r.f;
}

#1 楼

该答案使用指针广播进行类型优化，以节省空间。实际上，请继续使用联合（在ISO C99中安全，在C ++中作为GNU和MSVC扩展安全）或memcpy（在C和C ++中安全）。仅在MSVC或具有-fno-strict-aliasing的GNU兼容编译器中，这种指针广播才是安全的。初始逼近

打包位字段不仅在这里不是必需的，而且会使结果更糟。您可以执行以下操作：

uint64_t i = *(uint64_t*)&f;           // double to int
i = i - (1023 << 52);                  // remove IEEE754 bias
i = i >> 1;                            // divide exponent by 2
i = i + (1023 << 52);                  // add IEEE754 bias

如果更改操作顺序，则可以在单个常量中处理偏差：

uint64_t i = *(uint64_t*)&f;           // double to int
i = i >> 1;                            // divide exponent by 2
i = i + (1023 << 52) - (1023 << 51);   // remove/add IEEE754 bias

常量可以简化为：

uint64_t i = *(uint64_t*)&f;           // double to int
i = (i >> 1) + (1023 << 51);           // all of the above

请注意，移位时我并没有遮掩，所以指数是最右边的位落入尾数，尾数也会移动。这是一个功能，而不是错误。解释为整数时，故意将IEEE754格式选择为单调的，因此明确允许指数和尾数之间的残留。

按我的近似，映射如下（区间内的线性映射）：

original   number is    exponent   manitssa   orig interval
exponent  in interval   in sqrt   first bit     maps to
   0     [  1.0,  2.0 )    0         0       [  1.0,  1.5 )
   1     [  2.0,  4.0 )    0         1       [  1.5,  2.0 )
   2     [  4.0,  8.0 )    1         0       [  2.0,  3.0 )
   3     [  8.0, 16.0 )    1         1       [  3.0,  4.0 )
   4     [ 16.0, 32.0 )    2         0       [  4.0,  6.0 )

使用您的代码，映射是（在间隔内也是线性的）：

original   number is    exponent  orig interval
exponent  in interval   in sqrt     maps to
   0     [  1.0,  2.0 )    0     [  1.0,  2.0 )
   1     [  2.0,  4.0 )    0     [  1.0,  2.0 )
   2     [  4.0,  8.0 )    1     [  2.0,  4.0 )
   3     [  8.0, 16.0 )    1     [  2.0,  4.0 )
   4     [ 16.0, 32.0 )    2     [  4.0,  8.0 )

Newton- Raphson

给出一个很好的近似值，Newton-Raphson将每次迭代的有效位数翻倍（二次收敛）。上面的近似值提供了大约4位精度（最大误差：6％或〜1/16），因此单精度需要3个Newton-Raphson迭代，而双精度则需要4个迭代。

整个代码可能看起来像这样：

float sqrtf(float x)
{
    float y = x;
    // Approximation
    uint32_t* i = (uint32_t*)&x;
    *i = (*i >> 1) + (127 << 22);
    // Newton-Raphson
    x = (x + y/x) / 2;
    x = (x + y/x) / 2;
    x = (x + y/x) / 2;
    return x;
}

double sqrt(double x)
{
    double y = x;
    // Approximation
    uint64_t* i = (uint64_t*)&x;
    *i = (*i >> 1) + (1023 << 51);
    // Newton-Raphson
    x = (x + y/x) / 2;
    x = (x + y/x) / 2;
    x = (x + y/x) / 2;
    x = (x + y/x) / 2;
    return x;
}

次常态和其他危害

您已经对nan，inf，0和x < 0进行了特殊处理，这很好，但是我们也需要处理非正规数。您和我的代码在遇到这样的数字时都会失败。作为补救措施，我建议您通过2^200缩放次正规数，并通过2^-100缩放结果。我之所以选择2^200，是因为它可以安全地将所有次法线恢复到所有精度（包括四精度）的正常范围内，但其他数字也能很好地工作（但只有两个的幂才能保持全精度）。

还要注意我的代码依赖于IEEE754属性。它不适用于非IEEE754数字，尤其是x86的80位extended double格式。如果要处理此类数字，可以将它们转换为float或double进行近似计算，然后将其转换为extended double进行牛顿-拉夫森计算。

#2 楼

仅供参考，在IEEE754 sqrt中是“基本”操作，需要正确舍入（舍入误差<= 0.5ulp），与+-* /相同。硬件FPU（我认为）总是提供sqrt，如果它们提供其他操作，尤其是通常以类似方式实现（并且具有类似性能）的除法。每个涉及除法的NR迭代都不会更快。比硬件sqrt。我认为这就是为什么x86的近似指令用于倒数和rsqrt而不是sqrt的部分原因：1/sqrt(x)的牛顿迭代公式不涉及除法运算，因此在与mulps相比sqrtps非常慢的旧CPU上，值得使用sqrt(x) ~= x * approx_rsqrt(x)，还可以选择使用对rox_rsqrt结果进行牛顿迭代。也许有时仍然值得在不执行其他任何操作的循环中进行此操作，但通常您应该争取更多的计算强度（在寄存器中热数据时对数据进行更多处理；不要编写仅加载/ sqrt /存储的循环）。 br />
IDK，如果代码可能会收敛到不是正确舍入结果的值，但是您应该通过检查实现的sqrt函数来进行测试。

（以其他方式调用您的函数，例如在GCC中，除非您使用sqrt，否则-fno-builtin-sqrt是一个特殊名称。调用者将在具有x86 sqrtsd之类的目标上内嵌sqrt指令，仅对需要设置errno的输入调用libc函数。或永远不要使用-fno-math-errno进行编译，这总是一个好主意。尽管GNU C确实支持Math-errno，但并非所有libcs​​都支持Math-errno。这是一个非常糟糕的过时设计，没有人可以使用；这就是为什么我们需要fenv.h来检查线程粘性FP异常标志。）


输入错误的初始猜测<1.0

 /* Divide the exponent by 2 */
 r.o.e -= EXPONENT_BIAS;
 r.o.e = (r.o.e & 1) | (r.o.e >> 1);
 r.o.e += EXPONENT_BIAS;

您的指数位域类型是unsigned short 1，所以这是合乎逻辑的右移。

但是无论如何，逻辑右移是您初始近似值中的一个错误。

指数字段是2的补码整数（使用偏差编码，但这只是您已经处理过的编码/解码细节） 。小指数的逻辑右移（例如0.0001）将导致大的正指数从零移入。您最初对0.0001的猜测将接近DBL_MAX，远高于1.0。

您需要对指数进行算术右移，以使其更接近中间指数值（无偏）0或（有偏见）127。因此，使表示值的大小更接近1.0。

大多数C实现都明智地选择了带符号整数类型的>>是算术右移。 ISO C要求实现定义它是什么，但是不幸的是，它们让他们选择逻辑还是算术。 （不同于有符号的左移，它不能是UB）。但是幸运的是，您的代码是用GNU C（不是ISO C）编写的，例如在结构上使用__attribute__((packed))。 GNU C保证有符号整数上的>>是算术的。 （而且带符号的整数类型是2的补码！）因此，唯一可以编译代码的实现是>>上的int可以满足您的要求。


（GCC手册，整数实现）定义的行为）：
带符号的“ >>”通过符号扩展作用于负数。但是您可以放心使用int32_t（保证是2的补码整数类型），也可以仅使用int临时保留用于偏差/无偏计算的字段值。偏差/无偏差可能会越过零（无符号环绕），而不是有符号溢出。


或者更好，请使用@Rainer P的方法：首先不要无偏心，只有无符号偏移，然后再加上一半的偏置常数。 （请参阅他的答案，以了解如何简化到这一点）。在无偏之前进行移位意味着指数是无符号的，所以移位为零是好的。 （我认为即使对于偏零指数接近零（即接近最小可能的指数编码）的微小值，此方法也适用。 net / FloatConverter / IEEE754.html，手动移动指数，然后将1（带进位）加到第二高的指数位。例如2.3509887E-38（有偏指数= 0b0000'0010）变为2.1684043E-19（有偏指数= 0b0100'0001）。

但是当移位整个位模式时，请注意不要从q的符号位中移入1 -0.0。您正在将a == 0.0作为特殊情况进行处理，因此您可以在这里很好，因为-0.0 == 0.0为true。我猜你不能让0.0穿过主NR循环，因为你将被零除。


脚注1：

IDK为什么你请为您的指数位字段类型选择unsigned short而不是unsigned。位字段的基础类型可以比字段宽，而不会损害任何内容，并且某些编译器会在该类型的宽度上进行数学运算。因此，如果unsigned int足够宽，则是一个好主意。 int至少为16位。

使用uint64_t作为位域类型会导致某些编译器（尤其是MSVC）的32位代码运行缓慢，但我认为unsigned int对于指数总是很好的。对于有偏的指数，无符号确实有意义。

有趣的事实：位域类型只有int，而签名则留给实现。您可以使用signed int强制签名。 https://stackoverflow.com/questions/42527387/signed-bit-field-represetation


展开而不是检查c % 2 == 0


编译器仍然可以执行此操作，但是将循环展开2而不是使用计数器进行其他每次迭代的特殊操作可能更有意义。

实际的Newton迭代表达式足够紧凑，以至于与通过c++和if (c%2 == 0)逻辑进行操作相比，对人类重复读取一次可能更简单。

展开时，您可以将while(foo)条件作为if(!foo) break;处理，也可以简单地将其保留，只检查是否保留了每2个牛顿迭代循环一次，除非您需要检查两种可能性以捕获两个值之​​间的交替。

（初始估计很差，无论如何都要进行多次迭代，这可能是值得的。否则，可能不值得，特别是在超标量/乱序CPU上，它可以与没有很多ILP（指令级并行性）；主要是一个长的依赖链的主依赖链并行地检查分支条件。 br />

拉快速路径上的一些检查

您的版本在进行实际工作之前对a进行了4次单独检查。

    if(a < -(TYPE)0)
    {
        errno = EDOM;
        return NAME_SFX(nan)("");
    }
    if(a == (TYPE)0 || is_nan(a) || a > TYPE_MAX)
    {
        return a;
    }

带有一个或两个比较的特殊情况，并在快速路径之外的那个块中将它们分类。 sqrt(NaN)的性能远不如sqrt(1.234)的性能重要。

    if(! a > (TYPE)0)
    {
        if(a == (TYPE)0 || is_nan(a))
            return a;
        // else a < 0
        errno = EDOM;
        return NAME_SFX(nan)("");
    }
    if( a > TYPE_MAX)
    {
        return a;
    }

如果您需要专门处理次态（至少是为了生成初始近似值？），您可以可能使用if (! a>=DBL_MIN)。或者，您可以决定您的软浮点实现根本不处理次态。

IDK此代码的目的是什么；如果它是软浮动库的一部分，则仅使用FP add，mul和div似乎效率低下。您想一起优化多个操作，等等。

但是如果是用于硬件FPU，那么我认为实际上所有FPU都将内置sqrt。

所以我认为这只是一个学习练习。因此，我们将不会针对针对任何特定硬件或编译器进行优化进行详细介绍。

#3 楼

我不明白为什么要自己定义此常量：

#define MANTISSA_SIZE 52

鉴于我们已经假设FLT_RADIX为2，我们可以使用来自<float.h>的适当宏（用于DBL_MANT_DIG的double等）。


我认为这里存在整数溢出的危险：

/* Divide the exponent by 2 */
r.o.e -= EXPONENT_BIAS;
r.o.e = (r.o.e & 1) | (r.o.e >> 1);
r.o.e += EXPONENT_BIAS;

我们最好将其提取到足够大的临时区域中，并对其施加指数偏差：

int e = r.o.e - EXPONENT_BIAS;
e = (e & 1) | (e >> 1);
r.o.e = e + EXPONENT_BIAS;

转移到位，然后使用一半的偏差进行校正；我还没有检查是否等效。

#4 楼

模块化奇数校验

编译器仍然可以执行此操作，但是不能

if(c % 2 == 0)

if (!(c & 1))

？

因素

不是

r.f -= (r.f*r.f-a)/(2*r.f);

等同于：

r.f = (r.f + a/r.f)/2?

优化

根据所使用的编译器以及传递的标志，某些尝试进行的数学运算可能会被优化程序，并且在与硬件紧密耦合的实现中，这可能很重要。您检查过生产的组件了吗？您可以向我们展示您的构建标志吗？