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

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

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

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

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

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

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

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

Одним из способов программной реализации 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 в библиотеке векторных классов Agner Fog. Он использует другую функцию для разделения векторов поплавков, которая разделяется по-разному. Я пытался воспроизвести здесь скалярную версию

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).

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

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

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