fma() 是如何实现的

根据文档，math.h中有一个fma()函数。这非常好，我知道FMA是如何工作的以及使用它做什么。但是，我不太确定这在实践中是如何实施的？我最感兴趣的是x86和x86_64架构。

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

是否使用 FMA3 或 FMA4 指令？

当依赖精度时，是否有内在因素可以确保使用真正的 FMA？

实际实现因平台而异，但讲义非常广泛：

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

当依赖精度时，是否有内在因素可以确保使用真正的 FMA？

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

在软件中实现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;
}

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

我从论文 GPU 计算的扩展精度浮点数中的twoProd函数和 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};
}

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

不幸的是，

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

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

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

相关内容

最新更新

热门标签：