YSC*_*YSC 51 c c++ floating-point optimization ieee-754
我找到了一个相当奇怪但工作的平方根逼近floats; 我真的不明白.有人能解释一下为什么这段代码有效吗？
float sqrt(float f)
{
const int result = 0x1fbb4000 + (*(int*)&f >> 1);
return *(float*)&result;
}
我测试了一下它输出的值std::sqrt()约为1到3%.我知道Quake III的快速反平方根,我想这里有类似的东西(没有牛顿迭代),但我真的很感激它的工作原理.
Oli*_*rth 70
(*(int*)&f >> 1)右移位的按位表示f.这几乎将指数除以2,这大约相当于取平方根.1
为何几乎？在IEEE-754中,实际指数是e-127.2 要将此除以2,我们需要e/2 - 64,但上述近似值仅给出e/2 - 127.所以我们需要在结果指数上加上63.这是由该魔术常量(0x1fbb4000)的位30-23贡献的.
我想是已经选择了魔法常数的剩余部分来最小化尾数范围内的最大误差,或类似的东西.然而,尚不清楚它是通过分析,迭代还是启发式确定的.
值得指出的是,这种方法有点不便携.它(至少)做出以下假设:
float.float代表性的字节顺序.因此,应该避免它,除非你确定它在你的平台上提供了可预测的行为(事实上,它提供了一个有用的加速比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
Dav*_*lor 13
见为什么奥利弗·查尔斯沃思解释几乎工作.我正在解决评论中提出的问题.
由于有几个人已经指出了这种不可移植性,这里有一些方法可以使它更具可移植性,或者至少让编译器告诉你它是否不起作用.
首先,C++允许您std::numeric_limits<float>::is_iec559在编译时检查,例如在static_assert.您也可以检查一下sizeof(int) == sizeof(float),如果int是64位则不会是真的,但您真正想要做的是使用uint32_t,如果它存在,则总是正好是32位宽,将具有明确定义的行为和溢出,如果您的奇怪架构没有这样的整数类型,将导致编译错误.无论哪种方式,您还应该static_assert()使用相同大小的类型.静态断言没有运行时成本,如果可能,您应该始终以这种方式检查前提条件.
不幸的是,是否将a中的位转换float为a uint32_t和shift的测试是big-endian,little-endian或者都不能计算为编译时常量表达式.在这里,我将运行时检查放在依赖于它的代码部分,但您可能希望将其置于初始化中并执行一次.实际上,gcc和clang都可以在编译时优化此测试.
你不想使用不安全的指针转换,并且我在现实世界中有一些系统可能会因为总线错误而导致程序崩溃.转换对象表示的最大可移植方式是memcpy().在下面的示例中,我使用a来键入pun 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中实现type-pun,它是标准允许的,并将函数称为extern "C".我认为类型惩罚更优雅,类型安全并且与该程序的准功能风格一致memcpy().我也不认为你获得了太多,因为你仍然可以从假设的陷阱表示中得到未定义的行为.此外,clang ++ 3.9.1 -O -S能够静态分析类型 - 双关语版本,将变量优化为is_little_endian常量0x1,并消除运行时测试,但它只能将此版本优化为单指令存根.
但更重要的是,这些代码不能保证在每个编译器上都可以移植.例如,一些旧计算机甚至无法准确地处理32位内存.但在这些情况下,它应该无法编译并告诉你原因.没有任何编译器会突然间无缘无故地破坏大量的遗留代码.虽然标准在技术上允许这样做,并且仍然说它符合C++ 14,但它只会发生在与我们期望的完全不同的架构上.如果我们的假设是如此无效,以至于某些编译器会将a 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);
}
我测试的编译器也可以将其优化为折叠常数.Stroustrup的推理是[原文如此]:
reinterpret_cast从声明类型的对象访问一个不同类型的结果仍然是未定义的行为,但至少我们可以看到一些棘手的事情正在发生.
设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 % 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;
}