这个浮点平方根近似值是如何工作的？
How does this float square root approximation work?
我发现了一个相当奇怪但有效的 floats 平方根近似值；我真的不明白。谁能解释一下为什么这段代码有效？
float sqrt(float f)
{
const int result = 0x1fbb4000 + (*(int*)&f >> 1);
return *(float*)&result;
}
我测试了一下 it outputs values off of std::sqrt() by about 1 to 3%. I know of the Quake III's fast inverse square root，我猜它与这里类似（没有牛顿迭代），但我真的很感谢解释 它是如何工作的.
(*(int*)&f >> 1) 右移 f 的按位表示。这个几乎将指数除以2，大约相当于取平方根1
为什么几乎？在 IEEE-754 中，实际指数是 e - 127。2 要将其除以二，我们需要 e/2 - 64，但是上面的近似只给了我们 e/2 - 127。所以我们需要将 63 添加到结果指数中。这是由该魔法常量 (0x1fbb4000) 的第 30-23 位贡献的。
我想已经选择了魔法常数的剩余位来最小化尾数范围内的最大误差，或类似的东西。但是，尚不清楚它是通过分析、迭代还是启发式确定的。
值得指出的是，这种方法有些不可移植。它（至少）做出以下假设：
- 该平台使用单精度 IEEE-754
float。
float 表示的字节顺序。- 您不会受到未定义行为的影响，因为这种方法违反了 C/C++ 的 strict-aliasing rules。
因此应该避免使用它，除非您确定它在您的平台上提供可预测的行为（事实上，与 sqrtf 相比，它提供了有用的加速）。
1. sqrt(a^b) = (a^b)^0.5 = a^(b/2)
2。参见例如https://en.wikipedia.org/wiki/Single-precision_floating-point_format#Exponent_encoding
令 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% 以内，但对于次正规数则非常差。 YMMV
Worst:1.401298e-45 211749.20%
Average:0.63%
Worst:1.262738e-38 3.52%
Average:0.02%
注意参数为 +/-0.0 时，结果不为零。
printf("% e % e\n", sqrtf(+0.0), sqrt_apx(0.0)); // 0.000000e+00 7.930346e-20
printf("% e % e\n", 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;
}
请参阅 Oliver Charlesworth 对 almost 为何有效的解释。我正在解决评论中提出的问题。
由于几个人已经指出了它的不可移植性，这里有一些方法可以使它更具可移植性，或者至少让编译器告诉你它是否不起作用。
首先，C++ 允许您在编译时检查 std::numeric_limits<float>::is_iec559，例如在 static_assert 中。您还可以检查 sizeof(int) == sizeof(float)，如果 int 是 64 位，这将不正确，但您真正想要做的是使用 uint32_t，如果它存在，它将始终是32 位宽，将具有明确定义的移位和溢出行为，如果您的怪异体系结构没有这种整数类型，则会导致编译错误。无论哪种方式，您还应该 static_assert() 类型具有相同的大小。静态断言没有 运行 时间成本，如果可能，您应该始终以这种方式检查您的先决条件。
不幸的是，将 float 中的位转换为 uint32_t 和移位的测试是大端、小端还是两者都不是编译时常量表达式.这里，我把运行-time check放在依赖它的代码部分，但你可能想把它放在初始化中，做一次。实际上，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++ Core Guidelines 中推荐了 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);
}
我测试的编译器也可以将其优化为折叠常量。 Stroustrup 的推理是 [原文]:
Accessing the result of an reinterpret_cast to a different type from the objects declared type is still undefined behavior, but at least we can see that something tricky is going on.
更新
来自评论：C++20 引入了 std::bit_cast，它将对象表示转换为具有 未指定 的不同类型，而不是 未定义, 行为。这并不能保证您的实现将使用与此代码预期相同的 float 和 int 格式，但它不会让编译器全权委托任意中断您的程序，因为存在技术上未定义的行为在它的一行中。它还可以为您提供 constexpr 转换。
我发现了一个相当奇怪但有效的 floats 平方根近似值；我真的不明白。谁能解释一下为什么这段代码有效？
float sqrt(float f)
{
const int result = 0x1fbb4000 + (*(int*)&f >> 1);
return *(float*)&result;
}
我测试了一下 it outputs values off of std::sqrt() by about 1 to 3%. I know of the Quake III's fast inverse square root，我猜它与这里类似（没有牛顿迭代），但我真的很感谢解释 它是如何工作的.
(*(int*)&f >> 1) 右移 f 的按位表示。这个几乎将指数除以2，大约相当于取平方根1
为什么几乎？在 IEEE-754 中，实际指数是 e - 127。2 要将其除以二，我们需要 e/2 - 64，但是上面的近似只给了我们 e/2 - 127。所以我们需要将 63 添加到结果指数中。这是由该魔法常量 (0x1fbb4000) 的第 30-23 位贡献的。
我想已经选择了魔法常数的剩余位来最小化尾数范围内的最大误差，或类似的东西。但是，尚不清楚它是通过分析、迭代还是启发式确定的。
值得指出的是，这种方法有些不可移植。它（至少）做出以下假设：
- 该平台使用单精度 IEEE-754
float。 float表示的字节顺序。- 您不会受到未定义行为的影响，因为这种方法违反了 C/C++ 的 strict-aliasing rules。
因此应该避免使用它，除非您确定它在您的平台上提供可预测的行为（事实上，与 sqrtf 相比，它提供了有用的加速）。
1. sqrt(a^b) = (a^b)^0.5 = a^(b/2)
2。参见例如https://en.wikipedia.org/wiki/Single-precision_floating-point_format#Exponent_encoding
令 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% 以内，但对于次正规数则非常差。 YMMV
Worst:1.401298e-45 211749.20%
Average:0.63%
Worst:1.262738e-38 3.52%
Average:0.02%
注意参数为 +/-0.0 时，结果不为零。
printf("% e % e\n", sqrtf(+0.0), sqrt_apx(0.0)); // 0.000000e+00 7.930346e-20
printf("% e % e\n", 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;
}
请参阅 Oliver Charlesworth 对 almost 为何有效的解释。我正在解决评论中提出的问题。
由于几个人已经指出了它的不可移植性，这里有一些方法可以使它更具可移植性，或者至少让编译器告诉你它是否不起作用。
首先，C++ 允许您在编译时检查 std::numeric_limits<float>::is_iec559，例如在 static_assert 中。您还可以检查 sizeof(int) == sizeof(float)，如果 int 是 64 位，这将不正确，但您真正想要做的是使用 uint32_t，如果它存在，它将始终是32 位宽，将具有明确定义的移位和溢出行为，如果您的怪异体系结构没有这种整数类型，则会导致编译错误。无论哪种方式，您还应该 static_assert() 类型具有相同的大小。静态断言没有 运行 时间成本，如果可能，您应该始终以这种方式检查您的先决条件。
不幸的是，将 float 中的位转换为 uint32_t 和移位的测试是大端、小端还是两者都不是编译时常量表达式.这里，我把运行-time check放在依赖它的代码部分，但你可能想把它放在初始化中，做一次。实际上，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++ Core Guidelines 中推荐了 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);
}
我测试的编译器也可以将其优化为折叠常量。 Stroustrup 的推理是 [原文]:
Accessing the result of an
reinterpret_castto a different type from the objects declared type is still undefined behavior, but at least we can see that something tricky is going on.
更新
来自评论：C++20 引入了 std::bit_cast，它将对象表示转换为具有 未指定 的不同类型，而不是 未定义, 行为。这并不能保证您的实现将使用与此代码预期相同的 float 和 int 格式，但它不会让编译器全权委托任意中断您的程序，因为存在技术上未定义的行为在它的一行中。它还可以为您提供 constexpr 转换。