Resultado de punto flotante diferente con optimización habilitada - ¿error del compilador?

La salida debe ser: 4.5 4.6 Eso es lo que la salida sería si tuviera una precisión infinita, o si estuviera trabajando con un dispositivo que utiliza una representación de coma flotante basada en decimales en lugar de basada en binarios. La mayoría de las computadoras usan el estándar binario de coma flotante IEEE.

Como Maxim Yegorushkin ya señaló en su respuesta, parte del problema es que internamente su computadora está utilizando una representación de punto flotante de 80 bits. Esto es sólo parte del problema, sin embargo. La base del problema es que cualquier número de la forma n. nn5 no tiene una representación flotante binaria exacta. Esos casos de esquina son siempre números inexactos.

Si realmente desea que su redondeo sea capaz de redondear de manera confiable estos casos de esquina, necesita un algoritmo de redondeo que aborde el hecho de que n.n5, n.nn5, o n.nnn5, etc. (pero no n.5) es siempre inexacta. Encuentre el caso de esquina que determina si algún valor de entrada redondea hacia arriba o hacia abajo y devuelve el valor redondeado hacia arriba o hacia abajo basado en una comparación con este caso de esquina. Y usted necesita tener cuidado de que un compilador de optimización no pondrá ese caso de esquina encontrado en un registro de precisión extendido.

Ver ¿Cómo Redondea Excel con éxito los números flotantes a pesar de que son imprecisos? para tal algoritmo.

O simplemente puede vivir con el hecho de que los casos de esquina a veces redondearán erróneamente.

Diferentes compiladores tienen diferentes configuraciones de optimización. Algunas de esas configuraciones de optimización más rápidas no mantienen reglas estrictas de coma flotante de acuerdo con IEEE 754. Visual Studio tiene una configuración específica, /fp:strict, /fp:precise, /fp:fast, donde /fp:fast viola el estándar sobre lo que se puede hacer. Puede encontrar que este indicador es lo que controla la optimización en tales configuraciones. También puede encontrar una configuración similar en GCC que cambia el comportamiento.

Si este es el caso entonces lo único que es diferente entre los compiladores es que GCC buscaría el comportamiento de coma flotante más rápido por defecto en optimizaciones más altas, mientras que Visual Studio no cambia el comportamiento de coma flotante con niveles de optimización más altos. Por lo tanto, podría no ser necesariamente un error real, sino el comportamiento previsto de una opción que no sabía que estaba encendiendo.

Para aquellos que no pueden reproducir el error: no descomentar los stmt de depuración comentados, afectan el resultado.

Esto implica que el problema está relacionado con las sentencias debug. Y parece que hay un error de redondeo causado por cargar los valores en registros durante las sentencias de salida, por lo que otros encontraron que puede arreglar esto con -ffloat-store

Otra pregunta:

Me preguntaba, ¿debería encender siempre -ffloat-store opción?

Para ser frívolo, debe haber una razón por la que algunos programadores no activan -ffloat-store, de lo contrario la opción no existiría (del mismo modo, debe haber una razón por la que algunos programadores activan). No recomendaría siempre encenderlo o siempre apagarlo. Activarlo evita algunas optimizaciones, pero desactivarlo permite el tipo de comportamiento que está obteniendo.

Pero, generalmente, hay algún desajuste entre binarios números de coma flotante (como los que usa la computadora) y números decimales de coma flotante (con los que la gente está familiarizada), y ese desajuste puede causar un comportamiento similar al que obtiene (para ser claros, el comportamiento que está obteniendo es no causado por este desajuste, pero un comportamiento similar puede ser). La cosa es, ya que usted ya tiene cierta vaguedad cuando se trata de punto flotante, no puedo decir que -ffloat-store lo hace mejor o peor.

En su lugar, usted puede quiero buscar en otras soluciones para el problema que está tratando de resolver (desafortunadamente, Koenig no apunta al papel real, y realmente no puedo encontrar un lugar "canónico" obvio para él, así que tendré que enviarte a Google ).

Si no estás redondeando para fines de salida, probablemente miraría std::modf() (en cmath) y std::numeric_limits<double>::epsilon() (en limits). Pensando en la función original round(), creo que sería más limpio reemplazar la llamada a std::floor(d + .5) con una llamada para esta función:

// this still has the same problems as the original rounding function
int round_up(double d)
{
    // return value will be coerced to int, and truncated as expected
    // you can then assign the int to a double, if desired
    return d + 0.5;
}

Creo que eso sugiere la siguiente mejora:

// this won't work for negative d ...
// this may still round some numbers up when they should be rounded down
int round_up(double d)
{
    double floor;
    d = std::modf(d, &floor);
    return floor + (d + .5 + std::numeric_limits<double>::epsilon());
}

Una nota simple: std::numeric_limits<T>::epsilon() se define como "el número más pequeño añadido a 1 que crea un número no igual a 1."Por lo general, es necesario utilizar un epsilon relativo (es decir, escalar epsilon de alguna manera para tener en cuenta el hecho de que está trabajando con números que no sean "1"). La suma de d, .5 y std::numeric_limits<double>::epsilon() debe estar cerca de 1, por lo que agrupar esa suma significa que std::numeric_limits<double>::epsilon() será aproximadamente el tamaño correcto para lo que lo estamos haciendo. En todo caso, std::numeric_limits<double>::epsilon() será demasiado grande (cuando la suma de los tres es menor que uno) y puede hacer que redondeemos algunos números hacia arriba cuando no deberíamos.

Hoy en día, usted debe considerar std::nearbyint().

