我的 fma() 坏了吗？
Is my fma() broken?
在使用 double fma(double x, double y, double z); 时，我希望在下面标有 '?' 的输出行中出现非零值 d。它似乎在内部仅使用long double精度而不是指定的无限精度。
The fma functions compute (x × y) + z, rounded as one ternary operation: they compute the value (as if) to infinite precision and round once to the result format, according to the current rounding mode. §7.12.13.1 2 (my emphasis)
那么我的 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)
这是 Cygwin 的错。或者更准确地说，它使用的 newlib C 库。它 explicitly says 它甚至不尝试使 fma() 模拟正确。
自 2015 年以来，GNU C 库对几乎所有 fma 变体都有正确的仿真。有关详细信息以及用于实现此功能的补丁，请参阅源软件错误 13304。
如果效率不是问题，那么我会简单地使用例如
#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 仿真），我建议他们尝试向上游推送一个补丁， link 到 GNU C library discussion on correct fma emulation.
我想知道的是，是否可以仅检查结果的低 M 位（在舍入中舍弃）以确定 ULP 的正确值在结果中，并相应地调整使用直接 a×b+c 操作获得的结果，使用 nextafter()；而不是使用多精度算法来实现整个操作。
编辑：不，因为加法可能会溢出，丢弃一个额外的位作为丢弃部分的 MSB。仅出于这个原因，我们确实需要进行整个操作。另一个原因是，如果 a×b 和 c 符号不同，那么我们不加法，而是从较大的幅度中减去较小的幅度（结果使用较大的符号），这可能会导致清除较大的几个高位，这会影响整个结果的哪些位在舍入中被丢弃。
但是，对于 x86 和 x86-64 架构上的 IEEE-754 Binary64 double，我相信使用 64 位（整数）寄存器和 128 位乘积的 fma 仿真仍然非常可行。我将试验一种表示，其中 64 位寄存器中的低 2 位用于舍入决策位（LSB 是所有丢弃位的逻辑或），53 位用于尾数，一个进位，剩下 8未使用和忽略的高位。当无符号整数尾数转换为（64 位）双精度数时执行舍入。如果这些实验产生了什么有用的东西，我会在这里描述它们。
初步发现：fma() 32 位系统上的仿真速度很慢。 387 FPU 上的 80 位东西在这里基本上没有用，在 32 位系统上实现 53×53 位乘法（和位移）只是..不值得付出努力。 glibc fma() 仿真代码 link 在我看来已经足够好了。
其他发现：处理 non-finite 值是 令人讨厌的 。 （次正规只是有点烦人，需要特殊处理（因为尾数中的隐式 MSB 为零）。）如果三个参数中的任何一个是 non-finite（无穷大或某种形式的 NaN），则返回 a*b + c （未融合）是唯一明智的选择。处理这些情况需要额外的分支，这会减慢仿真速度。
最终决定：以优化方式处理的案例数量（而不是使用 glibc 仿真中使用的多精度 "limb" 方法）大到足以使这种方法不值得付出努力。如果每个分支都是 64 位的，则 a、b 和 c 中的每一个分布在大多数 2 个肢体，a×b 在三个肢体上。 （对于 32 位肢体，分别只有 3 个肢体和 5 个肢体。）取决于 a×b 和 c 有相同或不同的符号，只有两种根本不同的情况要处理——在不同符号的情况下，加法变成减法（从大到小，结果得到与较大值相同的符号） .
简而言之，多精度方法更好。所需的实际精度非常有限，甚至不需要动态分配。如果a和b的尾数的乘积可以高效计算，多精度部分仅限于保存乘积和处理addition/subtraction.最后舍入可以通过将结果转换为 53 位尾数、指数和两个额外的低位来完成（较高的是舍入中丢失的最高有效位，较低的是舍入中丢失的其余位的或）四舍五入）。本质上，关键操作可以使用整数（或 SSE/AVX 寄存器）完成，并且从 55 位尾数到双精度的最终转换根据当前规则处理舍入。
在使用 double fma(double x, double y, double z); 时，我希望在下面标有 '?' 的输出行中出现非零值 d。它似乎在内部仅使用long double精度而不是指定的无限精度。
The
fmafunctions compute (x×y) +z, rounded as one ternary operation: they compute the value (as if) to infinite precision and round once to the result format, according to the current rounding mode. §7.12.13.1 2 (my emphasis)
那么我的 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)
这是 Cygwin 的错。或者更准确地说，它使用的 newlib C 库。它 explicitly says 它甚至不尝试使 fma() 模拟正确。
自 2015 年以来，GNU C 库对几乎所有 fma 变体都有正确的仿真。有关详细信息以及用于实现此功能的补丁，请参阅源软件错误 13304。
如果效率不是问题，那么我会简单地使用例如
#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 仿真），我建议他们尝试向上游推送一个补丁， link 到 GNU C library discussion on correct fma emulation.
我想知道的是，是否可以仅检查结果的低 M 位（在舍入中舍弃）以确定 ULP 的正确值在结果中，并相应地调整使用直接 a×b+c 操作获得的结果，使用 nextafter()；而不是使用多精度算法来实现整个操作。
编辑：不，因为加法可能会溢出，丢弃一个额外的位作为丢弃部分的 MSB。仅出于这个原因，我们确实需要进行整个操作。另一个原因是，如果 a×b 和 c 符号不同，那么我们不加法，而是从较大的幅度中减去较小的幅度（结果使用较大的符号），这可能会导致清除较大的几个高位，这会影响整个结果的哪些位在舍入中被丢弃。
但是，对于 x86 和 x86-64 架构上的 IEEE-754 Binary64 double，我相信使用 64 位（整数）寄存器和 128 位乘积的 fma 仿真仍然非常可行。我将试验一种表示，其中 64 位寄存器中的低 2 位用于舍入决策位（LSB 是所有丢弃位的逻辑或），53 位用于尾数，一个进位，剩下 8未使用和忽略的高位。当无符号整数尾数转换为（64 位）双精度数时执行舍入。如果这些实验产生了什么有用的东西，我会在这里描述它们。
初步发现：fma() 32 位系统上的仿真速度很慢。 387 FPU 上的 80 位东西在这里基本上没有用，在 32 位系统上实现 53×53 位乘法（和位移）只是..不值得付出努力。 glibc fma() 仿真代码 link 在我看来已经足够好了。
其他发现：处理 non-finite 值是 令人讨厌的 。 （次正规只是有点烦人，需要特殊处理（因为尾数中的隐式 MSB 为零）。）如果三个参数中的任何一个是 non-finite（无穷大或某种形式的 NaN），则返回 a*b + c （未融合）是唯一明智的选择。处理这些情况需要额外的分支，这会减慢仿真速度。
最终决定：以优化方式处理的案例数量（而不是使用 glibc 仿真中使用的多精度 "limb" 方法）大到足以使这种方法不值得付出努力。如果每个分支都是 64 位的，则 a、b 和 c 中的每一个分布在大多数 2 个肢体，a×b 在三个肢体上。 （对于 32 位肢体，分别只有 3 个肢体和 5 个肢体。）取决于 a×b 和 c 有相同或不同的符号，只有两种根本不同的情况要处理——在不同符号的情况下，加法变成减法（从大到小，结果得到与较大值相同的符号） .
简而言之，多精度方法更好。所需的实际精度非常有限，甚至不需要动态分配。如果a和b的尾数的乘积可以高效计算，多精度部分仅限于保存乘积和处理addition/subtraction.最后舍入可以通过将结果转换为 53 位尾数、指数和两个额外的低位来完成（较高的是舍入中丢失的最高有效位，较低的是舍入中丢失的其余位的或）四舍五入）。本质上，关键操作可以使用整数（或 SSE/AVX 寄存器）完成，并且从 55 位尾数到双精度的最终转换根据当前规则处理舍入。