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

109

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 ++?

Oso
fuente
51
Para todos los nuevos usuarios de SO: así es como se hace una pregunta. +1
diez cuatro
1
FWIW, obtengo la salida correcta con g ++ 4.5.0 usando MinGW.
Steve Blackwell
2
ideone usa 4.3.4 ideone.com/b8VXg
Daniel A. White
5
Debe tener en cuenta que es poco probable que su rutina funcione de manera confiable con todo tipo de resultados. A diferencia de redondear un doble a un número entero, esto es vulnerable al hecho de que no todos los números reales se pueden representar, por lo que debe esperar obtener más errores como este.
Jakub Wieczorek
2
Para aquellos que no pueden reproducir el error: no eliminen los comentarios de los stmts de depuración comentados, afectan el resultado.
n. 'pronombres' m.

Respuestas:

91

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.

Maxim Egorushkin
fuente
20
Creo que esta es la respuesta. La constante 4.55 se convierte en 4.54999999999999 que es la representación binaria más cercana en 64 bits; multiplique por 10 y redondee de nuevo a 64 bits y obtendrá 45,5. Si omite el paso de redondeo manteniéndolo en un registro de 80 bits, obtendrá 45.4999999999999.
Mark Ransom
Gracias, ni siquiera conozco esta opción. Pero me preguntaba, ¿debería activar siempre la opción -float-store? Debido a que la versión g ++ que probé se envía con CentOS / Redhat 5 y CentOS / Redhat 6. He compilado muchos de mis programas en estas plataformas, me preocupa que eso cause errores inesperados dentro de mis programas.
Bear
5
@ Bear, la declaración de depuración probablemente hace que el valor se vacíe de un registro a la memoria.
Mark Ransom
2
@Bear, normalmente su aplicación debería beneficiarse de una precisión extendida, a menos que opere en valores extremadamente pequeños o grandes cuando se espera que un flotante de 64 bits se exceda o se desborde y produzca inf. No existe una buena regla general, las pruebas unitarias pueden darte una respuesta definitiva.
Maxim Egorushkin
2
@bear Como regla general, si necesita resultados que sean perfectamente predecibles y / o exactamente lo que un humano obtendría haciendo las sumas en papel, entonces debe evitar el punto flotante. -float-store elimina una fuente de imprevisibilidad, pero no es una fórmula mágica.
Plugwash
10

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.

David Hammen
fuente
6

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.

Perrito
fuente
4
Hay un -ffast-mathinterruptor para GCC que, y no está activado por ninguno de los -Oniveles de optimización desde la cita: "puede dar como resultado una salida incorrecta para programas que dependen de una implementación exacta de las reglas / especificaciones IEEE o ISO para funciones matemáticas".
Mat
@Mat: Lo intenté -ffast-mathy algunas otras cosas en mi g++ 4.4.3y todavía no puedo reproducir el problema.
NPE
Agradable: con -ffast-mathobtengo 4.5en ambos casos para niveles de optimización superiores a 0.
Kerrek SB
(Corrección: obtengo 4.5con -O1y -O2, pero no con -O0y -O3en GCC 4.4.3, pero con -O1,2,3en GCC 4.6.1.)
Kerrek SB
4

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().

Max Lybbert
fuente
Un "épsilon relativo" se llama 1 ulp (1 unidad en el último lugar). x - nextafter(x, INFINITY)está relacionado con 1 ulp para x (pero no lo use; estoy seguro de que hay casos de esquina y lo acabo de inventar). El ejemplo de referencia de cpp para epsilon() tiene un ejemplo de escalado para obtener un error relativo basado en ULP .
Peter Cordes
2
Por cierto, la respuesta de 2016 a -ffloat-storees: no use x87 en primer lugar. Utilice matemáticas SSE2 (binarios de 64 bits, o -mfpmath=sse -msse2para hacer binarios viejos y crujientes de 32 bits), porque SSE / SSE2 tiene provisionales sin precisión adicional. doubley las floatvariables en los registros XMM están realmente en formato IEEE de 64 o 32 bits. (A diferencia de x87, donde los registros son siempre de 80 bits y el almacenamiento en memoria se redondea a 32 o 64 bits)
Peter Cordes
3

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.

tmandry
fuente
1

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!

calandoa
fuente
0

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.

cdcdcd
fuente