resultados inesperados con modo de redondeo específica

Question

En primer lugar, lo siento por el mal título:/

estoy tratando de reproducir los resultados de un trabajo sobre el cálculo de los valores propios de una matriz simétrica tridiagonales. Estoy determinando los límites superiores e inferiores de algunos valores redondeando a más y menos a infinito, respectivamente.

En lugar de cambiar el modo de redondeo cada vez, solo uso el 'truco': fl⁻ (y) = -fl⁺ (-y), donde fl⁻ (y) es el valor de y cuando se usa el signo menos modo de redondeo infinito y fl⁺ (y) es el valor de y cuando se usa el modo de redondeo a más infinito. lo tanto, tengo el siguiente fragmento de código en C:

fesetround(FE_UPWARD); 
first = - (-d[i] + x); 
second = (- ((e[i-1]*e[i-1])/a_inf)); 
a_inf = first + second; 

first = d[i] - x; 
second = - ((e[i-1]*e[i-1])/a_sup); 
a_sup = first + second; 

y trabaja muy bien a excepción de un ejemplo en el que A_INF me da el resultado correcto, pero a_sup da un resultado erróneo, aunque ambas primera y segunda variables parece tener los mismos valores

Sin embargo, si lo hago de esta manera:

fesetround(FE_UPWARD); 
    first = - (-d[i] + x); 
    second = (- ((e[i-1]*e[i-1])/a_inf)); 

    fesetround(FE_DOWNWARD); 
    first = - (-d[i] + x); 
    second = (- ((e[i-1]*e[i-1])/a_sup)); 

puedo obtener los resultados correctos. Entonces, si uso el truco fl⁻ (y) = -fl⁺ (-y), obtengo los resultados correctos, si cambio el modo de redondeo y uso la expresión original obtengo resultados incorrectos. ¿Alguna idea de por qué?

En ambos casos, las variables primero y segundo valores son los siguientes:

first 1.031250000000000e+07, second -1.031250000000000e+07 
first 1.031250000000000e+07, second -1.031250000000000e+07 

Y los valores correctos para A_INF y a_sup son -1.862645149230957e-09 y + 1.862645149230957e-09, respectivamente, pero en el primer caso a_sup = 0, lo que está mal

lo que supongo que está sucediendo es una especie de cancelación catastrófica, pero no tengo ni idea de cómo resolverlo en este caso ...

Gracias de ¡avanzar!

Answer 1

El primer problema que tiene es que solo está utilizando el 'truco' para obtener el redondeo de first hacia -inf, y no second, por lo que todavía se redondea hacia + inf.

El segundo problema es que C no ofrece ningún tipo de garantía acerca de cómo realiza los cálculos en coma flotante. En particular, el compilador es libre de reordenar y volver a asociar las operaciones como lo considere apropiado para el rendimiento, incluso si tales reordenamientos pueden cambiar el comportamiento de redondeo del programa. Así, cuando decimos:

first = - (-d[i] + x); 

el compilador puede reorganizar y hacer que un solo restar, en lugar de negar, agregar, anular, que invierte la dirección del redondeo. Ahora, a veces puede hacer que las cosas funcionen de la manera que espera al deshabilitar toda la optimización, pero incluso con eso, no hay garantía.