我的fma（）坏了吗？
在使用double fma(double x, double y, double z); 我希望下面的输出行中标有'?'的非零d 。 内部似乎只使用long double精度而不是指定的无限精度 。
fma函数计算（x×y）+z，舍入为一个三元运算：根据当前舍入模式，它们将值（如同）计算为无限精度并对结果格式舍入一次。 §7.12.13.12（我的重点）
那么我的fma()坏了，或者我在代码或编译选项中如何错误地使用它？
#include #include #include int main(void) { // Invoking: Cygwin C Compiler // gcc -std=c11 -O0 -g3 -pedantic -Wall -Wextra -Wconversion -c -fmessage-length=0 // -v -MMD -MP -MF"xd" -MT"xo" -o "xo" "../xc" printf("FLT_EVAL_METHOD %d\n", FLT_EVAL_METHOD); for (unsigned i = 20; i = 27 && a != 1.0) == !d) ? "?" : ""; printf("i:%2u a:%21.13ac:%21.13ad:%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)
这是Cygwin的错。 或者更确切地说，它使用的newlib C库。 它明确表示它甚至不会尝试使fma()仿真正确。
自2015年以来，GNU C库对几乎所有fma变体都有正确的模拟。有关详细信息以及用于实现此function的补丁，请参阅源软件错误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仿真），我建议他们尝试上游推送补丁，并链接到GNU C库讨论正确的fma仿真 。
我想知道的是，是否有可能只检查结果的低M位（在舍入中丢弃）以确定结果中ULP的正确值，并调整使用直接a × b获得的结果+ c操作相应，使用nextafter() ; 而不是使用multiprecision算法来实现整个操作。
编辑：不，因为添加可能会溢出，丢弃额外的位作为丢弃部分的MSB。 仅仅因为这个原因，我们确实需要完成整个操作。 另一个原因是，如果a × b和c具有不同的符号，那么我们从幅度越大（使用更大的符号的结果）减去幅度越小而不是加法，这可能导致清除更大的几个高位，并且影响舍入中丢弃整个结果的哪些位。
但是，对于x86和x86-64架构上的IEEE-754 Binary64 double ，我确实认为使用64位（整数）寄存器和128位产品的fma仿真仍然是非常可行的。 我将试验一个表示，其中64位寄存器中的低2位用于舍入决策位（LSB是所有丢弃位的逻辑或），53位用于尾数，一个进位，留下8位未使用和忽略的高位。 当无符号整数尾数转换为（64位）双精度时执行舍入。 如果这些实验产生任何有用的东西，我将在这里描述它们。
初步调查结果：32位系统上的fma()仿真速度很慢。 387 FPU上的80位内容在这里基本没用，在32位系统上实现53×53位乘法（和位移）只是不值得的。 在我看来，与上面链接的glibc fma()仿真代码已经足够了。
其他发现：处理非有限值是令人讨厌的 。 （次正规只是稍微烦人，需要特殊处理（因为尾数中的隐式MSB为零）。）如果三个参数中的任何一个是非有限的（无穷大或某种forms的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的最终转换根据当前规则处理舍入。