Как реализована функция fma()

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

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

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

Есть ли что-то, чтобы убедиться, что используется настоящая FMA, когда на точность положена?

Источник

3 ответа

Решение

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

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

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

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

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

Есть ли что-то, чтобы убедиться, что используется настоящая FMA, когда на точность положена?

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

Источник

Предложение FMA Z-бозона, основанное на алгоритме Деккера, к сожалению, неверно. В отличие от двух продуктов Dekker, в более общем случае FMA величина c не известна относительно терминов продукта, и, следовательно, могут произойти неправильные отмены.

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

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

Источник

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

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 функция в статье Расширенные прецизионные числа с плавающей точкой для вычислений на GPU и из mul_sub_x функция в векторной библиотеке классов Агнера Фога. Он использует другую функцию для разделения векторов поплавков, которая разделяется по-разному. Я попытался воспроизвести скалярную версию здесь

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 кроме как на машине с аппаратным FMA (в этом случае я использую _mm_fmsub_ss тем не мение).

Источник
Другие вопросы по тегам