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

El siguiente código funciona en Visual Studio 2008 con y sin optimización. Pero solo funciona en g ++ sin optimización (O0).

#include <cstdlib>
#include <iostream>
#include <cmath>

double round(double v, double digit)
{
    double pow = std::pow(10.0, digit);
    double t = v * pow;
    //std::cout << "t:" << t << std::endl;
    double r = std::floor(t + 0.5);
    //std::cout << "r:" << r << std::endl;
    return r / pow;
}

int main(int argc, char *argv[])
{
    std::cout << round(4.45, 1) << std::endl;
    std::cout << round(4.55, 1) << std::endl;
}

La salida debe ser:

4.5
4.6

Pero g ++ con optimización ( O1- O3) dará como resultado:

4.5
4.5

Si agrego la volatilepalabra clave antes de t, funciona, entonces, ¿podría haber algún tipo de error de optimización?

Prueba en g ++ 4.1.2 y 4.4.4.

Aquí está el resultado en ideone: http://ideone.com/Rz937

Y la opción que pruebo en g ++ es simple:

g++ -O2 round.cpp

El resultado más interesante, incluso si activo la /fp:fastopción en Visual Studio 2008, el resultado sigue siendo correcto.

Otra pregunta:

Me preguntaba, ¿debería activar siempre la -ffloat-storeopción?

Porque la versión g ++ que probé se envía con CentOS / Red Hat Linux 5 y CentOS / Redhat 6 .

Compilé muchos de mis programas en estas plataformas y me preocupa que causen errores inesperados dentro de mis programas. Parece un poco difícil investigar todo mi código C ++ y bibliotecas usadas si tienen tales problemas. ¿Cualquier sugerencia?

¿Alguien está interesado en por qué incluso /fp:fastencendido, Visual Studio 2008 todavía funciona? ¿Parece que Visual Studio 2008 es más confiable en este problema que g ++?

Los procesadores Intel x86 utilizan internamente precisión extendida de 80 bits, mientras doubleque normalmente tiene un ancho de 64 bits. Los diferentes niveles de optimización afectan la frecuencia con la que los valores de punto flotante de la CPU se guardan en la memoria y, por lo tanto, se redondean de precisión de 80 bits a precisión de 64 bits.

Utilice la -ffloat-storeopción gcc para obtener los mismos resultados de punto flotante con diferentes niveles de optimización.

Alternativamente, use el long doubletipo, que normalmente tiene 80 bits de ancho en gcc para evitar el redondeo de precisión de 80 bits a 64 bits.

man gcc lo dice todo:

   -ffloat-store
       Do not store floating point variables in registers, and inhibit
       other options that might change whether a floating point value is
       taken from a register or memory.

       This option prevents undesirable excess precision on machines such
       as the 68000 where the floating registers (of the 68881) keep more
       precision than a "double" is supposed to have.  Similarly for the
       x86 architecture.  For most programs, the excess precision does
       only good, but a few programs rely on the precise definition of
       IEEE floating point.  Use -ffloat-store for such programs, after
       modifying them to store all pertinent intermediate computations
       into variables.

En las compilaciones x86_64, los compiladores usan registros SSE para floaty doublede forma predeterminada, de modo que no se usa precisión extendida y este problema no ocurre.

gccLa opción del compilador-mfpmath controla eso.

La salida debería ser: 4.5 4.6 Eso es lo que sería la salida si tuviera una precisión infinita, o si estuviera trabajando con un dispositivo que usara una representación de punto flotante basada en decimal en lugar de basada en binario. Pero no lo eres. La mayoría de las computadoras utilizan el estándar binario de coma flotante IEEE.

Como ya señaló Maxim Yegorushkin en su respuesta, parte del problema es que internamente su computadora está usando una representación de punto flotante de 80 bits. Sin embargo, esto es solo parte del problema. 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 pueda 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) siempre es inexacto. Encuentre el caso de la esquina que determina si algún valor de entrada se redondea hacia arriba o hacia abajo y devuelva el valor redondeado hacia arriba o hacia abajo según una comparación con este caso de esquina. Y debe tener cuidado de que un compilador de optimización no coloque ese caso de esquina encontrado en un registro de precisión extendido.

Consulte ¿Cómo Excel redondea correctamente los números flotantes aunque sean imprecisos?para tal algoritmo.

O simplemente puede vivir con el hecho de que las cajas de las esquinas a veces se redondean erróneamente.

