Comparación relativa de números de coma flotante

Tengo una función numérica que f(x, y)devuelve un número de coma flotante doble que implementa alguna fórmula y quiero verificar que sea correcta contra las expresiones analíticas para todas las combinaciones de los parámetros xy yque me interese. ¿Cuál es la forma correcta de comparar lo calculado y números analíticos de coma flotante?

Digamos que los dos números son ay b. Hasta ahora me he asegurado de que tanto los errores absolutos ( abs(a-b) < eps) como los relativos ( abs(a-b)/max(abs(a), abs(b)) < eps) sean menores que eps. De esa forma, detectará imprecisiones numéricas incluso si los números son, digamos, alrededor de 1e-20.

Sin embargo, hoy descubrí un problema, el valor numérico ay el valor analítico bfueron:

In [47]: a                                                                     
Out[47]: 5.9781943146790832e-322

In [48]: b                                                                     
Out[48]: 6.0276008792632078e-322

In [50]: abs(a-b)                                                              
Out[50]: 4.9406564584124654e-324

In [52]: abs(a-b) / max(a, b)                                                  
Out[52]: 0.0081967213114754103

Entonces, el error absoluto [50] es (obviamente) pequeño, pero el error relativo [52] es grande. Entonces pensé que tenía un error en mi programa. Al depurar, me di cuenta, que estos números son normales . Como tal, escribí la siguiente rutina para hacer la comparación relativa adecuada:

real(dp) elemental function rel_error(a, b) result(r)
real(dp), intent(in) :: a, b
real(dp) :: m, d
d = abs(a-b)
m = max(abs(a), abs(b))
if (d < tiny(1._dp)) then
    r = 0
else
    r = d / m
end if
end function

Donde tiny(1._dp)devuelve 2.22507385850720138E-308 en mi computadora. Ahora todo funciona y simplemente obtengo 0 como error relativo y todo está bien. En particular, el error relativo anterior [52] es incorrecto, simplemente es causado por una precisión insuficiente de los números denormales. ¿Es rel_errorcorrecta mi implementación de la función? ¿Debo comprobar que abs(a-b)es menor que minúsculo (= denormal) y devolver 0? ¿O debería comprobar alguna otra combinación, como max(abs(a), abs(b))?

Me gustaría saber cuál es la forma "adecuada".

Puede comprobar directamente para denormals utilizando isnormal()partir de math.h(C99 o temprano, POSIX.1 o posterior). En Fortran, si el módulo ieee_arithmeticestá disponible, puede usarlo ieee_is_normal(). Para ser más precisos acerca de la igualdad difusa, debe considerar la representación en coma flotante de denormales y decidir qué quiere decir para que los resultados sean lo suficientemente buenos.

Más aún, para creer que cualquiera de los resultados es correcto, debe asegurarse de no perder demasiados dígitos en un paso intermedio. La computación con denormales generalmente no es confiable y debe evitarse haciendo que su algoritmo se vuelva a escalar internamente. Para asegurarse de que su escala interna fue exitosa, le recomiendo activar las excepciones de punto flotante usando feenableexcept()C99 o el ieee_arithmeticmódulo en Fortran.

Aunque puede hacer que su aplicación capte la señal que se genera en excepciones de coma flotante, todos los núcleos que he intentado restablecen el indicador de hardware para fetestexcept()que no devuelva un resultado útil. Cuando se ejecuta con -fp_trap, los programas PETSc (por defecto) imprimirán un seguimiento de la pila cuando se genera un error de coma flotante, pero no identificarán la línea ofensiva. Si se ejecuta en un depurador, el depurador conserva el indicador de hardware y rompe la expresión ofensiva. Puede verificar el motivo exacto llamando fetestexceptdesde el depurador donde el resultado es un OR bit a bit de los siguientes indicadores (los valores pueden variar según la máquina, vea fenv.h; estos valores son para x86-64 con glibc).

• FE_INVALID = 0x1
• FE_DIVBYZERO = 0x4
• FE_OVERFLOW = 0x8
• FE_UNDERFLOW = 0x10
• FE_INEXACT = 0x20

Donald Knuth tiene una propuesta para un algoritmo de comparación de coma flotante en el volumen 2 "Algoritmos seminéricos" de "El arte de la programación de computadoras". Fue implementado en C por Th. Belding (ver paquete fcmp ) y está disponible en GSL .

Los números desnormalizados óptimamente redondeados pueden tener un alto error relativo. (Vaciar eso a cero mientras todavía lo llama un error relativo es engañoso).

Pero cerca de cero, calcular los errores relativos no tiene sentido.

Por lo tanto, incluso antes de alcanzar números desnormalizados, probablemente debería cambiar a una precisión absoluta (es decir, la que desea garantizar en este caso).

$y$$x$$|y-x|\le absacc+relacc*\max(|x|,|y|)$

Entonces, los usuarios de su código saben con precisión cuánta precisión tienen realmente.