Cómo verificar las dependencias de los flotadores

Question

Quiero determinar (en C++) si un número flotante es el inverso multiplicativo de otro número flotante. El problema es que tengo que usar una tercera variable para hacerlo. Por ejemplo, este código:

float x=5,y=0.2; 
if(x==(1/y)) cout<<"They are the multiplicative inverse of eachother"<<endl; 
else cout<<"They are NOT the multiplicative inverse of eachother"<<endl; 

la salida: "no son ...", que es un error y este código:

float x=5,y=0.2,z; 
z=1/y; 
if(x==z) cout<<"They are the multiplicative inverse of eachother"<<endl; 
else cout<<"They are NOT the multiplicative inverse of eachother"<<endl; 

la salida: "que son ..." que está justo .
¿por qué está pasando esto?

Answer 1

El flotador de precisión Problema

Usted tiene dos problemas aquí, pero ambos vienen de la misma raíz

No se puede comparar con precisión los flotadores. No puedes restarlos o dividirlos con precisión. No puede contar nada precisamente para ellos. Cualquier operación con ellos podría (y casi siempre) traer algún error en el resultado. Incluso a=0.2f no es una operación precisa. Las razones más profundas de eso están muy bien explicadas por los autores de las otras respuestas aquí. _{(Mi agradecimiento y votos para ellos por eso.)}

Aquí viene su primer y más simple error. Nunca se debe, nunca se , Nunca, Nunca, NUNCA uso en ellos == o su equivalente en cualquier idioma. En lugar de a==b, use Abs(a-b)<HighestPossibleError en su lugar.

---

Pero este no es el único problema en su tarea.

Abs(1/y-x)<HighestPossibleError tampoco funcionará. Al menos, no funcionará con la suficiente frecuencia. ¿Por qué?

Tomemos el par x = 1000 e y = 0.001. Tomemos el error relativo "inicial" de y para 10 ^-6.

_{(Error relativo = error/valor).}

Los errores relativos de los valores se suman a la multiplicación y división.

1/y es aproximadamente 1000. Su error relativo es el mismo 10 ^-6. ("1" no tiene errores)

Eso hace un error absoluto = 1000 * 10 ^-6 = 0.001. Cuando resta x más tarde, ese error será todo lo que quede. (Se suman los errores absolutos al sumar y restar, y el error de x es insignificante.) Seguramente, no cuenta con errores tan grandes, seguramente HighestPossibleError se establecerá más bajo y su programa arrojaría un buen par de x, y

Entonces, las siguientes dos reglas para operaciones de flotante: trate de no dividir mayores valorador por uno menor y Dios lo salve de restar los valores cercanos después de eso.

Existen dos formas sencillas de solucionar este problema.

• Por fundador lo de x, y tiene el mayor valor abs y dividir 1 por la mayor uno y sólo más tarde para restar el menor uno.

• Si desea comparar 1/y against x, mientras se está trabajando todavía con letras, no los valores, y sus operaciones que no hay errores, se multiplican ambos lados de la comparación por y y tiene 1 against x*y. _{(Por lo general, debe verificar los signos en esa operación, pero aquí utilizamos valores abs, por lo tanto, está limpio.)} La comparación de resultados no tiene división en absoluto.

En un camino más corto:

1/y V x <=> y*(1/y) V x*y <=> 1 V x*y 

ya sabemos que la comparación como 1 against x*y debe hacerse de modo:

