我的fma()坏了吗？

Question

在使用中,double fma(double x, double y, double z);我希望d下面标有的输出行中有一个非零值'?'.它似乎在内部只使用long double精度,而不是无限精度的规定.

的fma函数计算(x× y)+ z,四舍五入一个三元操作:它们计算的值(仿佛)到无限精度和圆一次的结果的格式,根据当前的舍入模式.§7.12.13.12(我的重点)

那我的fma()坏了,或者我在代码或编译选项中如何错误地使用它？

#include <float.h>
#include <math.h>
#include <stdio.h>

int main(void) {
  // Invoking: Cygwin C Compiler
  // gcc -std=c11 -O0 -g3 -pedantic -Wall -Wextra -Wconversion -c -fmessage-length=0 
  //   -v -MMD -MP -MF"x.d" -MT"x.o" -o "x.o" "../x.c"

  printf("FLT_EVAL_METHOD %d\n", FLT_EVAL_METHOD);
  for (unsigned i = 20; i < 55; i++) {
    volatile double a = 1.0 + 1.0 / pow(2, i);
    volatile double b = a;
    volatile double c = a * b;
    volatile double d = fma(a, b, -c);
    volatile char *nz = ((i >= 27 && a != 1.0) == !d) ? "?" : "";
    printf("i:%2u a:%21.13a c:%21.13a d:%10a %s\n", i, a, c, d, nz);
  }
  return 0;
}

产量

FLT_EVAL_METHOD 2
i:20 a: 0x1.0000100000000p+0 c: 0x1.0000200001000p+0 d:    0x0p+0 
i:21 a: 0x1.0000080000000p+0 c: 0x1.0000100000400p+0 d:    0x0p+0 
i:22 a: 0x1.0000040000000p+0 c: 0x1.0000080000100p+0 d:    0x0p+0 
i:23 a: 0x1.0000020000000p+0 c: 0x1.0000040000040p+0 d:    0x0p+0 
i:24 a: 0x1.0000010000000p+0 c: 0x1.0000020000010p+0 d:    0x0p+0 
i:25 a: 0x1.0000008000000p+0 c: 0x1.0000010000004p+0 d:    0x0p+0 
i:26 a: 0x1.0000004000000p+0 c: 0x1.0000008000001p+0 d:    0x0p+0 
i:27 a: 0x1.0000002000000p+0 c: 0x1.0000004000000p+0 d:   0x1p-54 
i:28 a: 0x1.0000001000000p+0 c: 0x1.0000002000000p+0 d:   0x1p-56 
i:29 a: 0x1.0000000800000p+0 c: 0x1.0000001000000p+0 d:   0x1p-58 
i:30 a: 0x1.0000000400000p+0 c: 0x1.0000000800000p+0 d:   0x1p-60 
i:31 a: 0x1.0000000200000p+0 c: 0x1.0000000400000p+0 d:   0x1p-62 
i:32 a: 0x1.0000000100000p+0 c: 0x1.0000000200000p+0 d:    0x0p+0 ?
i:33 a: 0x1.0000000080000p+0 c: 0x1.0000000100000p+0 d:    0x0p+0 ?
i:34 a: 0x1.0000000040000p+0 c: 0x1.0000000080000p+0 d:    0x0p+0 ?
...
i:51 a: 0x1.0000000000002p+0 c: 0x1.0000000000004p+0 d:    0x0p+0 ?
i:52 a: 0x1.0000000000001p+0 c: 0x1.0000000000002p+0 d:    0x0p+0 ?
i:53 a: 0x1.0000000000000p+0 c: 0x1.0000000000000p+0 d:    0x0p+0 
i:54 a: 0x1.0000000000000p+0 c: 0x1.0000000000000p+0 d:    0x0p+0 

版本信息

gcc -v

