c++ - standard - ¿Cómo funciona esta aproximación de raíz cuadrada flotante?
standard 754 ieee (4)
Encontré una aproximación de raíz cuadrada bastante extraña pero funcional para
float
s;
Realmente no lo entiendo.
¿Alguien puede explicarme por qué funciona este código?
float sqrt(float f)
{
const int result = 0x1fbb4000 + (*(int*)&f >> 1);
return *(float*)&result;
}
Lo probé un poco y
genera valores fuera de
std::sqrt()
en aproximadamente 1 a 3%
.
Sé de la
raíz cuadrada inversa rápida
del Quake III y supongo que es algo similar aquí (sin la iteración de Newton), pero realmente agradecería una explicación de
cómo funciona
.
(nota: lo he etiquetado como c y c ++ ya que es válido-ish (ver comentarios) Código C y C ++)
Agregar un arnés de prueba wiki para probar todo el
float
.
La aproximación está dentro del 4% para muchos
float
, pero muy pobre para los números por debajo de lo normal.
YMMV
Worst:1.401298e-45 211749.20%
Average:0.63%
Worst:1.262738e-38 3.52%
Average:0.02%
Tenga en cuenta que con un argumento de +/- 0.0, el resultado no es cero.
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
Código de prueba
#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;
}
Deje y = sqrt (x),
De las propiedades de los logaritmos se desprende que log (y) = 0.5 * log (x) (1)
Interpretar un
float
normal como un entero da INT (x) = Ix = L * (log (x) + B - σ) (2)
donde L = 2 ^ N, N el número de bits del significado, B es el sesgo del exponente, y σ es un factor libre para ajustar la aproximación.
La combinación de (1) y (2) da: Iy = 0.5 * (Ix + (L * (B - σ)))
Que está escrito en el código como
(*(int*)&x >> 1) + 0x1fbb4000;
Encuentre el σ para que la constante sea igual a 0x1fbb4000 y determine si es óptima.
Vea la explicación de Oliver Charlesworth de por qué esto casi funciona. Estoy abordando un problema planteado en los comentarios.
Dado que varias personas han señalado la no portabilidad de esto, aquí hay algunas formas en que puede hacerlo más portátil, o al menos hacer que el compilador le diga si no funcionará.
Primero, C ++ le permite verificar
std::numeric_limits<float>::is_iec559
en tiempo de compilación, como en un
static_assert
.
También puede verificar que
sizeof(int) == sizeof(float)
, que no será cierto si
int
es de 64 bits, pero lo que realmente quiere hacer es usar
uint32_t
, que si existe siempre tendrá exactamente 32 bits de ancho , tendrá un comportamiento bien definido con cambios y desbordamiento, y provocará un error de compilación si su arquitectura extraña no tiene ese tipo integral.
De cualquier manera, también debe
static_assert()
que los tipos tengan el mismo tamaño.
Las aserciones estáticas no tienen costo de tiempo de ejecución y siempre debe verificar sus condiciones previas de esta manera si es posible.
Desafortunadamente, la prueba de si la conversión de los bits en un
float
a
uint32_t
y el desplazamiento es big-endian, little-endian o ninguno de los dos puede calcularse como una expresión constante en tiempo de compilación.
Aquí, puse la verificación de tiempo de ejecución en la parte del código que depende de él, pero es posible que desee ponerlo en la inicialización y hacerlo una vez.
En la práctica, tanto gcc como clang pueden optimizar esta prueba en tiempo de compilación.
No desea utilizar el puntero inseguro, y hay algunos sistemas en los que he trabajado en el mundo real en los que podría bloquear el programa con un error de bus.
La forma más portátil de convertir representaciones de objetos es con
memcpy()
.
En mi ejemplo a continuación, escribo juegos de palabras con una
union
, que funciona en cualquier implementación existente.
(Los abogados de idiomas se oponen, pero ningún compilador exitoso romperá ese código heredado en
silencio
.) Si debe hacer una conversión de puntero (ver más abajo) hay
alignas()
.
Pero, independientemente de cómo lo haga, el resultado estará definido por la implementación, por lo que verificamos el resultado de convertir y cambiar un valor de prueba.
De todos modos, no es probable que lo use en una CPU moderna, aquí hay una versión de C ++ 14 que verifica esos supuestos no portátiles:
#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;
}
Actualizar
Aquí hay una definición alternativa de
reinterpret<T,U>()
que evita la escritura de tipos.
También puede implementar el juego de palabras en C moderno, donde está permitido por estándar, y llamar a la función como
extern "C"
.
Creo que la
memcpy()
es más elegante, segura y consistente con el estilo casi funcional de este programa que
memcpy()
.
Tampoco creo que ganes mucho, porque aún podrías tener un comportamiento indefinido de una representación hipotética de trampa.
Además, clang ++ 3.9.1 -O -S es capaz de analizar estáticamente la versión de
is_little_endian
, optimizar la variable
is_little_endian
a la constante
0x1
y eliminar la prueba de tiempo de ejecución, pero solo puede optimizar esta versión a una sola- trozo de instrucciones.
Pero lo que es más importante, no se garantiza que este código funcione de manera portátil en cada compilador.
Por ejemplo, algunas computadoras antiguas ni siquiera pueden direccionar exactamente 32 bits de memoria.
Pero en esos casos, no debería compilarse y decirte por qué.
Ningún compilador de repente va a romper una gran cantidad de código heredado sin ninguna razón.
Aunque el estándar técnicamente da permiso para hacerlo y aún dice que se ajusta a C ++ 14, solo sucederá en una arquitectura muy diferente de lo que esperamos.
Y si nuestras suposiciones son tan inválidas que algún compilador va a convertir un juego de palabras entre un
float
y un entero sin signo de 32 bits en un error peligroso, dudo mucho que la lógica detrás de este código se mantenga si usamos
memcpy()
lugar.
Queremos que ese código falle en el momento de la compilación y que nos diga por qué.
#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;
Sin embargo, Stroustrup et al., En las
Pautas principales de C ++
, recomiendan un
reinterpret_cast
lugar:
#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);
}
Los compiladores que probé también pueden optimizar esto a una constante doblada. El razonamiento de Stroustrup es [sic]:
Acceder al resultado de un
reinterpret_casta un tipo diferente de los objetos declarados tipo sigue siendo un comportamiento indefinido, pero al menos podemos ver que algo complicado está sucediendo.
(*(int*)&f >> 1)
desplaza a la derecha la representación bit a bit de
f
.
Esto
casi
divide el exponente por dos, lo que es aproximadamente equivalente a sacar la raíz cuadrada.
1
¿Por qué
casi
?
En IEEE-754, el exponente real es
e - 127
.
2
Para dividir esto entre dos, necesitaríamos
e / 2 - 64
, pero la aproximación anterior solo nos da
e / 2 - 127
.
Por lo tanto, necesitamos agregar 63 al exponente resultante.
Esto es contribuido por los bits 30-23 de esa constante mágica (
0x1fbb4000
).
Me imagino que los bits restantes de la constante mágica se han elegido para minimizar el error máximo en todo el rango de mantisa, o algo así. Sin embargo, no está claro si se determinó de forma analítica, iterativa o heurística.
Vale la pena señalar que este enfoque es algo no portátil. Hace (al menos) los siguientes supuestos:
-
La plataforma utiliza IEEE-754 de precisión simple para
float. -
La endianidad de la representación
float. - Que no se verá afectado por un comportamiento indefinido debido al hecho de que este enfoque viola las reglas de alias estricto de C / C ++.
Por lo tanto, debe evitarse a menos que esté seguro de que proporciona un comportamiento predecible en su plataforma (y de hecho, proporciona una aceleración útil frente a
sqrtf
).
1. sqrt (a ^ b) = (a ^ b) ^ 0.5 = a ^ (b / 2)
2. Ver, por ejemplo, https://en.wikipedia.org/wiki/Single-precision_floating-point_format#Exponent_encoding