Как реализован fma()

Question

Согласно documentation, в math.h есть функция fma(). Это очень приятно, и я знаю, как работает FMA и для чего его использовать. Однако я не уверен, как это реализовано на практике? Меня больше всего интересуют архитектуры x86 и x86_64.

Есть ли инструкция с плавающей точкой (не вектор) для FMA, возможно, как определено в IEEE-754 2008?

Используется ли инструкция FMA3 или FMA4?

Есть ли внутренняя цель, чтобы убедиться, что используется реальная FMA, когда на точность полагаются?

Answer 1

Фактическая реализация зависит от платформы, но говорить очень широко:

• Если вы сообщите ваш компилятор для целевой машины с инструкциями аппаратного FMA (PowerPC, ARM с VFPv4 или AArch64, Intel Haswell или AMD Bulldozer и далее), компилятор может заменить вызовы на fma(), просто отбросив соответствующую инструкцию в ваш код. Это не гарантируется, но в целом это хорошая практика. В противном случае вы получите вызов в математическую библиотеку и:

• При работе на процессоре с аппаратным FMA эти инструкции должны использоваться для реализации функции. Однако, если у вас установлена ​​более старая версия вашей операционной системы или более старая версия математической библиотеки, она может не воспользоваться этими инструкциями.

• Если вы работаете на процессоре, у которого нет аппаратного FMA, или вы используете более старую (или просто не очень хорошую) математическую библиотеку, тогда вместо этого будет использоваться программная реализация FMA. Это может быть реализовано с использованием умных трюков с плавающей запятой с расширенной точностью или с целочисленной арифметикой.

• Результат функции fma() всегда должен быть правильно закруглен (т. Е. «Реальная fma»). Если это не так, это ошибка в математической библиотеке вашей системы. К сожалению, fma() является одной из наиболее сложных функций математической библиотеки для правильной реализации, поэтому многие реализации имеют ошибки. Пожалуйста, сообщите об этом поставщику вашей библиотеки, чтобы они исправлены!

Есть ли внутренняя, чтобы убедиться, что реальный FMA используется, когда точность полагаться?

Учитывая хороший компилятор, это не обязательно; достаточно использовать функцию fma() и сообщить компилятору, какую архитектуру вы планируете использовать. Однако компиляторы не идеальны, поэтому вам может понадобиться использовать _mm_fmadd_sd() и связанные с ними функции на x86 (но сообщать об ошибке вашему поставщику компилятора!)

Answer 2

Один из способов реализации FMA в программном обеспечении - это разделение значительных на высокие и бит.Я использую Dekker's algorithm

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}; 
} 

После того, как вы разделяете поплавок можно вычислить a*b-c с одним округлением, как этот

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; 
} 

В основном это вычитает c от (ahi,alo)*(bhi,blo) = (ahi*bhi + ahi*blo + alo*bhi + alo*blo).

Я получил эту идею от функции twoProd в статье Extended-Precision Floating-Point Numbers for GPU Computation и от функции mul_sub_x в Agner Fog's vector class library. Он использует другую функцию для расщепления векторов поплавков, которые расщепляются по-разному. Я попытался воспроизвести скалярную версию здесь

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}; 
} 

В любом случае с помощью split или split2 в fmsub хорошо согласуется с fma(a,b,-c) из математической библиотеки в Glibc. По какой-то причине моя версия значительно быстрее, чем fma, за исключением машины, имеющей аппаратное обеспечение (в этом случае я использую _mm_fmsub_ss).

Answer 3

Условное обозначение FMA Boson по алгоритму Деккера, к сожалению, неверно. В отличие от двухпроцессов Деккера, в более общем случае FMA величина c неизвестна относительно термов продукта, и, следовательно, могут иметь место неправильные отмены.

Итак, хотя два продукта Dekker могут быть значительно ускорены с помощью аппаратного FMA, вычисление суммы ошибок для двух продуктов Dekker составляет , а не - надежная реализация FMA.

Для правильной реализации потребуется либо использовать алгоритм суммирования с более высокой, чем двойной точностью, либо добавить термины в порядке убывания.