¿Está rota mi fma()?

10 minutos de lectura

avatar de usuario
chux – Reincorporar a Monica

en el uso double fma(double x, double y, double z); Esperaría un distinto de cero d en las líneas de salida a continuación marcadas con '?'. Eso aparece para uso interno solamente long double precisión en lugar de precisión infinita Como se especificó.

los fma funciones calculan (x × y) + zredondeado como una operación ternaria: calculan el valor (como si) para precisión infinita y redondee una vez al formato de resultado, según el modo de redondeo actual. §7.12.13.1 2 (énfasis mío)

también lo es mi fma() roto, o ¿cómo lo estoy usando incorrectamente en el código o en las opciones de compilación?

#include <float.h>
#include <math.h>
#include <stdio.h>

int main(void) {
  // Invoking: Cygwin C Compiler
  // gcc -std=c11 -O0 -g3 -pedantic -Wall -Wextra -Wconversion -c -fmessage-length=0 
  //   -v -MMD -MP -MF"x.d" -MT"x.o" -o "x.o" "../x.c"

  printf("FLT_EVAL_METHOD %d\n", FLT_EVAL_METHOD);
  for (unsigned i = 20; i < 55; i++) {
    volatile double a = 1.0 + 1.0 / pow(2, i);
    volatile double b = a;
    volatile double c = a * b;
    volatile double d = fma(a, b, -c);
    volatile char *nz = ((i >= 27 && a != 1.0) == !d) ? "?" : "";
    printf("i:%2u a:%21.13a c:%21.13a d:%10a %s\n", i, a, c, d, nz);
  }
  return 0;
}

Producción

FLT_EVAL_METHOD 2
i:20 a: 0x1.0000100000000p+0 c: 0x1.0000200001000p+0 d:    0x0p+0 
i:21 a: 0x1.0000080000000p+0 c: 0x1.0000100000400p+0 d:    0x0p+0 
i:22 a: 0x1.0000040000000p+0 c: 0x1.0000080000100p+0 d:    0x0p+0 
i:23 a: 0x1.0000020000000p+0 c: 0x1.0000040000040p+0 d:    0x0p+0 
i:24 a: 0x1.0000010000000p+0 c: 0x1.0000020000010p+0 d:    0x0p+0 
i:25 a: 0x1.0000008000000p+0 c: 0x1.0000010000004p+0 d:    0x0p+0 
i:26 a: 0x1.0000004000000p+0 c: 0x1.0000008000001p+0 d:    0x0p+0 
i:27 a: 0x1.0000002000000p+0 c: 0x1.0000004000000p+0 d:   0x1p-54 
i:28 a: 0x1.0000001000000p+0 c: 0x1.0000002000000p+0 d:   0x1p-56 
i:29 a: 0x1.0000000800000p+0 c: 0x1.0000001000000p+0 d:   0x1p-58 
i:30 a: 0x1.0000000400000p+0 c: 0x1.0000000800000p+0 d:   0x1p-60 
i:31 a: 0x1.0000000200000p+0 c: 0x1.0000000400000p+0 d:   0x1p-62 
i:32 a: 0x1.0000000100000p+0 c: 0x1.0000000200000p+0 d:    0x0p+0 ?
i:33 a: 0x1.0000000080000p+0 c: 0x1.0000000100000p+0 d:    0x0p+0 ?
i:34 a: 0x1.0000000040000p+0 c: 0x1.0000000080000p+0 d:    0x0p+0 ?
...
i:51 a: 0x1.0000000000002p+0 c: 0x1.0000000000004p+0 d:    0x0p+0 ?
i:52 a: 0x1.0000000000001p+0 c: 0x1.0000000000002p+0 d:    0x0p+0 ?
i:53 a: 0x1.0000000000000p+0 c: 0x1.0000000000000p+0 d:    0x0p+0 
i:54 a: 0x1.0000000000000p+0 c: 0x1.0000000000000p+0 d:    0x0p+0 

Información de la versión

gcc -v

