comparar floats o doubles

[swapd0]

Hola

Estoy haciendo una libreria de geometria, se supone que no debes poner cosas como a == b cuando a y b son floats o doubles. Asi que despues de buscar por google he hecho pruebas con la versión de comparar usando enteros (http://www.lomont.org/Math/Papers/2005/CompareFloat.pdf).

Pero tengo el siguiente problema. Con estos valores tan proximos la diferencia (diff) da un valor muy grande, mientras que en otros casos (normalmente cuando un valor no es 0.0) parece que va bien.

a = 1.2831915740410711e-14
b = 0.0

diff vale 4399142829531438451

¿Mando a la mierda esta forma de comparar y uso otra? ¿Algun otro truco para hacer que las comparaciones sean robustas? He pensado en "trasladar" al 0.0 antes de hacer la comparacion para ver si asi hay menos problemas de precision. El codigo seria b-=a; a=0.0;

Código [Seleccionar] Expand


static bool near_equal(double a, double b, boost::int64_t maxDiff = EPSILON_DIFF)
{
boost::int64_t ai = *reinterpret_cast<boost::int64_t *>(&a);
boost::int64_t bi = *reinterpret_cast<boost::int64_t *>(&b);
boost::int64_t sign = (((boost::uint64_t)(ai ^ bi)) >> 63) - 1;

BOOST_ASSERT( (sign == 0) || (sign == 0xffffffffffffffff) );

boost::int64_t diff = (((0x8000000000000000 - ai) & (~sign)) | (ai & sign)) - bi;
boost::int64_t v1 = maxDiff - diff;
boost::int64_t v2 = maxDiff + diff;

return (v1 | v2) >= 0;
}


[TrOnTxU]

Hola, voy a contestar pero no me hagas mucho caso porque no piloto mucho de IEEE 754 ^^

Con lo poco que entiendo del tema, lo de restar y compara con cero: NOOOOOO! (mira "Infernal zero" y "Catastrophic cancellation" en el artículo que pongo más abajo).

Creo que el problema lo tienes en el mismo algoritmo. Con flotantes de simple precision no te deberia funcionar tampoco cerca de cero.

La razón es que es una "optimización" de la propuesta de Bruce Dawson (Valve Software) para comparar flotantes basandose en la ULP (Unit in the Last Position). Que no tendria que ver en mi opinión con epsilon, como pones de constante. (aunque realmente no se que vale esa epsilon).

El mismo Bruce Dawson expone los problemas de su tan "conocida" tecnica o aproximación: http://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/.
Aunque yo te recomiendo que te leas todos sus articulos sobre flotantes si estas indagando sobre ello.

Mi "corrección" con unos unittest comparandolo con el "original" esta en floats:

Código [Seleccionar] Expand


#define FLT_EPSILON_REDEF     1.192092896e-07F        /* smallest such that 1.0+FLT_EPSILON != 1.0 */

// Dawson "optimizado"
static bool is_near_equal_float(float a, float b, I32 maxDiff = 1)
{
I32 ai = *((I32*)(&a));
I32 bi = *((I32*)(&b));
I32 sign = (((I32)(ai ^ bi)) >> 31) - 1;

TY_ASSERT( (sign == 0) || (sign == 0xffffffffffffffff) );

I32 diff = (((0x80000000 - ai) & (~sign)) | (ai & sign)) - bi;
I32 v1 = maxDiff - diff;
I32 v2 = maxDiff + diff;

return (v1 | v2) >= 0;
}

// Nuevo Dawson, no tan "optimizado", pero funciona cerca de cero
static bool is_near_equal_float_new_dawson(float a, float b, I32 maxUlpsDiff = 1, float maxRelDiff = FLT_EPSILON_REDEF )
{
// Check if the numbers are really close -- needed
// when comparing numbers near zero.
float absDiff = fabs(a - b);
if (absDiff <= maxRelDiff)
return true;

I32 ai = *((I32*)(&a));
I32 bi = *((I32*)(&b));
I32 sign = (((I32)(ai ^ bi)) >> 31) - 1;

TY_ASSERT( (sign == 0) || (sign == 0xffffffffffffffff) );

I32 diff = (((0x80000000 - ai) & (~sign)) | (ai & sign)) - bi;
I32 v1 = maxUlpsDiff - diff;
I32 v2 = maxUlpsDiff + diff;

return (v1 | v2) >= 0;
}


TY_TEST( StratosFloatQuestion, GlobalTest )
{
// El mismo numero
TY_CHECK( is_near_equal_float(0.1f,0.1f) );
TY_CHECK( is_near_equal_float_new_dawson(0.1f,0.1f) );

// Lejos de cero .. no hay problema
TY_CHECK( is_near_equal_float(0.1f,0.100000000001f) );
TY_CHECK( is_near_equal_float_new_dawson(0.1f,0.100000000001f) );

// Más de lo mismo, mumeros más grandes, menos precision, ... pero sigue funcionando
TY_CHECK( is_near_equal_float(10000.1f,10000.100000001f) );
TY_CHECK( is_near_equal_float_new_dawson(10000.1f,10000.100000001f) );

// mmm ... cerca de cero ... si implementas algo de Dawson ... haz caso a Dawson
TY_CHECK( is_near_equal_float(0.000000000011f,0.00000000001f) == false );
TY_CHECK( is_near_equal_float_new_dawson(0.000000000011f,0.00000000001f) == true );

// Y con cero ya ni te cuento ^^ ... viva Dawson!!!!
TY_CHECK( is_near_equal_float(0.0f,0.00000000001f) == false );
TY_CHECK( is_near_equal_float_new_dawson(0.0f,0.00000000001f) == true );
}


Como ves he "tocado" algunas cosas para que fuera en mi framework, pero el trabajo de vuelta lo veo sencillo.
Podrias pasarlo a doubles y hacer las comparaciones haber si es el problema.

El problema aqui es el fabs y el if, si lo que quieres es eliminar branchy code no parece una buena solución, pero sabiendo que el 0x80000000  debiera ser una mascara para el signo igual se puede hacer otro tipo de "optimización".
De todas formas creo que tambien depende de si la arquitectura en la que estas trabajando tiene que realizar cambios de registros de flotantes a enteros, ... en fin que no me preguntes de como iria más rápido, tendria que probarlo en varias maquinas.


Lo que me lleva a mi pregunta: ¿para que utilizas comparaciones "optimizadas"? ¿lo necesitas en tiempo real? ¿es para preprocesar algo?

Bueno, en realidad son tres, lo siento ^^

Si es para runtime, prueba a utilizar el Dawson "original", pero ten cuidado con los numeros cerca de cero, de todas formas yo que tu me miraria el articulo. Como dice Dawson en él: "There is no silver bullet. You have to choose wisely."


Espero haberte ayudado, un saludo.


EDIT: despues del "tostón" de post he encontrado la repsuesta corta ^^:

Código [Seleccionar] Expand

// Check if the numbers are really close -- needed
// when comparing numbers near zero.
float absDiff = fabs(a - b);
if (absDiff <= maxRelDiff)
return true;


esto es lo que he cmabiado   ^_^'

[swapd0]

Ok, probare esta tarde a ver como va.

El EPSILON_DIFF en las primeras pruebas ponía 5, ahora vale 50.

El problema de la técnica de usar enteros es que (por lo que he visto por ahi) no hay forma de calcular bien este EPSILON_DIFF. Si lo vas a usar en una aplicación donde los números van a estar en un rango conocido pues lo pones fijo y a correr. Si el rango es desconocido y vas a calcular números pequeños con otros mas grandes pues ni idea.

[TrOnTxU]

Te remito al articulo que te he puesto antes, la "tecnica" original en la que esta basado tu algoritmo (el de los enteros como lo llamas) es de Dawson y el lo explica muy bien, pero repito lo dicho:

No hay una formula mágica!!!!

Si lo que necesitas es una funcion que "resuelva todo" y que encima este "optimizada" (como la optimizacion del pdf que pones) es muy dificil, por no decir imposible (ya que es lo que dicen los que saben del tema del IEEE, no yo ^^).

Lo que he puesto yo no es lo más rapido del mundo ni de lejos, ya que estas haciendo el fabs y la resta, con lo que la solucion podria parecer "redudante" con el algoritmo que tenias.

Lo que te ofrece respecto a la comparacion con un delta "float", es más precisión "relativa". Pero si quieres velocidad no es la solución, simplemente separa comparaciones cerca y lejos de cero. Con dos funciones por ejemplo.

Un saludo


EDIT: copio esto del articulo que no paro de mencionar porque me parece muy interesante y explicativo:

There is no silver bullet. You have to choose wisely.

+ If you are comparing against zero, then relative epsilons and ULPs based comparisons are usually meaningless. You'll need to use an absolute epsilon, whose value might be some small multiple of FLT_EPSILON and the inputs to your calculation. Maybe.

+ If you are comparing against a non-zero number then relative epsilons or ULPs based comparisons are probably what you want. You'll probably want some small multiple of FLT_EPSILON for your relative epsilon, or some small number of ULPs. An absolute epsilon could be used if you knew exactly what number you were comparing against.

+ If you are comparing two arbitrary numbers that could be zero or non-zero then you need the kitchen sink. Good luck and God speed.

EDIT2: Estaba repasando el algoritmo que he puesto para ver si podia elimnar el "branch code", y he encontrado que tenia una cast mal y no funcionaba con negativos ... sorry ^^. El cambio es esta linea:

Código [Seleccionar] Expand

I32 sign = (((I32)(ai ^ bi)) >> 31) - 1; // Parece que bien, pero NOOOO!!
por esta otra:

Código [Seleccionar] Expand

I32 sign = (((U32)(ai ^ bi)) >> 31) - 1; // Esta si :D

ups! ... chorrada pero importante

[swapd0]

Jope con las comparaciones con floats y doubles, es bastante mas complicado de lo que parece.

EDIT: despues del "tostón" de post he encontrado la repsuesta corta ^^:

Código [Seleccionar] Expand

// Check if the numbers are really close -- needed
// when comparing numbers near zero.
float absDiff = fabs(a - b);
if (absDiff <= maxRelDiff)
return true;

Parece que con esto va mucho mejor. La version anterior que tenia hacia la comparacion restando los dos numeros y comparandolos con un epsilon pero en casos donde las coordenadas eran mas grandes el epsilon no funcionaba bien por falta de precision.

De paso le añadiré el código que has puesto de prueba al mío para asegurarme de que todo va bien.

Muchas gracias.

[TrOnTxU]

Jope con las comparaciones con floats y doubles, es bastante mas complicado de lo que parece.

Como dicen los expertos: "No importa cuan complicada creas que es la matematica de punto flotante, siempre lo será más."

En fin ... que remedio, es lo que hay ^^