the*_*ine 5 floating-point instruction-set ieee-754 fma
根据文档,有一个fma()功能math.h.这非常好,我知道FMA如何工作以及如何使用它.但是,我不太确定这在实践中如何实施？我最感兴趣的是x86和x86_64架构.
是否存在FMA的浮点(非向量)指令,可能是IEEE-754 2008定义的？
是使用FMA3还是FMA4指令？
在依赖精度的情况下,是否存在确保使用真实FMA的内在因素？
实际实施因平台而异,但讲得非常广泛:
如果您告诉编译器使用硬件FMA指令(PowerPC,带有VFPv4或AArch64,Intel Haswell或AMD Bulldozer及其后的ARM)定位计算机,编译器可以fma( )通过将适当的指令放入代码中来替换调用.这不是保证,但通常是良好的做法.否则,您将调用数学库,并且:
在具有硬件FMA的处理器上运行时,应使用这些指令来实现该功能.但是,如果您的操作系统版本较旧,或者数学库版本较旧,则可能无法利用这些说明.
如果您运行的是没有硬件FMA的处理器,或者您使用的是较旧的(或者不是很好的)数学库,那么将使用FMA的软件实现.这可以使用巧妙的扩展精度浮点技巧或整数运算来实现.
该fma( )函数的结果应始终正确舍入(即"真正的外汇").如果不是,那就是系统数学库中的错误.不幸的是,fma( )正确实现的数学库函数是一个比较困难的,所以很多实现都有bug.请将它们报告给您的图书馆供应商,以便修复它们!
在依赖精度的情况下,是否存在确保使用真实FMA的内在因素？
鉴于良好的编译器,这不应该是必要的; 它应该足以使用该fma( )函数并告诉编译器您要定位的架构.但是,编译器并不完美,因此您可能需要_mm_fmadd_sd( )在x86上使用相关的内在函数(但是将错误报告给编译器供应商!)
在软件中实现 FMA 的一种方法是将有效位分为高位和低位。我使用Dekker 算法
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论文Extended-Precision Floating-Point Numbers for GPU Computation 中的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无论如何都使用）。
小智 5
不幸的是，Z 玻色子基于 Dekker 算法的 FMA 建议是不正确的。与 Dekker 的 twoProduct 不同，在更一般的 FMA 情况下，c 的大小相对于乘积项是未知的，因此可能会发生错误的抵消。
因此，虽然 Dekker 的 twoProduct 可以通过硬件 FMA 大大加速，但 Dekker 的 twoProduct 的误差项计算并不是一个健壮的 FMA 实现。
正确的实现需要使用高于双精度的求和算法，或者以数量级递减的顺序添加项。