const float HighestPossibleError=1e-10; 
if(Abs(x*y-1.0)<HighestPossibleError){... 

Eso es todo.

---

P.S. Si realmente lo necesita todo en una sola línea, utilice:

if(Abs(x*y-1.0)<1e-10){... 

Pero es mal estilo. Yo no lo aconsejaría.

P.P.S. En su segundo ejemplo, el compilador optimiza el código para que configure z a 5 antes de ejecutar cualquier código. Entonces, verificar 5 contra 5 funciona incluso para carrozas.

Answer 2

Deberá definir con precisión lo que significa que dos aproximaciones sean inversas multiplicativas. De lo contrario, no sabrá qué es lo que se supone que debe probar.

0.2 no tiene una representación binaria exacta. Si almacena números que no tienen una representación exacta con precisión limitada, no obtendrá respuestas que sean exactamente correctas.

Lo mismo ocurre en decimal. Por ejemplo, 1/3 no tiene representación decimal exacta. Puede almacenarlo como .333333. Pero luego tienes un problema. ¿Son 3 y .333333 inversos multiplicativos? Si los multiplicas, obtienes .999999. Si desea que la respuesta sea "sí", deberá crear una prueba de inversos multiplicativos que no sea tan simple como multiplicar y probar la igualdad en 1.

Lo mismo ocurre con el binario.

Answer 3

El problema es que 0.2 no se puede representar exactamente en binario, debido a su expansión binaria tiene un número infinito de dígitos:

1/5: 0.0011001100110011001100110011001100110011... 

Esto es similar a cómo 1/3 no se puede representar exactamente en decimal. Desde x se almacena en una float que tiene un número finito de bits, estas cifras se consiguen cortar en algún momento, por ejemplo:

x: 0.0011001100110011001100110011001 

El problema surge porque las CPU a menudo utilizan una mayor precisión internamente, por lo que cuando 'acaba de calcular 1/y, el resultado tendrá más dígitos, y cuando carga x para compararlos, x se extenderá para que coincida con la precisión interna de la CPU.

1/y: 0.0011001100110011001100110011001100110011001100110011 
    x: 0.0011001100110011001100110011001000000000000000000000 

Así que cuando haces una comparación directa bit por bit, son diferentes.

En el segundo ejemplo, sin embargo, almacenar el resultado en una variable significa que se trunca antes de hacer la comparación, por lo que la comparación de ellos en esta precisión, son iguales:

x: 0.0011001100110011001100110011001 
    z: 0.0011001100110011001100110011001 

Muchos compiladores tienen interruptores le puede permitir que los valores intermedios se trunquen en cada paso para coherencia, sin embargo, el consejo habitual es evitar hacer comparaciones directas entre valores de punto flotante y en su lugar comprobar si difieren en menos de un valor épsilon, que es lo que Gangnus is suggesting.

Answer 4

Lo que llama la atención es que cualquiera que sea la regla de redondeo, espera que el resultado de las dos versiones sea el mismo (ya sea dos veces incorrecto o dos veces a la derecha).

Lo más probable es que en el primer caso tenga lugar una promoción a mayor precisión en los registros FPU al evaluar x == 1/y, mientras que z = 1/y almacena realmente el resultado de precisión simple.

Otros colaboradores han explicado por qué 5 == 1/0.2 puede fallar, no necesito repetir eso.

Answer 5

Las discusiones en otras respuestas son excelentes y por eso no repetiré ninguna de ellas, pero no hay ningún código. Aquí hay un poco de código para verificar realmente si un par de flotantes da exactamente 1.0 cuando se multiplican.

El código hace algunas suposiciones/afirmaciones (que normalmente se cumplen en la plataforma x86):
- float 's son binario de 32 bits (AKA single precision) IEEE-754
- ya sea int' s o long 's son de 32 bits (I decidió no depender de la disponibilidad de uint32_t)
- memcpy() copias flotadores a Ints/largos de tal manera que se convierte en 8873283.0f 0x4B076543 (es decir, cierto "orden de bits" se espera)

Una suposición adicional es esta :
- recibe los flotantes reales que se multiplicarían * (es decir,multiplicación de los flotadores no usaría valores de precisión más altos que el hardware de matemáticas/biblioteca puede utilizar internamente)

#include <stdio.h> 
#include <string.h> 
#include <limits.h> 
#include <assert.h> 

#define C_ASSERT(expr) extern char CAssertExtern[(expr)?1:-1] 

#if UINT_MAX >= 0xFFFFFFFF 
typedef unsigned int uint32; 
#else 
typedef unsigned long uint32; 
#endif 
typedef unsigned long long uint64; 

C_ASSERT(CHAR_BIT == 8); 
C_ASSERT(sizeof(uint32) == 4); 
C_ASSERT(sizeof(float) == 4); 

int ProductIsOne(float f1, float f2) 
{ 
    uint32 m1, m2; 
    int e1, e2, s1, s2; 
    int e; 
    uint64 m; 

    // Make sure floats are 32-bit IEE754 and 
    // reinterpreted as integers as we expect 
    { 
    static const float testf = 8873283.0f; 
    uint32 testi; 
    memcpy(&testi, &testf, sizeof(testf)); 
    assert(testi == 0x4B076543); 
    } 

    memcpy(&m1, &f1, sizeof(f1)); 
    s1 = m1 >= 0x80000000; 
    m1 &= 0x7FFFFFFF; 
    e1 = m1 >> 23; 
    m1 &= 0x7FFFFF; 
    if (e1 > 0) m1 |= 0x800000; 

    memcpy(&m2, &f2, sizeof(f2)); 
    s2 = m2 >= 0x80000000; 
    m2 &= 0x7FFFFFFF; 
    e2 = m2 >> 23; 
    m2 &= 0x7FFFFF; 
    if (e2 > 0) m2 |= 0x800000; 

    if (e1 == 0xFF || e2 == 0xFF || s1 != s2) // Inf, NaN, different signs 
    return 0; 

    m = (uint64)m1 * m2; 

    if (!m || (m & (m - 1))) // not a power of 2 
    return 0; 

    e = e1 + !e1 - 0x7F - 23 + e2 + !e2 - 0x7F - 23; 
    while (m > 1) m >>= 1, e++; 

    return e == 0; 
} 

const float testData[][2] = 
{ 
    { .1f, 10.0f }, 
    { 0.5f, 2.0f }, 
    { 0.25f, 2.0f }, 
    { 4.0f, 0.25f }, 
    { 0.33333333f, 3.0f }, 
    { 0.00000762939453125f, 131072.0f }, // 2^-17 * 2^17 
    { 1.26765060022822940E30f, 7.88860905221011805E-31f }, // 2^100 * 2^-100 
    { 5.87747175411143754E-39f, 1.70141183460469232E38f }, // 2^-127 (denormalized) * 2^127 
}; 

int main(void) 
{ 
    int i; 
    for (i = 0; i < sizeof(testData)/sizeof(testData[0]); i++) 
    printf("%g * %g %c= 1\n", 
      testData[i][0], testData[i][1], 
      "!="[ProductIsOne(testData[i][0], testData[i][1])]); 
    return 0; 
} 

de salida (ver en ideone.com):

0.1 * 10 != 1 
0.5 * 2 == 1 
0.25 * 2 != 1 
4 * 0.25 == 1 
0.333333 * 3 != 1 
7.62939e-06 * 131072 == 1 
1.26765e+30 * 7.88861e-31 == 1 
5.87747e-39 * 1.70141e+38 == 1