Compare doble a cero usando epsilon

Hoy, estaba buscando un código C ++ (escrito por otra persona) y encontré esta sección:

double someValue = ...
if (someValue <  std::numeric_limits<double>::epsilon() && 
    someValue > -std::numeric_limits<double>::epsilon()) {
  someValue = 0.0;
}

Estoy tratando de averiguar si esto tiene sentido.

La documentación para epsilon()dice:

La función devuelve la diferencia entre 1 y el valor más pequeño mayor que 1 que es representable [por un doble].

¿Esto se aplica también a 0, epsilon()es decir, el valor más pequeño es mayor que 0? ¿O hay números entre 0y 0 + epsilonque pueden ser representados por a double?

Si no, ¿no es la comparación equivalente a someValue == 0.0?

Suponiendo que IEEE doble de 64 bits, hay una mantisa de 52 bits y un exponente de 11 bits. Vamos a romperlo en pedazos:

1.0000 00000000 00000000 00000000 00000000 00000000 00000000 × 2^0 = 1

El número representable más pequeño mayor que 1:

1.0000 00000000 00000000 00000000 00000000 00000000 00000001 × 2^0 = 1 + 2^-52

Por lo tanto:

epsilon = (1 + 2^-52) - 1 = 2^-52

¿Hay algún número entre 0 y épsilon? Mucho ... Por ejemplo, el número mínimo positivo representable (normal) es:

1.0000 00000000 00000000 00000000 00000000 00000000 00000000 × 2^-1022 = 2^-1022

De hecho, hay (1022 - 52 + 1)×2^52 = 4372995238176751616números entre 0 y épsilon, que es el 47% de todos los números representables positivos ...

La prueba ciertamente no es lo mismo que someValue == 0. La idea de los números de coma flotante es que almacenan un exponente y un significado. Por lo tanto, representan un valor con un cierto número de cifras binarias significativas de precisión (53 en el caso de un IEEE doble). Los valores representables están mucho más densos cerca de 0 que cerca de 1.

Para usar un sistema decimal más familiar, suponga que almacena un valor decimal "a 4 cifras significativas" con exponente. Entonces el siguiente valor representable mayor que 1es 1.001 * 10^0, y epsilones 1.000 * 10^-3. Pero 1.000 * 10^-4también es representable, suponiendo que el exponente puede almacenar -4. Puede creer que un doble IEEE puede almacenar exponentes menores que el exponente de epsilon.

No puede distinguir solo con este código si tiene sentido o no usar epsilonespecíficamente como límite, debe mirar el contexto. Puede ser que epsilonsea ​​una estimación razonable del error en el cálculo que produjo someValue, y puede ser que no lo sea.

Hay números que existen entre 0 y épsilon porque épsilon es la diferencia entre 1 y el siguiente número más alto que se puede representar por encima de 1 y no la diferencia entre 0 y el siguiente número más alto que se puede representar por encima de 0 (si lo fuera, que el código haría muy poco): -

#include <limits>

int main ()
{
  struct Doubles
  {
      double one;
      double epsilon;
      double half_epsilon;
  } values;

  values.one = 1.0;
  values.epsilon = std::numeric_limits<double>::epsilon();
  values.half_epsilon = values.epsilon / 2.0;
}

Usando un depurador, detenga el programa al final de main y mire los resultados y verá que epsilon / 2 es distinto de epsilon, cero y uno.

Entonces, esta función toma valores entre +/- epsilon y los hace cero.

Con el siguiente programa se puede imprimir una aproximación de épsilon (diferencia más pequeña posible) alrededor de un número (1.0, 0.0, ...). Imprime el siguiente resultado:
epsilon for 0.0 is 4.940656e-324
epsilon for 1.0 is 2.220446e-16
un poco de reflexión deja en claro que el épsilon se hace más pequeño cuanto más pequeño es el número que usamos para observar su valor épsilon, porque el exponente puede ajustarse al tamaño de ese número.

#include <stdio.h>
#include <assert.h>
double getEps (double m) {
  double approx=1.0;
  double lastApprox=0.0;
  while (m+approx!=m) {
    lastApprox=approx;
    approx/=2.0;
  }
  assert (lastApprox!=0);
  return lastApprox;
}
int main () {
  printf ("epsilon for 0.0 is %e\n", getEps (0.0));
  printf ("epsilon for 1.0 is %e\n", getEps (1.0));
  return 0;
}

Supongamos que estamos trabajando con números de punto flotante de juguete que caben en un registro de 16 bits. Hay un bit de signo, un exponente de 5 bits y una mantisa de 10 bits.

