fma() 是如何实现的

Question

根据documentation，math.h中有一个fma()函数。这非常好，而且我知道 FMA 的工作原理以及它的用途。但是，我不太确定这在实践中是如何实现的？我最感兴趣的是x86 和x86_64 架构。

是否有用于 FMA 的浮点（非向量）指令，可能由 IEEE-754 2008 定义？

是使用FMA3还是FMA4指令？

在依赖精度的情况下，是否存在确保使用真正的 FMA 的内在因素？

Answer 1

实际实现因平台而异，但范围很广：

• 如果您告诉编译器以具有硬件 FMA 指令（PowerPC、带有 VFPv4 或 AArch64 的 ARM、Intel Haswell 或 AMD Bulldozer 及更高版本）的机器为目标，编译器可能替换调用fma( ) 只需将适当的指令放入您的代码中即可。这不能保证，但通常是很好的做法。否则你会接到数学库的电话，并且：

• 在具有硬件 FMA 的处理器上运行时，应使用这些指令来实现该功能。但是，如果您的操作系统版本较旧，或者数学库版本较旧，则可能无法利用这些说明。

• 如果您在没有硬件 FMA 的处理器上运行，或者您使用的是较旧（或不是很好）的数学库，则将使用 FMA 的软件实现。这可以使用巧妙的扩展精度浮点技巧或整数运算来实现。

• fma( ) 函数的结果应始终正确舍入（即“真正的 fma”）。如果不是，那是系统数学库中的错误。不幸的是，fma( ) 是更难正确实现的数学库函数之一，因此许多实现都有错误。请将它们报告给您的图书馆供应商，以便他们得到修复！

在依赖精度的情况下，是否存在确保使用真正的 FMA 的内在因素？

考虑到一个好的编译器，这不应该是必要的；使用fma( ) 函数并告诉编译器您的目标是什么架构就足够了。但是，编译器并不完美，因此您可能需要在 x86 上使用 _mm_fmadd_sd( ) 和相关的内在函数（但请向您的编译器供应商报告错误！）

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

这基本上是从(ahi,alo)*(bhi,blo) = (ahi*bhi + ahi*blo + alo*bhi + alo*blo) 中减去c。

我从论文Extended-Precision Floating-Point Numbers for GPU Computation 中的twoProd 函数和Agner Fog's vector class library 中的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};
}

在任何情况下，在fmsub 中使用split 或split2 与glibc 中数学库中的fma(a,b,-c) 非常一致。无论出于何种原因，我的版本都比 fma 快得多，除非在具有硬件 fma 的机器上（在这种情况下，我还是使用 _mm_fmsub_ss）。

Answer 3

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

因此，虽然 Dekker 的 twoProduct 可以通过硬件 FMA 大大加速，但 Dekker 的 twoProduct 的误差项计算不是稳健的 FMA 实现。

正确的实现需要使用高于双精度的求和算法，或者按数量级降序添加项。