Using built-in specs.
COLLECT_GCC=gcc
COLLECT_LTO_WRAPPER=/usr/lib/gcc/i686-pc-cygwin/5.3.0/lto-wrapper.exe
Target: i686-pc-cygwin
Configured with: /cygdrive/i/szsz/tmpp/gcc/gcc-5.3.0-5.i686/src/gcc-5.3.0/configure --srcdir=/cygdrive/i/szsz/tmpp/gcc/gcc-5.3.0-5.i686/src/gcc-5.3.0 --prefix=/usr --exec-prefix=/usr --localstatedir=/var --sysconfdir=/etc --docdir=/usr/share/doc/gcc --htmldir=/usr/share/doc/gcc/html -C --build=i686-pc-cygwin --host=i686-pc-cygwin --target=i686-pc-cygwin --without-libiconv-prefix --without-libintl-prefix --libexecdir=/usr/lib --enable-shared --enable-shared-libgcc --enable-static --enable-version-specific-runtime-libs --enable-bootstrap --enable-__cxa_atexit --with-dwarf2 --with-arch=i686 --with-tune=generic --disable-sjlj-exceptions --enable-languages=ada,c,c++,fortran,java,lto,objc,obj-c++ --enable-graphite --enable-threads=posix --enable-libatomic --enable-libcilkrts --enable-libgomp --enable-libitm --enable-libquadmath --enable-libquadmath-support --enable-libssp --enable-libada --enable-libjava --enable-libgcj-sublibs --disable-java-awt --disable-symvers --with-ecj-jar=/usr/share/java/ecj.jar --with-gnu-ld --with-gnu-as --with-cloog-include=/usr/include/cloog-isl --without-libiconv-prefix --without-libintl-prefix --with-system-zlib --enable-linker-build-id --with-default-libstdcxx-abi=gcc4-compatible
Thread model: posix
gcc version 5.3.0 (GCC) 

  • Si significa algo, esto aparece arcoíris y mariposas usando el sonido clang 3.8 dist en mi mac. No ? para ser encontrado.

    – WhozCraig

    10 de febrero de 2017 a las 18:51

  • @WhozCraig Parece que estoy experimentando el punto n. ° 4 de esta respuesta relacionada a diferencia de su buena plataforma.

    – chux – Reincorporar a Monica

    10 de febrero de 2017 a las 18:56


  • Bien que apesta. (como si tuviera que decírtelo). Sry, hombre.

    – WhozCraig

    10 de febrero de 2017 a las 18:58

  • @chux Usando gcc 6.3.1 ¡Todo se ve bien! no ? ser visto.

    – PeCosta

    10 de febrero de 2017 a las 19:54

  • Parece que es un problema con su biblioteca de matemáticas: ¿qué significa ldd en el ejecutable resultante dar?

    – Simón Byrne

    10 de febrero de 2017 a las 21:51

avatar de usuario
Animal Nominal

Es culpa de Cygwin. O más correctamente, el nueva biblioteca biblioteca C que utiliza. Eso dice explícitamente ni siquiera trata de conseguir fma() emulación correcta.

La biblioteca GNU C tiene una emulación correcta para casi todas las variantes de fma desde 2015. Para obtener detalles y los parches utilizados para implementar esto, consulte el error de código fuente 13304.

Si la eficiencia no es un problema, simplemente usaría, por ejemplo

#if defined(__CYGWIN__) && !defined(__FMA__) && !defined(__FMA3__) && !defined(__FMA4__)
#define fma(x, y, z)  fma_emulation(x, y, z)

double fma_emulation(double x, double y, double z)
{
    /* One of the implementations linked above */
}
#endif

Personalmente, no uso Windows en absoluto, pero si alguien lo hace (usa Windows y necesita la emulación fma), le sugiero que intente enviar un parche aguas arriba, con un enlace a la Discusión de la biblioteca GNU C sobre la emulación correcta de fma.


Lo que me pregunto es si sería posible examinar sólo la baja METRO bits del resultado (descartados en el redondeo) para determinar el valor correcto de la ULP en el resultado, y ajustar el resultado obtenido usando el sencillo a×B+C operación en consecuencia, utilizando nextafter(); en lugar de usar aritmética de multiprecisión para implementar toda la operación.

Editar: No, porque la adición puede desbordarse, eliminando un bit adicional como el MSB de la parte descartada. Solo por esa razón, necesitamos hacer toda la operación. Otra razón es que si a×B y C tienen diferentes signos, luego en lugar de sumar, restamos la magnitud más pequeña de la magnitud más grande (resultado usando el signo del más grande), lo que puede conducir a borrar varios bits altos en el más grande, y eso afecta qué bits del resultado completo se eliminan en el redondeo

Sin embargo, para IEEE-754 Binary64 double en arquitecturas x86 y x86-64, creo que la emulación fma usando registros de 64 bits (enteros) y productos de 128 bits sigue siendo bastante factible. Experimentaré con una representación en la que se utilizan 2 bits bajos en un registro de 64 bits para los bits de decisión de redondeo (siendo LSB un OR lógico de todos los bits perdidos), 53 bits utilizados para la mantisa y un bit de acarreo, dejando 8 bits altos no utilizados e ignorados. El redondeo se realiza cuando la mantisa entera sin signo se convierte en un doble (64 bits). Si estos experimentos arrojan algo útil, los describiré aquí.


Hallazgos iniciales: fma() la emulación en un sistema de 32 bits es lenta. Las cosas de 80 bits en la FPU 387 son básicamente inútiles aquí, e implementar la multiplicación de 53 × 53 bits (y el cambio de bits) en un sistema de 32 bits simplemente no vale la pena. la glibc fma() El código de emulación vinculado anteriormente es lo suficientemente bueno en mi opinión.

Hallazgos adicionales: el manejo de valores no finitos es desagradable. (Los subnormales son solo un poco molestos y requieren un manejo especial (ya que el MSB implícito en la mantisa es cero entonces).) Si alguno de los tres argumentos no es finito (infinito o alguna forma de NaN), entonces regresa a*b + c (no fusionado) es la única opción sana. El manejo de estos casos requiere bifurcaciones adicionales, lo que ralentiza la emulación.