La respuesta aceptada es correcta si está compilando para un objetivo x86 que no incluye SSE2. Todos los procesadores x86 modernos admiten SSE2, por lo que si puede aprovecharlo, debe:

-mfpmath=sse -msse2 -ffp-contract=off

Vamos a desglosar esto.

-mfpmath=sse -msse2. Esto realiza el redondeo mediante el uso de registros SSE2, que es mucho más rápido que almacenar todos los resultados intermedios en la memoria. Tenga en cuenta que este es ya el predeterminado en GCC para x86-64. De la wiki del CCG :

En más los procesadores x86 modernos que admiten SSE2, especificando las opciones del compilador -mfpmath=sse -msse2, aseguran que todas las operaciones float y double se realicen en registros SSE y se redondeen correctamente. Estas opciones no afectan al ICB y, por lo tanto, deben utilizarse siempre que sea posible para obtener resultados numéricos predecibles.

-ffp-contract=off. Sin embargo, controlar el redondeo no es suficiente para una coincidencia exacta. Las instrucciones FMA (multiplicación-suma fusionada) pueden cambiar el comportamiento de redondeo frente a sus contrapartes no fusionadas, por lo que necesitamos deshabilitarlo. Este es el valor predeterminado en Clang, no GCC. Como se explica en esta respuesta :

Un FMA tiene solo un redondeo (efectivamente mantiene una precisión infinita para el resultado interno de multiplicación temporal), mientras que un ADD + MUL tiene dos.

Al deshabilitar FMA, obtenemos resultados que coinciden exactamente con la depuración y la liberación, a costa de cierto rendimiento (y precisión). Todavía podemos aprovechar otros beneficios de rendimiento de SSE y AVX.

Personalmente, he tenido el mismo problema yendo hacia el otro lado - de gcc a VS. En la mayoría de los casos creo que es mejor evitar la optimización. La única vez que vale la pena es cuando se trata de métodos numéricos que implican grandes matrices de datos en coma flotante. Incluso después de desensamblar, a menudo me decepcionan las opciones de los compiladores. Muy a menudo es más fácil usar intrínsecos del compilador o simplemente escribir el ensamblaje usted mismo.

Profundicé más en este problema y puedo aportar más precisiones. Primero, las representaciones exactas de 4.45 y 4.55 según gcc en x84_64 son las siguientes (con libquadmath para imprimir la última precisión):

float 32:   4.44999980926513671875
double 64:  4.45000000000000017763568394002504646778106689453125
doublex 80: 4.449999999999999999826527652402319290558807551860809326171875
quad 128:   4.45000000000000000000000000000000015407439555097886824447823540679418548304813185723105561919510364532470703125

float 32:   4.55000019073486328125
double 64:  4.54999999999999982236431605997495353221893310546875
doublex 80: 4.550000000000000000173472347597680709441192448139190673828125
quad 128:   4.54999999999999999999999999999999984592560444902113175552176459320581451695186814276894438080489635467529296875

Como Maxim dijo anteriormente, el problema se debe al tamaño de 80 bits de los registros FPU.

Pero ¿por qué el problema nunca ocurre en Windows? en IA-32, la FPU x87 fue configurada para usar una precisión interna para la mantissa de 53 bits (equivalente hasta un tamaño total de 64 bits: double). Para Linux y Mac OS, se utilizó la precisión predeterminada de 64 bits (equivalente a un tamaño total de 80 bits: long double). Por lo tanto, el problema debería ser posible, o no, en estas diferentes plataformas cambiando la palabra de control de la FPU (suponiendo que la secuencia de instrucciones desencadenaría el error). El problema se informó a gcc como bug 323 (lea al menos el comentario 92! ).

Para mostrar la precisión de mantissa en Windows, puede compilar esto en 32 bits con VC++:

#include "stdafx.h"
#include <stdio.h>  
#include <float.h>  

int main(void)
{
    char t[] = { 64, 53, 24, -1 };
    unsigned int cw = _control87(0, 0);
    printf("mantissa is %d bits\n", t[(cw >> 16) & 3]);
}

Y en Linux/Cygwin:

#include <stdio.h>

int main(int argc, char **argv)
{
    char t[] = { 24, -1, 53, 64 };
    unsigned int cw = 0;
    __asm__ __volatile__ ("fnstcw %0" : "=m" (*&cw));
    printf("mantissa is %d bits\n", t[(cw >> 8) & 3]);
}

Tenga en cuenta que con gcc puede establecer la precisión de FPU con -mpc32/64/80, aunque se ignora en Cygwin. Pero tenga en cuenta que modificará el tamaño de la mantisa, pero no el exponente, dejando la puerta abierta a otros tipos de comportamiento diferente.

En la arquitectura x86_64, SSE se usa como dice tmandry , por lo que el problema no ocurrirá a menos que fuerce la antigua FPU x87 para la computación FP con -mfpmath=387, o a menos que compila en modo de 32 bits con -m32 (necesitará el paquete multilib). Podría reproducir el problema en Linux con diferentes combinaciones de banderas y versiones de gcc:

g++-5 -m32 floating.cpp -O1
g++-8 -mfpmath=387 floating.cpp -O1

Probé algunas combinaciones en Windows o Cygwin con VC++/gcc/tcc pero el error nunca apareció. Supongo que la secuencia de instrucción generada no es la misma.

Finalmente, tenga en cuenta que una forma exótica de evitar este problema con 4.45 o 4.55 sería usar _Decimal32/64/128, pero el soporte es realmente escaso... Me pasó mucho tiempo solo para poder hacer un printf con libdfp !