c++ - online - punto flotante metodos numericos

Comparación de punto flotante revisado

Question

Este tema ha aparecido muchas veces en StackOverflow, pero creo que esta es una nueva toma. Sí, he leído los artículos de Bruce Dawson y Lo que todo científico informático debe saber sobre la aritmética de punto flotante y esta buena respuesta

Como lo entiendo, en un sistema típico hay cuatro problemas básicos al comparar números de punto flotante para la igualdad:

1. Los cálculos de punto flotante no son exactos
2. Si ab es "pequeño" depende de la escala de a y b
3. Si ab es "pequeño" depende del tipo de a y b (por ejemplo, float, double, long double)
4. El punto flotante suele tener representaciones + -infinidas, NaN y desnormalizadas, cualquiera de las cuales puede interferir con una formulación ingenua.

Esta respuesta - aka. "El enfoque de Google" - parece ser popular. Se encarga de todos los casos difíciles. Y escala la comparación de manera muy precisa, verificando si dos valores están dentro de un número fijo de ULPs entre sí. Así, por ejemplo, un número muy grande compara "casi igual" al infinito.

Sin embargo:

• Es muy desordenado, en mi opinión.
• No es particularmente portátil, depende en gran medida de representaciones internas, utiliza una unión para leer los bits de un flotador, etc.
• Solo maneja la precisión simple y la precisión doble IEEE 754 (en particular, no x86 long double)

Quiero algo similar, pero usando C ++ estándar y manejo de dobles largos. Por "estándar", me refiero a C ++ 03 si es posible y C ++ 11 si es necesario.

Aquí está mi intento.

#include <cmath>
#include <limits>
#include <algorithm>
namespace {
// Local version of frexp() that handles infinities specially.
template<typename T>
T my_frexp(const T num, int *exp)
{
typedef std::numeric_limits<T> limits;
// Treat +-infinity as +-(2^max_exponent).
if (std::abs(num) > limits::max())
{
*exp = limits::max_exponent + 1;
return std::copysign(0.5, num);
}
else return std::frexp(num, exp);
}
}
template<typename T>
bool almostEqual(const T a, const T b, const unsigned ulps=4)
{
// Handle NaN.
if (std::isnan(a) || std::isnan(b))
return false;
typedef std::numeric_limits<T> limits;
// Handle very small and exactly equal values.
if (std::abs(a-b) <= ulps * limits::denorm_min())
return true;
// frexp() does the wrong thing for zero. But if we get this far
// and either number is zero, then the other is too big, so just
// handle that now.
if (a == 0 || b == 0)
return false;
// Break the numbers into significand and exponent, sorting them by
// exponent.
int min_exp, max_exp;
T min_frac = my_frexp(a, &min_exp);
T max_frac = my_frexp(b, &max_exp);
if (min_exp > max_exp)
{
std::swap(min_frac, max_frac);
std::swap(min_exp, max_exp);
}
// Convert the smaller to the scale of the larger by adjusting its
// significand.
const T scaled_min_frac = std::ldexp(min_frac, min_exp-max_exp);
// Since the significands are now in the same scale, and the larger
// is in the range [0.5, 1), 1 ulp is just epsilon/2.
return std::abs(max_frac-scaled_min_frac) <= ulps * limits::epsilon() / 2;
}

Afirmo que este código (a) maneja todos los casos relevantes, (b) hace lo mismo que la implementación de Google para IEEE-754 de precisión simple y doble, y (c) es C ++ perfectamente estándar.

Una o más de estas afirmaciones es casi ciertamente errónea. Aceptaré cualquier respuesta que demuestre eso, preferiblemente con una solución. Una buena respuesta debe incluir uno o más de:

• Entradas específicas que difieren en más de unidades de ulps en el último lugar, pero para las que esta función devuelve true (cuanto mayor sea la diferencia, mejor)
• Entradas específicas que difieren hasta un ulps unidades en el último lugar, pero para las que esta función devuelve falso (cuanto más pequeña sea la diferencia, mejor)
• Cualquier caso (s) que he perdido
• Cualquier forma en que este código se base en un comportamiento indefinido o saltos dependiendo del comportamiento definido por la implementación. (Si es posible, por favor cite una especificación relevante).
• Arreglos para cualquier problema que identifiques
• Cualquier forma de simplificar el código sin romperlo.

Tengo la intención de colocar una recompensa no trivial en esta pregunta.

Answer 1

"Casi igual" no es una buena función

4 no es un valor apropiado: la respuesta que señala a los estados "Por lo tanto, 4 debería ser suficiente para el uso normal", pero no contiene ninguna base para esa afirmación. De hecho, hay situaciones normales en las que los números calculados en punto flotante por diferentes medios pueden diferir en muchos ULP, aunque serían iguales si se calcularan con matemáticas exactas. Por lo tanto, no debe haber un valor predeterminado para la tolerancia; a cada usuario se le debe exigir que proporcione el suyo, con suerte basado en un análisis exhaustivo de su código.