Decisión final: la cantidad de casos que se manejarán de manera optimizada (en lugar de usar el enfoque de “limb” de precisión múltiple como se usa en la emulación glibc) es lo suficientemente grande como para que este enfoque no valga la pena. Si cada miembro es de 64 bits, cada uno de a, By C se distribuye en un máximo de 2 extremidades, y a×B sobre tres extremidades. (Con miembros de 32 bits, son solo 3 y 5 miembros, respectivamente). Dependiendo de si a×B y C tienen signos iguales o diferentes, solo hay dos casos fundamentalmente diferentes para manejar: en el caso de signos diferentes, la suma se convierte en resta (más pequeño de más grande, el resultado es el mismo signo que el valor más grande).

En resumen, el enfoque de multiprecisión es mejor. La precisión real necesaria está muy bien delimitada y ni siquiera necesita una asignación dinámica. Si el producto de las mantisas de a y B se puede calcular de manera eficiente, la parte de multiprecisión se limita a sostener el producto y manejar la suma/resta. El redondeo final se puede realizar convirtiendo el resultado a una mantisa de 53 bits, un exponente y dos bits extra bajos (siendo el mayor el bit más significativo perdido en el redondeo y el menor un OR del resto de los bits perdidos en el redondeo). el redondeo). Esencialmente, las operaciones clave se pueden realizar usando números enteros (o registros SSE/AVX), y la conversión final de una mantisa de 55 bits a una doble maneja el redondeo de acuerdo con las reglas actuales.

  • “es posible examinar solo los M bits bajos del resultado (descartados en el redondeo)” –> Yo no diría eso. Los bits descartados incluyen el inferior M pedacitos de la 2M producto y el M pedazos de c que puede estar muy a la “derecha” y contribuir a la ronda. Tal como lo veo, la suma precisa puede tener un ancho de bit de hasta (2*M + |Ae + Be – 2*Ce|) (o algo que depende en gran medida de la diferencia de exponente de AB frente a C. (Supongamos que el bit[0] es MSBit) El bit de la suma[M-1]poco[M] y el “o” de todos los bits menos significativos contribuyen a la ronda.

    – chux – Reincorporar a Monica

    11 de febrero de 2017 a las 17:24


  • Me siento como Snoopy sacudiendo su puño en el cygwin fma/Red-Barron

    – chux – Reincorporar a Monica

    11 de febrero de 2017 a las 17:28

  • @chux: No, no se deje engañar por el análogo de “precisión infinita”. Sólo tenemos que considerar los bits que afectar el redondeo. Hay seis casos básicos: ab y c tienen los mismos exponentes y los mismos signos; diferentes signos; ab tiene mayor exponente, pero mismo signo que c; diferentes signos; y ab tiene un exponente menor que c, pero mismo signo; y diferentes signos. Para los cuatro modos de redondeo de IEEE-754, solo necesitamos conocer el bit más significativo de la parte que se pierde en el redondeo, y podemos hacerlo en cada uno de los seis casos con registros de 2M-bit. ¿Debería dar más detalles sobre esto? (¿Es útil?)

    – Animal nominal

    11 de febrero de 2017 a las 18:08


  • Para redondeo IEEE predeterminado lazos para incluso se necesita el MSBit de la parte soltada (estamos de acuerdo en eso) y también si cualquier otro bit perdido es 1. (esto parece que no estamos de acuerdo). Si cualquiera de los otros bits descartados es 1, entonces el valor es encima a la mitad y se redondea lejos de 0. Si todos los demás bits eliminados son 0, el valor se redondea a par, ya que es un Corbata. Todos los bits de caída pueden afectar en ese modo de redondeo.

    – chux – Reincorporar a Monica

    11 de febrero de 2017 a las 19:29


  • @chux: Entonces, mi idea no funcionó. Sin embargo, su nota sobre el efecto de los bits eliminados en el redondeo me recordó que si empaquetamos los 53 bits más significativos de la mantisa, el bit más significativo de los bits eliminados y un bit adicional que es el OR lógico del resto de los bits descartados, obtenemos un entero sin signo de 55 bits que se redondea correctamente cuando se convierte en un doble. Esto significa que un enfoque de precisión múltiple que utiliza un número pequeño (límite superior fijo) de extremidades debería manejar la emulación de manera eficiente. El enfoque de redondeo puede que ser novedoso, pero lo dudo.

    – Animal nominal

    12 de febrero de 2017 a las 11:42


¿Ha sido útil esta solución?

Esta web utiliza cookies propias y de terceros para su correcto funcionamiento y para fines analíticos y para mostrarte publicidad relacionada con sus preferencias en base a un perfil elaborado a partir de tus hábitos de navegación. Al hacer clic en el botón Aceptar, acepta el uso de estas tecnologías y el procesamiento de tus datos para estos propósitos. Configurar y más información
Privacidad