El valor de este número de coma flotante es la mantisa, interpretada como un valor decimal binario, multiplicado por dos por la potencia del exponente.

Alrededor de 1 el exponente es igual a cero. Entonces, el dígito más pequeño de la mantisa es una parte en 1024.

Cerca de la mitad del exponente es menos uno, por lo que la parte más pequeña de la mantisa es la mitad de grande. Con un exponente de cinco bits puede alcanzar 16 negativos, en cuyo punto la parte más pequeña de la mantisa vale una parte en 32 m. Y en el exponente negativo 16, el valor es de alrededor de una parte en 32k, ¡mucho más cercano a cero que el épsilon alrededor de uno que calculamos arriba!

Ahora, este es un modelo de punto flotante de juguete que no refleja todas las peculiaridades de un sistema de punto flotante real, pero la capacidad de reflejar valores más pequeños que épsilon es razonablemente similar con los valores de punto flotante real.

La diferencia entre Xy el siguiente valor de Xvaría según X.
epsilon()es solo la diferencia entre 1y el siguiente valor de 1.
La diferencia entre 0y el siguiente valor de 0no es epsilon().

En su lugar, puede usar std::nextafterpara comparar un valor doble con 0el siguiente:

bool same(double a, double b)
{
  return std::nextafter(a, std::numeric_limits<double>::lowest()) <= b
    && std::nextafter(a, std::numeric_limits<double>::max()) >= b;
}

double someValue = ...
if (same (someValue, 0.0)) {
  someValue = 0.0;
}

Creo que eso depende de la precisión de su computadora. Eche un vistazo a esta tabla : puede ver que si su épsilon está representado por el doble, pero su precisión es mayor, la comparación no es equivalente a

someValue == 0.0

Buena pregunta de todos modos!

No puede aplicar esto a 0, debido a las partes de mantisa y exponente. Debido al exponente, puede almacenar números muy pequeños, que son más pequeños que épsilon, pero cuando intenta hacer algo como (1.0 - "número muy pequeño") obtendrá 1.0. Epsilon es un indicador no de valor, sino de precisión de valor, que está en mantisa. Muestra cuántos dígitos decimales consecuentes correctos de número podemos almacenar.

Con el punto flotante IEEE, entre el valor positivo distinto de cero más pequeño y el valor negativo distinto de cero más pequeño, existen dos valores: cero positivo y cero negativo. Probar si un valor está entre los valores más pequeños distintos de cero es equivalente a probar la igualdad con cero; Sin embargo, la asignación puede tener un efecto, ya que cambiaría un cero negativo a un cero positivo.

Sería concebible que un formato de punto flotante pudiera tener tres valores entre los valores positivos y negativos finitos más pequeños: infinitesimal positivo, cero sin signo e infinitesimal negativo. No estoy familiarizado con ningún formato de punto flotante que de hecho funcione de esa manera, pero tal comportamiento sería perfectamente razonable y posiblemente mejor que el de IEEE (quizás no sea lo suficientemente mejor como para que valga la pena agregar hardware adicional para admitirlo, pero matemáticamente 1) / (1 / INF), 1 / (- 1 / INF) y 1 / (1-1) deben representar tres casos distintos que ilustran tres ceros diferentes). No sé si algún estándar C exigiría que los infinitesimales firmados, si existen, tendrían que comparar igual a cero. Si no lo hacen, un código como el anterior podría garantizar que, por ejemplo,

Entonces, digamos que el sistema no puede distinguir 1.000000000000000000000 y 1.000000000000000000001. eso es 1.0 y 1.0 + 1e-20. ¿Crees que todavía hay algunos valores que se pueden representar entre -1e-20 y + 1e-20?

Además, una buena razón para tener una función de este tipo es eliminar "denormals" (esos números muy pequeños que ya no pueden usar el "1" inicial implícito y tienen una representación FP especial). Por qué querrías hacer esto? Debido a que algunas máquinas (en particular, algunas Pentium 4s más antiguas) se vuelven muy, muy lentas al procesar denormals. Otros se vuelven un poco más lentos. Si su aplicación realmente no necesita estos números muy pequeños, vaciarlos a cero es una buena solución. Un buen lugar para considerar esto son los últimos pasos de cualquier filtro IIR o función de descomposición.

Ver también: ¿Por qué cambiar 0.1f a 0 ralentiza el rendimiento en 10x?

y http://en.wikipedia.org/wiki/Denormal_number