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( ) 函数的结果应始终正确舍入（即 "real fma"）。如果不是，那是你系统数学库中的错误。不幸的是，fma( ) 是更难正确实现的数学库函数之一，因此许多实现都存在错误。请将它们报告给您的图书馆供应商，以便他们得到修复！

Is there an intrinsic to make sure that a real FMA is used, when the precision is relied upon?

如果有一个好的编译器，这应该不是必需的；使用 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 and from the mul_sub_x function in Agner Fog's vector class library中的twoProd函数得到的。他使用不同的函数来拆分以不同方式拆分的浮点向量。我试图在此处重现标量版本

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 玻色子基于 Dekker 算法的 FMA 建议是不正确的。与 Dekker 的 twoProduct 不同，在更一般的 FMA 情况下，c 的大小相对于乘积项是未知的，因此可能会发生错误的取消。

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

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

fma() 是如何实现的

How is fma() implemented

floating-point

ieee-754

instruction-set

fma