Como ejemplo de por qué un valor predeterminado de 4 ULP es malo, considere 1./49*49-1 . El resultado matemáticamente exacto es 0, pero el resultado calculado (binario IEEE 754 de 64 bits) es -0x1p-53, un error que excede 1e307 ULP del resultado exacto y casi 1e16 ULP del resultado calculado.

A veces, ningún valor es apropiado: en algunos casos, la tolerancia no puede ser relativa a los valores que se comparan, ni una tolerancia relativa matemáticamente exacta ni una tolerancia ULP cuantificada. Por ejemplo, casi todos los valores de salida en una FFT se ven afectados por casi todos los valores de entrada, y el error en cualquier elemento está relacionado con la magnitud de otros elementos. A la rutina "casi igual" se le debe proporcionar un contexto adicional con información sobre el posible error.

"Casi igual" tiene malas propiedades matemáticas: Esto muestra una de las deficiencias de "casi igual": la escala cambia los resultados. El siguiente código imprime 1 y 0.

double x0 = 1.1;
double x1 = 1.1 + 3*0x1p-52;
std::cout << almostEqual(x0, x1) << "/n";
x0 *= .8;
x1 *= .8;
std::cout << almostEqual(x0, x1) << "/n";

Otra falla es que no es transitiva; almostEqual(a, b) y almostEqual(b, c) no implica almostEqual(a, c) .

Un error en casos extremos

almostEqual(1.f, 1.f/11, 0x745d17) devuelve incorrectamente 1.

1.f / 11 es 0x1.745d18p-4. Restar esto de 1 (que es 0x10p-4) produce 0xe.8ba2e8p-4. Dado que un ULP de 1 es 0x1p-23, es decir 0xe.8ba2e8p19 ULP = 0xe8ba2e.8 / 2 ULP (desplazado 20 bits y dividido por 2, generando 19 bits) = 0x745d17.4 ULP. Eso excede la tolerancia especificada de 0x745d17, por lo que la respuesta correcta sería 0.

Este error se debe al redondeo en max_frac-scaled_min_frac .

Un escape fácil de este problema es especificar que las ulps deben ser menores que .5/limits::epsilon . Luego, el redondeo se produce en max_frac-scaled_min_frac solo si la diferencia (incluso cuando se redondea) supera las ulps ; si la diferencia es menor que eso, la resta es exacta, por Sterbenz ''Lemma.

Hubo una sugerencia sobre el uso del long double para corregir esto. Sin embargo, el long double no corregiría esto. Considere la posibilidad de comparar 1 y -0x1p-149f con los valores máximos establecidos en 1 / limits :: epsilon. A menos que su significand tenga 149 bits, el resultado de la resta se redondea a 1, que es menor o igual a 1 / limits :: epsilon ULP. Sin embargo, la diferencia matemática claramente supera 1.

Nota menor

El factor * limits::epsilon / 2 expresión factor * limits::epsilon / 2 convierte factor al tipo de punto flotante, lo que provoca errores de redondeo para valores grandes de factor que no son exactamente representables. Probablemente, la rutina no está diseñada para ser usada con valores tan grandes (millones de ULP en flotación), por lo que debe especificarse como un límite en la rutina en lugar de un error.

Answer 2

PRUEBAS FUTURAS : si alguna vez desea extender dicha comparación a float decimal http://en.wikipedia.org/wiki/Decimal64_floating-point_format en el futuro, y suponiendo que ldexp () y frexp () manejarán ese tipo con la raíz correcta , luego con voz estricta, 0.5 en return std::copysign(0.5, num); debe reemplazarse por T(1)/limits::radix() - o std::ldexp(T(1),-1) o algo ... (No pude encontrar una constante conveniente en std :: numeric_limits)

EDITAR Como comentó Nemo, las suposiciones de que ldexp y frexp usarían el FLOAT_RADIX correcto son falsas, se quedan con 2 ...

Por lo tanto, una versión portátil de prueba de futuro también debe usar

• std::scalbn(x,n) lugar de std::ldexp(x,n)

• exp=std::ilogb(std::abs(x)),y=std::scalbn(x,-exp) lugar de y=frexp(x,&exp)

• ahora que arriba y en es [1, FLOAT_RADIX) en lugar de [T (1) / Float_Radix, 1), devuelva copysign(T(1),num) lugar de 0.5 para el caso infinito de my_frexp, y pruebe los ulps*limits::epsilon() lugar de ulps * epsilon () / 2

Eso también requiere un estándar> = C ++ 11

Answer 3

Simplificación: puede evitar my_frexp descartando los casos no finitos primero todos juntos:

if( ! std::isfinite(a) || ! std::isfinite(b) )
return a == b;

Parece que isfinite está en C ++ 11 al menos

EDITAR Sin embargo, si la intención es tener limits::infinity() dentro de 1 ulp de limits::max()
entonces, la simplificación anterior no se mantiene, pero no debería my_frexp () devolver los limits::max_exponent+1 en *exp , en lugar de max_exponent + 2?