Как се изпълнява 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};
}

След като разделите float, можете да изчислите 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, за съжаление е неправилно. За разлика от twoProduct на Dekker, в по-общия случай на FMA големината на c не е известна спрямо условията на продукта и следователно могат да възникнат грешни анулации.

И така, докато twoProduct на Dekker може да бъде значително ускорен с хардуерен FMA, изчислението на термина за грешка на twoProduct на Dekker не е стабилна FMA реализация.

Правилното внедряване ще трябва или да използва алгоритъм за сумиране с по-висока от двойна точност, или да добави термините в низходящ ред на величина.