Los diferentes compiladores tienen diferentes configuraciones de optimización. Algunas de esas configuraciones de optimización más rápidas no mantienen reglas estrictas de punto flotante de acuerdo con IEEE 754 . Visual Studio tiene una configuración específica, /fp:strict, /fp:precise, /fp:fast,, donde /fp:fastviola la norma en 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, lo único diferente entre los compiladores es que GCC buscaría el comportamiento de punto flotante más rápido de forma predeterminada en optimizaciones más altas, mientras que Visual Studio no cambia el comportamiento de punto flotante con niveles de optimización más altos. Por lo tanto, puede que no sea necesariamente un error real, sino el comportamiento previsto de una opción que no sabía que estaba activando.

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

Esto implica que el problema está relacionado con las declaraciones de depuración. Y parece que hay un error de redondeo causado al cargar los valores en los registros durante las declaraciones de salida, razón por la cual otros descubrieron que puede solucionarlo con-ffloat-store

Otra pregunta:

Me preguntaba, ¿debería activar siempre la -ffloat-storeopción?

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

Pero, en general, existe cierta discrepancia entre los números de coma flotante binarios (como los que usa la computadora) y los números de coma flotante decimal (con los que la gente está familiarizada), y esa discrepancia puede causar un comportamiento similar al que obtiene (para ser claros, el comportamiento que obtiene no es causado por esta falta de coincidencia, sino similar puede ser comportamiento ). La cuestión es que, dado que ya tiene cierta vaguedad al tratar con el punto flotante, no puedo decir que -ffloat-storeeso lo mejore o empeore.

En su lugar, es posible que desee buscar otras soluciones al problema que está tratando de resolver (desafortunadamente, Koenig no apunta al documento real, y realmente no puedo encontrar un lugar "canónico" obvio para él, así que Tendré que enviarte a Google ).

Si no está redondeando para fines de salida, probablemente miraría std::modf()(in cmath) y std::numeric_limits<double>::epsilon()(in limits). Pensando en la round()función original , creo que sería más limpio reemplazar la llamada astd::floor(d + .5) con una llamada a 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 sumado a 1 que crea un número distinto de 1." Por lo general, es necesario utilizar un épsilon relativo (es decir, escalar épsilon de alguna manera para tener en cuenta el hecho de que está trabajando con números distintos de "1"). La suma de d, .5y std::numeric_limits<double>::epsilon()debe estar cerca de 1, por lo que agrupar esa suma significa questd::numeric_limits<double>::epsilon() será sobre el tamaño adecuado para lo que 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 cuando no deberíamos.

Hoy en día, deberías considerarlo std::nearbyint().

La respuesta aceptada es correcta si está compilando en un destino x86 que no incluye SSE2. Todos los procesadores x86 modernos admiten SSE2, por lo que, si puede aprovecharlo, debería:

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

Analicemos esto.

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

En procesadores x86 más modernos que admiten SSE2, especificar las opciones del compilador -mfpmath=sse -msse2garantiza que todas las operaciones flotantes y dobles se realicen en registros SSE y se redondeen correctamente. Estas opciones no afectan el ABI y, por lo tanto, deben usarse 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 (fusionada multiplicar-añadir) pueden cambiar el comportamiento de redondeo en comparación con sus contrapartes no fusionadas, por lo que debemos deshabilitarlo. Este es el predeterminado en Clang, no en GCC. Como se explica en esta respuesta :

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

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

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 dijo Maxim 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, el x87 FPU se configuró para usar una precisión interna para la mantisa de 53 bits (equivalente a 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. Entonces, el problema debería ser posible, o no, en estas diferentes plataformas cambiando la palabra de control de la FPU (asumiendo que la secuencia de instrucciones desencadenaría el error). El problema se informó a gcc como error 323 (¡lea al menos el comentario 92!).

Para mostrar la precisión de la mantisa 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 ten en cuenta que modificará el tamaño de la mantisa, pero no del exponente, dejando la puerta abierta a otros tipos de comportamientos diferentes.

En la arquitectura x86_64, SSE se usa como lo dice tmandry , por lo que el problema no ocurrirá a menos que fuerce el antiguo x87 FPU para computación FP con -mfpmath=387, oa menos que compile en modo de 32 bits con -m32(necesitará un 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 prevenir este problema con 4.45 o 4.55 sería usar _Decimal32/64/128, pero el soporte es realmente escaso ... ¡Pasé mucho tiempo solo para poder hacer un printf libdfp!

Personalmente, he tenido el mismo problema yendo en sentido contrario, 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 involucran grandes conjuntos de datos de punto flotante. Incluso después de desmontar, a menudo me decepcionan las opciones de los compiladores. Muy a menudo, es más fácil usar elementos intrínsecos del compilador o simplemente escribir el ensamblado usted mismo.