¿Qué comparación de punto flotante es más precisa y por qué?

Question

Estoy experimentando con diferentes implementaciones del método de Newton para calcular las raíces cuadradas. Una decisión importante es cuándo terminar el algoritmo.

Obviamente no va a hacer para utilizar la diferencia absoluta entre y*y y x donde y es la estimación actual de la raíz cuadrada de x, ya que para valores grandes de x puede que no sea posible representar su raíz cuadrada suficiente precisión.

Así que se supone que debo usar un criterio relativo. Ingenuamente hubiera usado algo como esto:

static int sqrt_good_enough(float x, float y) { 
    return fabsf(y*y - x)/x < EPS; 
} 

Y esto parece funcionar muy bien. Sin embargo, recientemente he empezado a leer Kernighan y los elementos del estilo de programación de Plauger y dan un programa Fortran para el mismo algoritmo en el capítulo 1, cuyos criterios de terminación, traducido en C, sería:

static int sqrt_good_enough(float x, float y) { 
    return fabsf(x/y - y) < EPS * y; 
} 

Tanto son matemáticamente equivalentes , pero ¿hay alguna razón para preferir una forma sobre la otra?

Answer 1

Los dos no son exactamente equivalentes matemáticamente, a menos que escriba fabsf (y * y - x)/(y * y) < EPS para el primero. (perdón por el error en mi comentario original)

Pero creo que el punto clave es hacer que la expresión aquí coincida con su fórmula para calcular y en la iteración de Newton. Por ejemplo, si su fórmula y es y = (y + x/y)/2, debe usar el estilo de Kernighan y Plauger. Si es y = (y * y + x)/(2 * y) debe usar (y * y - x)/(y * y) < EPS.

En general, los criterios de terminación deberían ser que abs (y (n + 1) - y (n)) es lo suficientemente pequeño (es decir, menor que y (n + 1) * EPS). Esta es la razón por la cual las dos expresiones deben coincidir. Si no coinciden exactamente, es posible que la prueba de terminación decida que el residuo no es lo suficientemente pequeño, mientras que la diferencia en y (n) es menor que el error de punto flotante, debido a diferentes escalas. El resultado sería un ciclo infinito porque y (n) ha dejado de cambiar y nunca se cumplen los criterios de terminación.

Por ejemplo, el siguiente código de Matlab es exactamente el mismo Newton solucionador como su primer ejemplo, pero corre siempre:

x = 6.800000000000002 
yprev = 0 
y = 2 
while abs(y*y - x) > eps*abs(y*y) 
    yprev = y; 
    y = 0.5*(y + x/y); 
end 

La versión de C/C++ de ella tiene el mismo problema.

Answer 2

Todavía no son equivalentes; el de abajo es matemáticamente equivalente a fabsf(y*y - x)/(y*y) < EPS. El problema que veo con los suyos es que si y*y se desborda (probablemente porque x es FLT_MAX y y se elige desafortunadamente), entonces la terminación puede que nunca ocurra. La siguiente interacción usa dobles.

>>> import math 
>>> x = (2.0 - 2.0 ** -52) * 2.0 ** 1023 
>>> y = x/math.sqrt(x) 
>>> y * y - x 
inf 
>>> y == 0.5 * (y + x/y) 
True 

EDIT: como un comentario (ahora suprimido) señaló, también es bueno para compartir las operaciones entre la iteración y la prueba de terminación.

EDIT2: ambos probablemente tengan problemas con el subnormal x. El professionals normaliza x para evitar las complicaciones de ambos extremos.