我找到了一个相当奇怪但有效的浮点数平方根近似值； 我真的不明白。 有人能解释一下为什么这段代码有效吗？
float sqrt(float f)
{
const int result = 0x1fbb4000 + (*(int*)&f >> 1);
return *(float*)&result;
}
我已经测试了一下，它的值与std::sqrt()输出的值相差约1到3％。我知道Quake III的快速反平方根算法，我猜这里也有类似的东西（没有牛顿迭代），但我真的很希望得到一个关于如何工作的解释。
(*(int*)&f >> 1) 右移 f 的位表示。这几乎将指数除以二，大约相当于取平方根。1
为什么是“几乎”？在IEEE-754中，实际的指数是 e - 127。2 要将其除以二，我们需要 e/2 - 64，但上述近似只给出了 e/2 - 127。因此，我们需要加上63到结果指数中。这是由那个神奇常数（0x1fbb4000）的30-23位贡献的。
我想剩下的神奇常数位被选择为最小化幅度范围内的最大误差或类似的东西。然而，不清楚它是通过分析、迭代还是启发式确定的。
float。float表示的字节顺序。sqrtf提供了有用的加速！），否则应避免使用此方法。
请看Oliver Charlesworth的解释，了解为什么这个方法几乎可行。我在回应评论中提出的问题。
由于有几个人指出这个方法不可移植，下面是一些方法可以使其更具可移植性，或者至少让编译器告诉你它是否可行。
首先，C++允许您在编译时检查std::numeric_limits<float>::is_iec559，例如在static_assert中。您还可以检查sizeof(int) == sizeof(float)，如果int是64位，则不成立，但您真正想要做的是使用uint32_t，如果存在，则始终恰好为32位宽，具有明确定义的移位和溢出行为，并且如果您的奇怪架构没有这样的整数类型，将导致编译错误。无论哪种方式，您还应该static_assert()类型具有相同的大小。静态断言没有运行时成本，如果可能的话，您应该始终以这种方式检查前提条件。
#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并消除运行时测试，但它只能将此版本优化为单指令存根。
#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;
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的推理是：
从声明类型不同的对象中访问
reinterpret_cast的结果仍然是未定义行为，但至少我们可以看到一些棘手的东西正在发生。
根据评论：C++20 引入了 std::bit_cast，它将对象表示转换为具有未指定而不是未定义的行为的不同类型。这并不保证您的实现将使用与此代码期望的 float 和 int 相同的格式，但它不会给编译器任意打破您的程序的机会，因为在其中一行中存在技术上未定义的行为。它还可以为您提供一个 constexpr 转换。
memcpy() 的版本，并解释了您提出的问题。 - Davislor令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，并确定它是否最优。
float时，尾数的MSbit未编码，仅在正常的float中假定为1。这会影响OP的float sqrt(float f)，但在INT(x)中没有考虑到。 - chux - Reinstate Monica添加一个维基测试工具来测试所有float。
对于许多float，该近似值在4%以内，但对于次正常数则非常差。结果可能因人而异
Worst:1.401298e-45 211749.20%
Average:0.63%
Worst:1.262738e-38 3.52%
Average:0.02%
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;
}
sqrtf()相比如何？它可能是一个快速的近似值吗？如果我们需要“还算可以”的近似值，但平均偏差为0.02％的平方根是可以接受的，那么在实时物理模拟中就会有价值。 - Deliothsqrt_apx() 显然更快。对于先进的处理器，可能会更快，由于优化代码和并行处理。需要设置特定的情况。请记住，sqrt_apx(0.0) 不是 0，这可能会造成实际问题。这一切都非常依赖于具体情况。也许您可以尝试模拟并发布您的结果？ - chux - Reinstate Monicaf 向上移动两个二进制数时，*(int*)&f >> 1 向上移动一个二进制数）。 - Pascal Cuoqfloat，但测试所有内容并不太耗时。使用double，你的想法更有优势，因为它具有2^64个组合。虽然B_1和B_2很重要，但是还需要评估与0x1fbb4000相加后导致幂指数变化的选择FP值。 - chux - Reinstate Monica
0x5f3759df的一种朋友。 - Eugene Sh.f的位表示向右移动一位大致相当于将指数除以二，这等同于取平方根。其他的部分可能是通过魔法来提高尾数范围内的精度。 - Oliver Charlesworth