Неужели мой fma() сломан?

Question

При использовании double fma(double x, double y, double z); я бы ожидал ненулевого d в выходных строках ниже, обозначенных '?'. Это появляется для внутреннего использования long double точность, а не бесконечная точность как указано.

В fma функции вычисления (x × y) + z, закругленная в качестве одного из трехкомпонентной операции: они вычисляют значение (как если бы), чтобы бесконечной точности и круглых один раз в формат результат, в соответствии с текущим округлением Режим. §7.12.13.1 2 (курсив мой)

Так мой 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, которую он использует. Это explicitly says он даже не пытается получить fma() эмуляция правильно.

Библиотека GNU C имеет правильную эмуляцию практически для всех вариантов fma с 2015 года. Подробности и исправления, используемые для реализации этого, см. В исходной ошибке 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), я хотел бы предложить, они пытаются подтолкнуть патч выше, со ссылкой на GNU C library discussion on correct fma emulation.

---

что я интересно, это будет ли это можно рассматривать только низкие M битые результат (отбрасывает в округлении), чтобы определить правильное значение ULP в результате, а также настроить результат, полученный с использованием простого a × b + c эксплуатация соответственно с использованием nextafter(); вместо использования арифметики multiprecision для реализации всей операции.

Редактировать: Нет, потому что добавление может переполняться, отбрасывая дополнительный бит в качестве MSB отбрасываемой части. По этой причине нам нужно выполнить всю операцию. Другая причина заключается в том, что если a × b и c имеют разные знаки, то вместо добавления мы вычитаем меньше по величине от большего по величине (результат использует знак большего), что может привести к очистке нескольких старших бит в большем, и это влияет на то, какие биты всего результата выпадают в округлении.

Однако для IEEE-754 Binary64 double на архитектурах x86 и x86-64 я действительно верю, что эмуляция fma с использованием 64-разрядных (целочисленных) регистров и 128-битного продукта по-прежнему вполне возможна. Я буду экспериментировать с представлением, в котором для битов решения округления используются младшие 2 бита в 64-битном регистре (LSB является логическим ИЛИ всех сброшенных битов), 53 бит, используемых для мантиссы, и один бит переноса, оставляя 8 неиспользуемые и игнорируемые высокие биты. Округление выполняется, когда целочисленная мантисса без знака преобразуется в (64-разрядный) двойной. Если эти эксперименты дают что-нибудь полезное, я опишу их здесь.

---

Исходные данные: fma() эмуляция на 32-битной системе медленно. 80-битный материал на FPU 387 в основном бесполезен здесь, и реализация 53 × 53-битного умножения (и смещения битов) на 32-битной системе просто ... не стоит усилий. На мой взгляд, код эмуляции glibc fma() достаточно хорош.

Дополнительные выводы: Обработка не конечных значений: nasty. (Subnormals только немного раздражают, требуя специальной обработки (поскольку неявный MSB в мантиссе равен нулю тогда).) Если какой-либо из трех аргументов не является конечным (бесконечность или некоторая форма NaN), то возврат a*b + c (не скомпоновано) единственный нормальный вариант. Для обработки этих случаев требуется дополнительное разветвление, что замедляет эмуляцию.

Окончательное решение: количество обработанных случаев оптимизированным образом (а не использование подхода «конечности» с несколькими точками, используемого в эмуляции glibc) достаточно велико, чтобы этот подход не стоил усилий. Если каждая конечность 64-бит, каждый из , б и с распространяется на самых большие 2 конечностей, а × б в течение трех конечностей. (С 32-битными конечностями, то есть всего 3 и 5 конечностей соответственно). В зависимости от того, имеют ли a × × b и c имеют одинаковые или разные знаки, есть только два принципиально разных случая - в случае с разными знаками добавление превращается в вычитание (меньше от большего, результат получает тот же знак, что и большее значение).

Короче говоря, подход с множественными точками лучше. Необходимая точность очень ограничена и не требует даже динамического распределения. Если произведение мантиссы a и b может быть рассчитано эффективно, часть с высокой точностью ограничена удерживанием продукта и обработкой сложения/вычитания. Окончательное округление может быть выполнено путем преобразования результата в 53-битную мантиссу, экспоненту и два дополнительных младших бита (более высокий - самый старший бит, потерянный в округлении, а нижний - OR остальных остатков, потерянных в округление). По сути, ключевые операции могут выполняться с использованием целых чисел (или регистров SSE/AVX), а окончательное преобразование из 55-битной мантиссы для двойного обращения округления в соответствии с текущими правилами.