这种浮点平方根近似是如何工作的

(*(int*)&f >> 1)右移f的按位表示形式。 这几乎将指数除以 2，这大约相当于取平方根。¹

为什么几乎？在 IEEE-754 中，实际指数是e - 127。^{阿拉伯数字}要将其除以 2，我们需要e/2 - 64，但上面的近似值只给了我们e/2 - 127。 因此，我们需要将 63 添加到生成的指数中。 这是由该魔术常数(0x1fbb4000)的第30-23位贡献的。

我想魔术常数的剩余部分已被选择以最小化尾数范围内的最大误差，或类似的东西。 但是，目前尚不清楚它是通过分析，迭代还是启发式确定的。

值得指出的是，这种方法有点不可移植。它(至少)做出以下假设：

• 该平台使用单精度IEEE-754进行float。
• float表示的字节序。
• 由于此方法违反了 C/C++ 的严格别名规则，因此您将不受未定义行为的影响。

因此，应该避免使用它，除非您确定它在您的平台上提供了可预测的行为(事实上，它提供了有用的加速与sqrtf！

_{^b)^0.5 = a^(b/2)}

_{2. 例如，请参阅 https://en.wikipedia.org/wiki/Single-precision_floating-point_format#Exponent_encoding}

参见奥利弗·查尔斯沃思(Oliver Charlesworth)对为什么这几乎有效的解释。 我正在解决评论中提出的问题。

由于有几个人指出了它的不可移植性，这里有一些方法可以使它更具可移植性，或者至少让编译器告诉你它是否不起作用。

首先，C++允许您在编译时检查std::numeric_limits<float>::is_iec559，例如在static_assert中。 您还可以检查sizeof(int) == sizeof(float)，如果int是 64 位，则不正确，但您真正想做的是使用uint32_t，如果它存在，它将始终正好是 32 位宽，将具有明确定义的移位和溢出行为，如果您奇怪的架构没有这样的积分类型，则会导致编译错误。 无论哪种方式，您还应该static_assert()类型具有相同的大小。 静态断言没有运行时成本，如果可能，应始终以这种方式检查前提条件。

遗憾的是，将float中的位转换为uint32_t和移位是大端序、小端序还是两者都不是的测试不能计算为编译时常量表达式。 在这里，我将运行时检查放在依赖于它的代码部分，但您可能希望将其放入初始化中并执行一次。 在实践中，gcc 和 clang 都可以在编译时优化这个测试。

您不想使用不安全的指针强制转换，并且我在现实世界中工作过的一些系统可能会因总线错误而使程序崩溃。 转换对象表示的最大可移植方法是使用memcpy(). 在下面的示例中，我用union键入双关语，它适用于任何实际存在的实现。 (语言律师反对它，但没有一个成功的编译器会默默地破坏那么多遗留代码。 如果必须进行指针转换(见下文)，则alignas(). 但无论你怎么做，结果都将是实现定义的，这就是我们检查转换和移动测试值的结果的原因。

无论如何，并不是说您可能会在现代CPU上使用它，这里有一个滔滔不绝的C++14版本，可以检查那些不可移植的假设：

#include <cassert>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <limits>
#include <vector>
using std::cout;
using std::endl;
using std::size_t;
using std::sqrt;
using std::uint32_t;
template <typename T, typename U>
inline T reinterpret(const U x)
/* Reinterprets the bits of x as a T.  Cannot be constexpr
* in C++14 because it reads an inactive union member.
*/
{
static_assert( sizeof(T)==sizeof(U), "" );
union tu_pun {
U u = U();
T t;
};

const tu_pun pun{x};
return pun.t;
}
constexpr float source = -0.1F;
constexpr uint32_t target = 0x5ee66666UL;
const uint32_t after_rshift = reinterpret<uint32_t,float>(source) >> 1U;
const bool is_little_endian = after_rshift == target;
float est_sqrt(const float x)
/* A fast approximation of sqrt(x) that works less well for subnormal numbers.
*/
{
static_assert( std::numeric_limits<float>::is_iec559, "" );
assert(is_little_endian); // Could provide alternative big-endian code.

/* The algorithm relies on the bit representation of normal IEEE floats, so
* a subnormal number as input might be considered a domain error as well?
*/
if ( std::isless(x, 0.0F) || !std::isfinite(x) )
return std::numeric_limits<float>::signaling_NaN();

constexpr uint32_t magic_number = 0x1fbb4000UL;
const uint32_t raw_bits = reinterpret<uint32_t,float>(x);
const uint32_t rejiggered_bits = (raw_bits >> 1U) + magic_number;
return reinterpret<float,uint32_t>(rejiggered_bits);
}
int main(void)
{  
static const std::vector<float> test_values{
4.0F, 0.01F, 0.0F, 5e20F, 5e-20F, 1.262738e-38F };

for ( const float& x : test_values ) {
const double gold_standard = sqrt((double)x);
const double estimate = est_sqrt(x);
const double error = estimate - gold_standard;

cout << "The error for (" << estimate << " - " << gold_standard << ") is "
<< error;
if ( gold_standard != 0.0 && std::isfinite(gold_standard) ) {
const double error_pct = error/gold_standard * 100.0;
cout << " (" << error_pct << "%).";
} else
cout << '.';
cout << endl;
}
return EXIT_SUCCESS;
}

更新

这是避免类型双关语的reinterpret<T,U>()的另一种定义。 您还可以在标准允许的现代 C 中实现类型双关语，并将函数调用为extern "C"。 我认为类型双关语比memcpy()更优雅、类型安全且符合该程序的准功能风格。 我也不认为你获得太多好处，因为你仍然可能从假设的陷阱表示中获得未定义的行为。此外，clang++ 3.9.1 -O -S 能够静态分析类型双关版本，优化变量is_little_endian常量0x1，并消除运行时测试，但它只能将此版本优化为单指令存根。

但更重要的是，此代码不能保证在每个编译器上都可以移植。 例如，一些旧计算机甚至无法准确寻址 32 位内存。 但是在这些情况下，它应该无法编译并告诉您原因。 没有编译器会突然无缘无故地破坏大量遗留代码。 尽管该标准在技术上允许这样做，并且仍然说它符合C++14，但它只会发生在与我们预期的非常不同的架构上。 如果我们的假设是如此无效，以至于某些编译器会将float和 32 位无符号整数之间的类型双关语变成一个危险的错误，我真的怀疑如果我们只使用memcpy()，这段代码背后的逻辑会站得住脚。 我们希望该代码在编译时失败，并告诉我们原因。

#include <cassert>
#include <cstdint>
#include <cstring>
using std::memcpy;
using std::uint32_t;
template <typename T, typename U> inline T reinterpret(const U &x)
/* Reinterprets the bits of x as a T.  Cannot be constexpr
* in C++14 because it modifies a variable.
*/
{
static_assert( sizeof(T)==sizeof(U), "" );
T temp;

memcpy( &temp, &x, sizeof(T) );
return temp;
}
constexpr float source = -0.1F;
constexpr uint32_t target = 0x5ee66666UL;
const uint32_t after_rshift = reinterpret<uint32_t,float>(source) >> 1U;
extern const bool is_little_endian = after_rshift == target;

然而，Stroustrup等人在C++核心指南中建议采用reinterpret_cast：

#include <cassert>
template <typename T, typename U> inline T reinterpret(const U x)
/* Reinterprets the bits of x as a T.  Cannot be constexpr
* in C++14 because it uses reinterpret_cast.
*/
{
static_assert( sizeof(T)==sizeof(U), "" );
const U temp alignas(T) alignas(U) = x;
return *reinterpret_cast<const T*>(&temp);
}

我测试的编译器也可以将其优化为折叠常量。 斯特劳斯特鲁普的推理是[原文如此]：

访问与声明类型的对象不同的类型的reinterpret_cast结果仍然是未定义的行为，但至少我们可以看到正在发生一些棘手的事情。

更新

从注释中：C++20 引入了std::bit_cast，它将对象表示转换为具有未指定而非未定义行为的不同类型的对象表示。 这并不能保证您的实现将使用此代码所期望的相同格式的float和int，但它不会让编译器全权任意破坏您的程序，因为它的一行中存在技术上未定义的行为。 它还可以为您提供constexpr转换。

let y = sqrt(x)，

根据对数的性质，log(y) = 0.5 * log(x) (1)

将普通float解释为整数得到 INT(x) = Ix = L * (log(x) + B - σ) (2)

其中 L = 2^N，N 是有效位数，B 是指数偏差，σ 是调整近似的自由因子。

结合 (1) 和 (2) 得到：Iy = 0.5 * (Ix + (L * (B - σ)))

在代码中写成(*(int*)&x >> 1) + 0x1fbb4000;

找到σ，使常量等于 0x1fbb4000 并确定它是否最优。

添加一个 wiki 测试工具来测试所有float。

许多float的近似值在4%以内，但次正规数的近似值非常差。 扬子晚报

Worst:1.401298e-45 211749.20%
Average:0.63%
Worst:1.262738e-38 3.52%
Average:0.02%

请注意，当参数为 +/-0.0 时，结果不为零。

printf("% e % en", sqrtf(+0.0), sqrt_apx(0.0));  //  0.000000e+00  7.930346e-20
printf("% e % en", sqrtf(-0.0), sqrt_apx(-0.0)); // -0.000000e+00 -2.698557e+19

测试代码

#include <float.h>
#include <limits.h>
#include <math.h>
#include <stddef.h>
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
float sqrt_apx(float f) {
const int result = 0x1fbb4000 + (*(int*) &f >> 1);
return *(float*) &result;
}
double error_value = 0.0;
double error_worst = 0.0;
double error_sum = 0.0;
unsigned long error_count = 0;
void sqrt_test(float f) {
if (f == 0) return;
volatile float y0 = sqrtf(f);
volatile float y1 = sqrt_apx(f);
double error = (1.0 * y1 - y0) / y0;
error = fabs(error);
if (error > error_worst) {
error_worst = error;
error_value = f;
}
error_sum += error;
error_count++;
}
void sqrt_tests(float f0, float f1) {
error_value = error_worst = error_sum = 0.0;
error_count = 0;
for (;;) {
sqrt_test(f0);
if (f0 == f1) break;
f0 = nextafterf(f0, f1);
}
printf("Worst:%e %.2f%%n", error_value, error_worst*100.0);
printf("Average:%.2f%%n", error_sum / error_count);
fflush(stdout);
}
int main() {
sqrt_tests(FLT_TRUE_MIN, FLT_MIN);
sqrt_tests(FLT_MIN, FLT_MAX);
return 0;
}