Using built-in specs.
COLLECT_GCC=gcc
COLLECT_LTO_WRAPPER=/usr/lib/gcc/i686-pc-cygwin/5.3.0/lto-wrapper.exe
Target: i686-pc-cygwin
Configured with: /cygdrive/i/szsz/tmpp/gcc/gcc-5.3.0-5.i686/src/gcc-5.3.0/configure --srcdir=/cygdrive/i/szsz/tmpp/gcc/gcc-5.3.0-5.i686/src/gcc-5.3.0 --prefix=/usr --exec-prefix=/usr --localstatedir=/var --sysconfdir=/etc --docdir=/usr/share/doc/gcc --htmldir=/usr/share/doc/gcc/html -C --build=i686-pc-cygwin --host=i686-pc-cygwin --target=i686-pc-cygwin --without-libiconv-prefix --without-libintl-prefix --libexecdir=/usr/lib --enable-shared --enable-shared-libgcc --enable-static --enable-version-specific-runtime-libs --enable-bootstrap --enable-__cxa_atexit --with-dwarf2 --with-arch=i686 --with-tune=generic --disable-sjlj-exceptions --enable-languages=ada,c,c++,fortran,java,lto,objc,obj-c++ --enable-graphite --enable-threads=posix --enable-libatomic --enable-libcilkrts --enable-libgomp --enable-libitm --enable-libquadmath --enable-libquadmath-support --enable-libssp --enable-libada --enable-libjava --enable-libgcj-sublibs --disable-java-awt --disable-symvers --with-ecj-jar=/usr/share/java/ecj.jar --with-gnu-ld --with-gnu-as --with-cloog-include=/usr/include/cloog-isl --without-libiconv-prefix --without-libintl-prefix --with-system-zlib --enable-linker-build-id --with-default-libstdcxx-abi=gcc4-compatible
Thread model: posix
gcc version 5.3.0 (GCC) 

Answer 1

这是Cygwin的错.或者更确切地说,它使用的newlib C库.它明确表示它甚至不试图让fma()仿真正确.

自2015年以来,GNU C库对几乎所有fma变体都有正确的模拟.有关详细信息以及用于实现此功能的补丁,请参阅源软件错误13304.

如果效率不是问题,那么我只想使用eg

#if defined(__CYGWIN__) && !defined(__FMA__) && !defined(__FMA3__) && !defined(__FMA4__)
#define fma(x, y, z)  fma_emulation(x, y, z)

double fma_emulation(double x, double y, double z)
{
    /* One of the implementations linked above */
}
#endif

我个人根本不使用Windows,但如果有人这样做(使用Windows并需要fma仿真),我建议他们尝试向上游推送补丁,并提供关于正确fma仿真的GNU C库讨论的链接.

我想知道的是,是否有可能只检查结果的低M位(在舍入中丢弃)以确定结果中ULP的正确值,并调整使用直接a × b获得的结果+ c操作相应,使用nextafter(); 而不是使用multiprecision算法来实现整个操作.

编辑:不,因为添加可能会溢出,丢弃一个额外的位作为丢弃部分的MSB.仅仅因为这个原因,我们确实需要完成整个操作.另一个原因是,如果a × b和c具有不同的符号,那么我们从幅度越大(使用更大的符号的结果)中减去幅度越小而不是相加,这可能导致清除更大的几个高位,并且影响舍入中丢弃整个结果的哪些位.

但是,对于doublex86和x86-64架构上的IEEE-754 Binary64 ,我确实认为使用64位(整数)寄存器和128位产品进行fma仿真仍然是非常可行的.我将尝试一种表示,其中64位寄存器中的低2位用于舍入决策位(LSB是所有丢弃位的逻辑或),53位用于尾数,以及一个进位,留下8位未使用和忽略的高位.当无符号整数尾数转换为(64位)双精度时执行舍入.如果这些实验产生任何有用的东西,我将在这里描述它们.

初步调查结果:fma()32位系统上的仿真速度很慢.387 FPU上的80位内容在这里基本没用,在32位系统上实现53×53位乘法(和位移)只是不值得的.fma()在我看来,与上面链接的glibc 仿真代码已经足够了.

其他发现:处理非有限值是令人讨厌的.(次正规只是稍微烦人,需要特殊处理(因为尾数中的隐式MSB为零).)如果三个参数中的任何一个是非有限的(无穷大或某种形式的NaN),则返回a*b + c(未融合)是只有理智的选择.处理这些情况需要额外的分支,这会减慢仿真速度.

最终决定:以优化方式处理的案例数量(而不是使用glibc仿真中使用的多精度"肢体"方法)足以使这种方法不值得付出努力.如果每个肢体是64位,则a,b和c中的每一个分布在最多2个肢体上,并且a × b在三个肢体上.(对于32位肢体,分别只有3和5个肢体.)根据a × b和c是否具有相同或不同的符号,只有两种根本不同的情况需要处理 - 在不同的标志情况下,加法变为减法(从较大的较小,结果得到与较大值相同的符号).

简而言之,多精度方法更好.所需的实际精度非常有限,甚至不需要动态分配.如果可以有效地计算a和b的尾数的乘积,则多精度部分限于保持产品并处理加法/减法.最后的舍入可以通过将结果转换为53位尾数,指数和两个额外的低位来实现(较高的是舍入中丢失的最高位,而较低位是其余位丢失的OR).四舍五入).本质上,关键操作可以使用整数(或SSE/AVX寄存器)完成,最终从55位尾数转换为double,根据当前规则处理舍入.