如何实现快速逆平方根没有未定义的行为

How to implement fast inverse sqrt without undefined behavior?

本文关键字：平方根 未定义 何实现 实现      更新时间：2023-10-16

根据我对严格混叠规则的理解，这段快速平方根反比的代码将导致c++中未定义的行为:

float Q_rsqrt( float number )
{
    long i;
    float x2, y;
    const float threehalfs = 1.5F;
    x2 = number * 0.5F;
    y  = number;
    i  = * ( long * ) &y; // type punning
    i  = 0x5f3759df - ( i >> 1 );
    y  = * ( float * ) &i;
    y  = y * ( threehalfs - ( x2 * y * y ) );
    return y;
}

这段代码真的会导致UB吗?如果是，如何以符合标准的方式重新实现?如果不是，为什么不呢?

假设:在调用此函数之前，我们以某种方式检查了浮点数是IEEE 754 32位格式，sizeof(long)==sizeof(float)和平台是little-endian。

标准的兼容方式是std::memcpy。在您指定的假设下，这应该足够符合标准。任何合理的编译器都会把它转换成一堆寄存器移动，如果可能的话。此外，我们还可以使用c++ 11的static_assert和<cstdint>的固定宽度整数类型来减轻(或至少检查)您所做的一些假设。无论如何，端序是无关的，因为我们这里没有处理任何数组，如果整数类型是小端序的，浮点类型也是。

float Q_rsqrt( float number )
{
    static_assert(std::numeric_limits<float>::is_iec559, 
                  "fast inverse square root requires IEEE-comliant 'float'");
    static_assert(sizeof(float)==sizeof(std::uint32_t), 
                  "fast inverse square root requires 'float' to be 32-bit");
    float x2 = number * 0.5F, y = number;
    std::uint32_t i;
    std::memcpy(&i, &y, sizeof(float));
    i  = 0x5f3759df - ( i >> 1 );
    std::memcpy(&y, &i, sizeof(float));
    return y * ( 1.5F - ( x2 * y * y ) );
}

您应该使用memcpy。据我所知，这是唯一符合标准的方式，编译器足够聪明，可以用一个单词move指令替换调用。

我不认为您可以在没有严格标准意义上的UB的情况下做到这一点，仅仅是因为它依赖于float和long的特定值表示。标准没有指定这些(故意的)，并且给予实现自由去做他们认为在这方面合适的事情，特别地将"UB"标签应用于所有依赖于这些事情的行为。

这并不意味着这将是一种"海森堡虫"或"鼻恶魔"类型的UB;最有可能的是，实现将为这种行为提供定义。只要这些定义满足代码的先决条件，就没问题。这里的前提条件是:

  • 实现允许类型双关语。
  • sizeof(float) == sizeof(long) .
  • float的内部表示是这样的，如果这些比特被解释为long, 0x5f3759df和>> 1就具有这种破解功能所需的含义。

然而，以一种完全可移植的方式重写它是不可能的(即没有UB) -例如，对于使用不同浮点运算实现的平台，你会怎么做?只要你依赖于特定平台的细节，你就依赖于标准中未定义但由平台定义的行为。

不。该算法首先假设IEEE754布局，这不是一个有效的标准兼容假设。特定的常数0x5f3759df对于任何其他布局来说都是完全错误的。甚至存在这样一个常数的假设也是有缺陷的。我认为只要你颠倒浮点数内的字段顺序，它就会失效。

我也不确定Q_rsqrt(-0.0f)是否正确工作，即使假设IEEE754。两个零(正负)应该相等，sqrt(x)只有在x<-0.0f

时才未定义
相关文章：