如何浮点平方根逼近工作？
我发现了一个相当奇怪的但是工作的float的平方根逼近; 我真的不明白。 有人可以解释我为什么这个代码的作品？
float sqrt(float f) { const int result = 0x1fbb4000 + (*(int*)&f >> 1); return *(float*)&result; }
我已经testing了一下， 它输出的值从std::sqrt()大约1到3％ 。 我知道Quake III的快速倒数平方根 ，我猜这是类似的东西（没有牛顿迭代），但是我真的很感激它的解释。
（nota：我已经标记了c和c ++，因为它们都是有效的，ish（参见注释）C和C ++代码）
(*(int*)&f >> 1)右移f的按位表示。 这几乎将指数除以2，这大致相当于取平方根。 1
为什么几乎 ？ 在IEEE-754中，实际指数是e-127 。 2除以二，我们需要e / 2 – 64 ，但上面的近似只给我们e / 2 – 127 。 所以我们需要在所得到的指数上加上63。 这是由魔术常数（ 0x1fbb4000 ）的位30-23所贡献的。
我会想象魔术常数的其余部分已经被select来最小化尾数范围内的最大误差，或类似的东西。 然而，目前还不清楚这是否是分析性的，迭代式的或启发式的。
值得指出的是，这种方法有些不可移植。 它（至less）做出以下假设：
- 该平台使用单精度IEEE-754
float。 -
float表示的字节顺序。 - 由于这种方法违反了C / C ++的严格别名规则，你将不会受到未定义行为的影响。
因此，除非您确定它在您的平台上提供了可预测的行为（事实上，它提供了一个有用的加速比与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
请参阅奥利弗·查尔斯沃斯（Oliver Charlesworth）的解释，为什么这几乎成功 我正在处理评论中提出的一个问题。
既然有几个人指出了这个不可移植性，这里有一些方法可以使它更加便携，或者至less让编译器告诉你它是否不起作用。
首先，C ++允许你在编译时检查std::numeric_limits<float>::is_iec559 ，比如static_assert 。 你也可以检查sizeof(int) == sizeof(float) ，如果int是64位，那么这将不是真的，但是你真正想要做的是使用uint32_t ，如果它存在的话总是正好32位宽，如果你的怪异架构没有这样的整数types，将会有明确定义的行为和移位溢出，并且会导致编译错误。 无论哪种方式，你也应该static_assert() ，types具有相同的大小。 静态断言没有运行时间的成本，你应该总是这样检查你的先决条件，如果可能的话。
不幸的是，是否将float的位转换为uint32_t并进行移位是大端，小端或者都不能作为编译时常量expression式进行计算。 在这里，我将运行时检查放在依赖于它的代码部分，但是您可能需要将其放入初始化中并执行一次。 实际上，gcc和clang都可以在编译时优化这个testing。
你不想使用不安全的指针强制转换，并且有一些我在现实世界中工作过的系统可能会因为总线错误而导致程序崩溃。 memcpy()是转换对象表示的最方便的方法。 在我的下面的例子中，我用一个uniontypes双关键字来处理任何实际存在的实现。 （语言律师反对，但没有一个成功的编译器会永远破坏那么多的遗留代码。）如果你必须做一个指针转换（参见下面），那就是alignas() 。 但是，无论如何，结果都是实现定义的，这就是为什么我们检查转换和移动testing值的结果。
不pipe怎样，你可能不会在现代的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>()的另一个定义，避免了types双关。 你也可以在现代C语言中实现types双关语，在标准允许的情况下，调用函数作为extern "C" 。 我认为打字比memcpy()更优雅，types安全，符合这个程序的准function风格。 我也不认为你获得了很多，因为你仍然可以从一个假设的陷阱表示中得到未定义的行为。 此外，铿锵++ 3.9.1-S能够静态分析typesis_little_endian版本，优化variablesis_little_endian为常量0x1 ，并消除运行时testing，但它只能优化这个版本，指令存根。
但更重要的是，这个代码并不能保证在每个编译器上都可移植。 例如，一些旧电脑甚至无法正确地处理32位内存。 但在这些情况下，它应该不能编译并告诉你为什么。 没有任何编译器会突然打破大量遗留代码。 虽然标准在技术上允许这样做，并且仍然认为它符合C ++ 14，但它只会发生在与我们预期完全不同的架构上。 如果我们的假设是如此无效，以至于一些编译器将float和32位无符号整数之间的types转换为危险的bug，我真的怀疑这个代码背后的逻辑将坚持，如果我们只是使用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); }
我testing过的编译器也可以优化这个折叠常量。 Stroustrup的推理是[原文如此]：
将
reinterpret_cast的结果从对象声明types中获取到一个不同的types仍然是未定义的行为，但是至less我们可以看到一些棘手的事情正在发生。
让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;
findσ，使常量等于0x1fbb4000，并确定它是否是最优的。
添加一个wikitesting工具来testing所有的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 % 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
testing代码
#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; }
