¿Está rota mi fma()?

Question 1

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)

Question 2

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.

¿Ha sido útil esta solución?