Sobre optimizaciones y precisión numérica

Recibí un correo diciendo que el siguiente código en PSeInt:
    a<-rc(13);
    Escribir a<a;
evidenciaba un error de lógica, ya que el resultado daba Verdadero, cuando debiera dar Falso (un número no puede ser menor que sí mismo). Pruebo el código en GNU/Linux y da Falso, pero el usuario dijo que usaba Windows. Entonces voy a Windows y efectivamente da Verdadero. Así que abro ZinjaI y me dispongo a Depurar en Windows, pero he aquí la sorpresa: en el depurador da Falso. Cuando se trata de variables sin inicializar, o manejo de direcciones de memoria, no es raro que una versión Debug dé mejor que una Release, pero en este caso estaba seguro que no se trataba de eso. Estaba a punto de ponerme a analizar el código objeto generado en cada caso (o sea, leer ensamblador, o tal vez saltar por la ventana, estaba decidiendo), pero antes noté algo interesante.

Aislé el problema y lo reproduje fuera de PSeInt. Observen este código:
    inline double StrToDbl(const string &s) { ... } 
    int main(int argc, char *argv[]) { 
        string s1="3.6055512754639891";
        string s2="3.6055512754639891";
        double d1=StrToDbl(s1);
        double d2=StrToDbl(s2);
//        cerr<<d1<<endl;
        cout<<(d1<d2)<<endl;
        return 0; 
    }

Si lo ejecuto compilado sin optimizaciones muestra "0" (o sea, lo que corresponde). Si lo compilo con el primer nivel de optimizaciones (-O1 en gcc) muestra "1". Y por ahora solo ocurre con el número que puse (raiz de 13, probé con algunos otros y no pasó). Conclusión: WTF!

Probé reemplazar -O1 por la lista de opciones de optimización individuales que se supone que están en ese grupo (un montón de -fxxx según el manual de gcc) para aislar la conflictiva, pero el problema así no se manifestaba. Hasta que vi un detalle fundamental: si descomento la linea del cerr el problema no aparece. Y eso da la pauta de por donde pasa el problema. Una optimización muy conocida y fácil de entender es la siguiente: si una variable se usa en una parte muy chica del código, y no ocupa mucho, a veces se puede colocar en un registro del micro y nunca pasar a la memoria ram. Esto lo vemos al depurar código optimizado cuando queremos inspeccionar una variable y el resultado es "<value optimized out>". Con el cerr comentado, eso es lo que pasa con uno de los doubles. Con el cerr descomentado, las variables se usan para más cosas y más lejos (el cerr es bastante grande en ensamblador), y entonces así no aplica esta optimización.

Pero la siguiente pregunta debería ser: ¿no se supone que un programa optimizado tiene que dar exactamente la misma salida? O sea, se supone que optimizar es hacer que ande más rápido, no cambiar el resultado, entonces el compilador solo debería hacer cambios que no alteren el resultado. Y además, ¿cómo es que esto altera el resultado?

La historia es así: los flotantes siempre traen dolores de cabeza. Sabemos que la PC no puede representar los flotantes con precisión. Por ejemplo, no se puede guardar un simple 0.4 en un float. El 0.4 escrito en binario es periódico, al guardarlo en el float se trunca y queda 0.[muchos-nueves]. Al mostrarlo en pantalla se produce cierto redondeo y entonces parece un 0.4, pero no, no lo es. Y si al 0.4 lo guardo en un double, tiene muchos más nueves que si lo guardo en un float. Y para terminar de complicarla, resulta que la parte del micro que hace las operaciones de punto flotante, en muchas arquitecturas trabaja en realidad con un tamaño que no es ni el del float (32 bits), ni el del double (64 bits), sino uno mayor (80 bits), que después se "trunca" para encajar en un float o double (en realidad redondea). Entonces internamente el micro los almacena y procesa con mayor precisión. El punto es que si queda en un registro del micro, y no pasa a la memoria, nunca sufre ese truncamiento. Por eso, en el código del ejemplo d1 es diferente de d2, porque d2 fue optimizado y quedó en un registro (sin truncar), mientras que d1 no.

Se soluciona diciéndole a gcc que no deje los flotantes en registros (con -ffloat-store). De hecho, el manual de gcc, en la descripción de este argumento aclara que hay que usarlo para que el resultado se atenga a lo que dice la IEEE sobre cómo deben ser los flotantes, pero nadie lee por gusto las mil y un opciones de gcc en el manual, así que uno no lo sabe hasta que no le toca (y aún así hay que encontrar la que toca).

Y para cerrar, comento que esto reabre una vieja duda mía. ¿Debe PSeInt hacerse cargo de los errores numéricos en las comparaciones? Se sabe por ejemplo que nunca hay que comparar flotantes por igual, sino ver si la diferencia esta dentro del epsilon de máquina o algún truco similar (porque eso tampoco es 100% correcto, pero no es el punto), para evitar problemas por el error de redondeo en los cálculos. Pero eso solo lo sabe alguien a quien le importa esa diferencia (alguien que programa cosas relacionadas al cálculo numérico por ejemplo), y tal vez no debería molestar al alumno que está tratando de aprender a usar variables y estructuras de control. Entonces, tal vez PSeInt deba introducir ese epsilon en cada comparación sin que el alumno se de cuenta, para ahorrarle un problema más. Desde hace un tiempo lo incluye solo en las comparaciones por igual, que es donde más se manifiesta, pero para ser coherente debería incluirlo por todos lados, no incluirlo nunca, o trabajar con precisión arbitraria (con bibliotecas como la gmp), pero esto último haría al intérprete más lento y no se notaría la mejora en el 99% de los casos. Como siempre, se aceptan sugerencias.

Actualización(04/07/2013): una alternativa para obligar al compilador a leer y escribir una variable en memoria y no dejarla en un registro es declararla con la palabra clave "volatile". Dado que -ffloat-store puede afectar notablemente la performance si el programa hace cálculos de punto flotante intensivamente (en un ejemplo que medí de generación de mayas las tiempos varían entre 25% y 35% aproximadamente según las demás opciones de compilación y el s.o.), puede ser conveniente usar este segundo mecanismo para hacer que solo ciertas variables sean pasadas obligatoriamente a memoria