floating point - ¿Cómo se implementa fma()?

floating-point ieee-754

Question

De acuerdo con la documentación , hay una función fma() en math.h Eso es muy bueno, y sé cómo funciona FMA y para qué usarlo. Sin embargo, no estoy tan seguro de cómo esto se implementa en la práctica? Estoy muy interesado en las x86 y x86_64 .

¿Hay alguna instrucción de punto flotante (no vectorial) para FMA, tal como se define en IEEE-754 2008?

¿Se utilizan las instrucciones FMA3 o FMA4?

¿Existe un elemento intrínseco para garantizar que se utilice un FMA real cuando se confía en la precisión?

Answer 1

La implementación real varía de una plataforma a otra, pero hablando muy ampliamente:

• Si le dice a su compilador que apunte a una máquina con instrucciones FMA de hardware (PowerPC, ARM con VFPv4 o AArch64, Intel Haswell o AMD Bulldozer en adelante), el compilador puede reemplazar las llamadas a fma( ) simplemente dejando caer la instrucción apropiada en su código. Esto no está garantizado, pero generalmente es una buena práctica. De lo contrario, recibirá una llamada a la biblioteca de matemáticas y:

• Cuando se ejecuta en un procesador que tiene hardware FMA, esas instrucciones se deben usar para implementar la función. Sin embargo, si tiene una versión anterior de su sistema operativo o una versión anterior de la biblioteca matemática, es posible que no tome ventaja de esas instrucciones.

• Si está ejecutando un procesador que no tiene hardware FMA, o está utilizando una biblioteca matemática anterior (o no muy buena), entonces se utilizará una implementación de software de FMA. Esto podría implementarse usando astutos trucos de punto flotante de precisión extendida o con aritmética de enteros.

• El resultado de la función fma( ) siempre debe redondearse correctamente (es decir, una "fma real"). Si no es así, es un error en la biblioteca matemática de su sistema. Desafortunadamente, fma( ) es una de las funciones de biblioteca matemática más difíciles de implementar correctamente, por lo que muchas implementaciones tienen errores. ¡Por favor infórmenos a su proveedor de la biblioteca para que se solucionen!

¿Existe un elemento intrínseco para garantizar que se utilice un FMA real cuando se confía en la precisión?

Dado un buen compilador, esto no debería ser necesario; debería bastar con usar la función fma( ) y decirle al compilador a qué arquitectura se dirige. Sin embargo, los compiladores no son perfectos, por lo que puede necesitar usar _mm_fmadd_sd( ) e intrínsecos relacionados en x86 (¡pero informe el error al proveedor del compilador!)

Answer 2

Una forma de implementar FMA en el software es dividir lo significativo en bits altos y bajos. Yo uso el algoritmo de Dekker

typedef struct { float hi; float lo; } doublefloat;  
doublefloat split(float a) {
    float t = ((1<<12)+1)*a;
    float hi = t - (t - a);
    float lo = a - hi;
    return (doublefloat){hi, lo};
}

Una vez que divide el flotador puede calcular a*bc con un solo redondeo como este

float fmsub(float a, float b, float c) {
    doublefloat as = split(a), bs = split(b);
    return ((as.hi*bs.hi - c) + as.hi*bs.lo + as.lo*bs.hi) + as.lo*bs.lo;
}

Esto básicamente resta c de (ahi,alo)*(bhi,blo) = (ahi*bhi + ahi*blo + alo*bhi + alo*blo) .

twoProd esta idea de la función twoProd en los Números de punto flotante de precisión extendida en papel para Computación de GPU y de la función mul_sub_x en la biblioteca de clases vectoriales de Agner Fog . Él usa una función diferente para dividir vectores de flotadores que se divide de manera diferente. Traté de reproducir una versión escalar aquí

typedef union {float f; int i;} u;
doublefloat split2(float a) {
    u lo, hi = {a};
    hi.i &= -(1<<12);
    lo.f = a - hi.f;
    return (doublefloat){hi.f,lo.f};
}

En cualquier caso, usar split o split2 en fmsub concuerda bien con fma(a,b,-c) de la biblioteca matemática en glibc. Por alguna razón, mi versión es significativamente más rápida que fma excepto en una máquina que tiene hardware fma (en cuyo caso uso _mm_fmsub_ss todos modos).

Answer 3

La sugerencia FMA de Z boson basada en el algoritmo de Dekker es desafortunadamente incorrecta. A diferencia de los dos productos de Dekker, en el caso más general de FMA, la magnitud de c no se conoce en relación con los términos del producto y, por lo tanto, pueden producirse cancelaciones incorrectas.

Entonces, mientras que los dos productos de Dekker se pueden acelerar en gran medida con un FMA de hardware, el cálculo del término de error de los dos productos de Dekker no es una implementación robusta de FMA.

Una implementación correcta necesitaría usar un algoritmo de suma con una precisión superior a la doble, o agregar los términos en orden de magnitud decreciente.