Porte de código numérico powerpc a intel da diferentes resultados usando float

Question

Mi problema esencial es cómo hacer que la aritmética con flotadores en x86 se comporte como un PowerPC, yendo de Classic MacOS (CodeWarrior) a Windows (VS 2008).

El código en cuestión, de los cuales hay muchos, tiene una pila de algoritmos que son altamente iterativos y numéricamente muy sensibles.

una línea compleja típica es:

Ims_sd = sqrt((4.0*Ams*sqr(nz)-8.0*(Ams+Dms)*nz+12.0*sqr(Ams))/
     (4.0*sqr(Ams)*(sqr(nz)-1)) - 
     sqr(Ims_av))*sqrt(nz-1); 

Está escrito usando un typedef'd float como el tipo de base.

Cambiando a double obtiene resultados muy similares en ambas plataformas, pero desafortunadamente los números no son aceptables por lo que no podemos tomar ese camino tan fácil.

El código de Mac se compila utilizando CodeWarrior y simplemente apagando la generación de las instrucciones FMADD & FMSUB tuvo un efecto drástico en los números creados. Por lo tanto, mi punto de partida fue buscar las opciones de Visual Studio (2008) que parecían más similares, asegurándome de que el fusionado add fuera en uso. Sospechamos que la clave está en el comportamiento del compilador al asignar almacenamiento intermedio en los cómputos

Actualmente, los mejores resultados se obtienen con una combinación de habilitar SSE2 y /fp:fast. La habilitación de funciones intrínsecas hace que los valores se desvíen más de los valores de Mac.

La documentación del conmutador /fp dice que solo /fp:strict desactiva el comportamiento de agregar fusionado.

MSDN habla de vincular FP10.OBJ "antes de LIBC.LIB, LIBCMT.LIB o MSVCRT.LIB". a garantizan una precisión de 64 bits. Al parecer, he logrado esto al especificar FP10.OBJ en el campo de entrada del enlazador (la salida del enlazador detallado lo muestra antes de MSVCRTD.lib).

también he ajustado precisión de 64 bits mediante la invocación de

_controlfp_s(&control_word, _PC_64, MCW_PC); 

en DllMain.

Tenga en cuenta que el problema es no debido a diferencias en el manejo de excepciones de punto flotante entre plataformas, ni es debido a la forma (una delicia) que PowerPC permite la división por cero enteros (que acaban de volver a cero) ya que estas áreas ya han sido auditadas y dirigido, gracias enormemente al PC-Lint. El programa se ejecuta y produce resultados algo plausibles, pero no lo suficientemente buenos.

ACTUALIZACIÓN:

un interesante comentario de un amigo: Una posibilidad es que el PPC tiene un gran número de registros temporales que pueden almacenar valores intermedios de 64 bits mientras que el código x86 puede tener que descargar y volver a cargar el FPU (truncando a 4 bytes y perdiendo precisión).

Esta es la razón por la que SSE2 funciona mejor ya que (IIRC) tiene más registros y más posibilidades de preservar los valores intermedios.

Una posibilidad: ¿se puede compilar su código como de 64 bits? El modo x64 también tiene más registros para intermedios y mejores instrucciones de FP, por lo que puede estar más cerca del PPC en diseño y ejecución.

En realidad, las pruebas iniciales con una compilación de 64 bits se acercaron, como sugirió que podría ser (primero pensé que se sobrepasó, pero eso se debió a una configuración de modelado incorrecta).

Resolución Final

Estoy seguro de que cualquier persona interesada en este tema es lo suficientemente obsesivo que les gustaría saber cómo se resolvió todo en el final. El software está terminado y brinda resultados numéricos consistentes. Nunca pudimos obtener todos los algoritmos para entregar resultados idénticos a la Mac, pero estaban lo suficientemente cerca como para ser estadísticamente aceptable. Dado que el procesamiento es guiado por un usuario experto que selecciona las áreas de interés y que la contribución del usuario es en parte reactiva a la evolución del modelo, el director científico lo consideró aceptable (¡esta no fue una decisión de la noche a la mañana!). Las diferencias numéricas restantes están dentro de los límites de lo que determina los diferentes resultados clínicos, por lo que no se han observado diagnósticos diferentes con las pruebas.

Answer 1

La cuestión del determinismo de coma flotante en múltiples plataformas parece ser un tema muy espinoso y cuanto más se profundiza en ello, peor parece ser.

