fma()是如何实现的

Question

根据文档，math.h中有一个fma()函数。这很好，我知道FMA是如何工作的，也知道它的用途。然而，我不太清楚这是如何在实践中实施的？我主要对x86和x86_64体系结构感兴趣。

FMA是否有浮点(非矢量)指令，也许是由IEEE-754 2008定义的？

是否使用FMA3或FMA4指令？

是否有一个内在的，以确保一个真正的FMA使用，当精度是依赖的？

Answer 1

各平台的实际实施情况各不相同，但从广义上讲：

• 如果您告诉编译器以硬件FMA指令(PowerPC、ARM与VFPv4或AArch64、Intel Haswell或AMD等)为目标，编译器可以通过将适当的指令丢弃到代码中来替换对fma( )的调用。这不是保证，但通常是良好做法。否则，您将接到对数学库的调用，并且：
• 在具有硬件FMA的处理器上运行时，应该使用这些指令来实现该功能。但是，如果您有旧版本的操作系统，或者数学库的旧版本，它可能不会利用这些说明。
• 如果您正在没有硬件FMA的处理器上运行，或者您使用的是一个较旧的(或者说不是很好的)数学库，那么将使用FMA的软件实现。这可以通过巧妙的扩展精度浮点技巧来实现，或者用整数算法实现.
• fma( )函数的结果应该总是正确的四舍五入(即“真正的fma")。如果不是，那就是系统数学库中的一个bug。不幸的是，fma( )是更难实现的数学库函数之一，因此许多实现都有错误。请向您的图书馆供应商报告，以便他们得到修复！

是否有一个内在的，以确保一个真正的FMA使用，当精度是依赖的？

如果有一个好的编译器，这不应该是必要的；只要使用fma( )函数并告诉编译器您的目标是什么体系结构就足够了。然而，编译器并不完美，因此您可能需要在x86上使用x86和相关的本质(但是请将错误报告给编译器供应商！)

Answer 2

在软件中实现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计算中的扩展精度浮点数和Agner Fog向量类库中的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与glibc中的数学库中的fma(a,b,-c)非常吻合。无论出于什么原因，我的版本都比fma快得多，除了在有硬件fma的机器上(在这种情况下，我还是使用_mm_fmsub_ss )。

Answer 3

不幸的是，基于Dekker算法的Z玻色子FMA建议是错误的。与Dekker的twoProduct不同，在更一般的FMA情况下，相对于乘积项，c的大小是未知的，因此可能会出现错误的取消。

因此，虽然使用硬件FMA可以大大加速德克尔的twoProduct，但是Dekker的twoProduct的错误项计算并不是一个健壮的FMA实现。

正确的实现需要使用比双精度更高的求和算法，或者以递减数量级添加项。