Encontré this interesting article que analiza el problema en gran profundidad, podría arrojar algunas ideas.

Answer 2

No es una respuesta como tal, sino más texto (y formato) de lo que podría caber en un comentario. Al leer su pregunta, me parece que probablemente haya considerado todo esto, pero no nos lo haya dicho, por lo que puede tratarse de una charla irrelevante. Si es así, me disculpo.

¿Puede usted (¿lo hizo?) Hacer cumplir las reglas IEEE754 para la aritmética de coma flotante en las versiones originales o portadas del programa? Mi primera suposición es que las dos plataformas (combinación de hardware, o/s, bibliotecas) implementan diferentes enfoques para la aritmética fp.

Qué suposiciones (si las hay) sobre los tamaños predeterminados, en las dos plataformas, de algunos de los tipos fundamentales como ints y flotantes. El estándar C (y creo que el estándar C++) permite la dependencia de la plataforma para algunos de ellos (no puedo olvidarlo, soy un programador de Fortran).

Conjetura final - Me he acostumbrado (en mi mundo de Fortranny) a especificar constantes flotantes como su 4.0 con suficientes dígitos para especificar todos los dígitos (decimales) en la representación preferida, es decir, algo así como 4.000000000000000000000000. Sé que, en Fortran, una constante de flotación de 4 bytes, como 3.14159625, cuando se envía automáticamente a 8 bytes, no llena los bytes adicionales con los dígitos adicionales en la expresión decimal de pi. Esto puede estarle afectando.

Nada de esto realmente le ayuda a asegurarse de que la versión portada de su código produzca los mismos resultados, al igual que la versión original, solo identifique las fuentes de diferencia.

Por último, ¿su requisito es que la nueva versión produzca los mismos resultados que la versión anterior o que garantiza a sus clientes que la nueva versión produce respuestas precisas? Su pregunta deja abierta la posibilidad de que la versión anterior del programa fuera 'defectuosa' que la nueva, teniendo en cuenta todas las fuentes de error en los cómputos numéricos.

Answer 3

remito a GCC bug 323:

me gustaría dar la bienvenida a los nuevos miembros de la comunidad de errores 323, donde todas x87 errores de punto flotante en gcc vienen a morir!Todos los errores de coma flotante que utilizan el x87 son bienvenidos, a pesar de que muchos de ellos son fáciles de reparar y muchos no. ¡Todos somos una familia feliz, cometiendo el error atroz de querer precisar la FPU de propósito general más precisa del mercado!

El breve resumen es que es increíblemente tedioso llegar "verdadera" IEEE individuales de punto flotante/duplica en un x87 y sin pérdida de rendimiento significativa; usted sufre doble redondeo de denorm incluso si usa fldcw debido al rango de exponente reducido (IIRC, IEEE FP específicamente permite a las implementaciones hacer sus propias denuncias WRT). Es de suponer que se podría hacer algo como esto:

1. Ronda a infinito positivo, realizar la parte (ldresult1 conseguir) la operación, redondeo al más cercano, incluso, convertir a flotar (conseguir fresult1).
2. RTNI, realice el op, RTNE, convierta para flotar.
3. Si son iguales, excelente: tiene el resultado correcto de flotación RTNE. Si no, entonces (creo) fresult2 < fresult1, y además, fresult1 = nextafterf (fresult2, + inf), y hay dos posibilidades:
   • ldresult1 == ((doble largo) fresult1 + fresult2)/2. La respuesta "correcta" es fresult2.
   • ldresult2 == ((doble largo) fresult1 + fresult2)/2. La respuesta "correcta" es fresult1.

probablemente estoy equivocado en los detalles en alguna parte, pero esto es presumiblemente el dolor que tiene que pasar cuando se obtiene una denorm.

Y luego aparece el otro problema: estoy bastante seguro de que no hay garantía de que sqrt() devuelva el mismo resol en diferentes implementaciones (y muy seguro para las funciones trigonométricas); la única garantía que he visto es que el resultado es "dentro de 1 ulp" (presumiblemente del resultado redondeado correctamente). Depende en gran medida del algoritmo utilizado, y las CPU modernas tienen instrucciones para estos, por lo que sufrirá una importante penalización de rendimiento si intenta implementarlo en el software. Sin embargo, ISTR es una biblioteca de punto flotante "portátil" en algún lugar que se supone que debe lograr consistencia, pero no recuerdo el nombre OTTOMH.