Est-ce que mon fma () est cassé?

En utilisant double fma(double x, double y, double z); Je m’attendrais à un non-zéro d dans les lignes de sortie ci-dessous marqués d’un '?' . Il semble que, de manière interne, elle utilise uniquement la long double précision long double plutôt que la précision infinie spécifiée.

Les fonctions fma calculent ( x × y ) + z , arrondies comme une opération ternaire: elles calculent la valeur (comme si) avec une précision infinie et arrondissent une fois au format de résultat, en fonction du mode d’arrondi actuel. §7.12.13.1 2 (mon emphase)

Donc, mon fma() cassé, ou comment est-ce que je l’utilise incorrectement dans le code ou les options de compilation?

 #include  #include  #include  int main(void) { // Invoking: Cygwin C Comstackr // 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; } 

Sortie

 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 

Informations de version

 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) 

C’est la faute de Cygwin. Ou plus exactement, la bibliothèque newlib C qu’elle utilise. Il dit explicitement qu’il n’essaie même pas d’obtenir une émulation fma() correcte.

La bibliothèque GNU C a une émulation correcte pour presque toutes les variantes de fma depuis 2015. Pour plus de détails et pour les correctifs utilisés, reportez-vous au bogue sourceware 13304 .

Si l’efficacité n’est pas un problème, j’utiliserais simplement par exemple

 #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 

Personnellement, je n’utilise pas du tout Windows, mais si quelqu’un le fait (utilise Windows et a besoin de l’émulation fma), je suggérerais d’essayer de placer un correctif en amont, avec un lien vers la discussion de la bibliothèque GNU C sur l’émulation fma correcte .


Ce que je me demande, est de savoir s’il serait possible d’examiner uniquement les M bits les plus faibles du résultat (mis au rebut lors de l’arrondi) afin de déterminer la valeur correcte du paramètre ULP dans le résultat et d’ajuster le résultat obtenu à l’aide du simple a × b. + c opération en conséquence, en utilisant nextafter() ; plutôt que d’utiliser l’arithmétique multi-précision pour mettre en œuvre l’ensemble de l’opération.

Edit: Non, car l’ajout risque de déborder, laissant tomber un bit supplémentaire en tant que MSB de la partie mise au rebut. Pour cette seule raison, nous devons effectuer toute l’opération. Une autre raison est que si a × b et c ont des signes différents, alors, au lieu d’addition, nous soustrayons les plus petites magnitudes des plus grandes (résultats utilisant le signe du plus grand), ce qui peut entraîner la suppression de plusieurs bits élevés dans les plus grands. affecte quels bits de l’ensemble du résultat sont supprimés dans l’arrondi.

Cependant, pour le double IEEE-754 Binary64 sur les architectures x86 et x86-64, je pense qu’une émulation fma utilisant des registres 64 bits (entier) et un produit 128 bits est encore tout à fait réalisable. J’expérimenterai une représentation où 2 bits bas dans un registre de 64 bits sont utilisés pour les bits de décision d’arrondi (LSB étant un OU logique de tous les bits supprimés), 53 bits utilisés pour la mantisse et un bit de retenue, laissant 8 bits élevés non utilisés et ignorés. L’arrondi est effectué lorsque la mantisse entière non signée est convertie en un double (64 bits). Si ces expériences donnent des résultats utiles, je les décrirai ici.


fma() initiales: l’émulation fma() sur un système 32 bits est lente. Les éléments 80 bits de la FPU 387 sont fondamentalement inutiles, et implémenter la multiplication 53 × 53 bits (et le décalage de bits) sur un système 32 bits ne vaut tout simplement pas la peine. Le code d’émulation glibc fma() lié à ci-dessus est assez bon à mon avis.

Constatations supplémentaires: Le traitement des valeurs non finies est désagréable . (Les subnormals ne sont que légèrement gênants et nécessitent un traitement spécial (car le MSB implicite dans la mantisse vaut alors zéro).) Si l’un des trois arguments est non fini (infini ou une forme quelconque de NaN), il retourne alors a*b + c ( non fusionné) est la seule option saine. La gestion de ces cas nécessite des twigs supplémentaires, ce qui ralentit l’émulation.

Décision finale: le nombre de cas à traiter de manière optimisée (plutôt que d’utiliser l’approche multi-précision “membre” utilisée dans l’émulation glibc) est suffisamment important pour que cette approche ne vaille pas la peine. Si chaque membre a 64 bits, chacun de a , b et c est réparti sur au plus 2 membres et a × b sur trois membres. (Avec les membres 32 bits, cela ne représente que 3 et 5 membres, respectivement.) Selon que a × b et c ont des signes identiques ou différents, il n’ya que deux cas fondamentalement différents à traiter – dans le cas des signes différents, l’addition se transforme en soustraction (plus petit plus grand, résultat obtenant le même signe que la plus grande valeur).

En bref, l’approche multi-précision est préférable. La précision requirejse est très bien délimitée et ne nécessite pas même une allocation dynamic. Si le produit des mantisses de a et b peut être calculé efficacement, la partie multi-précision est limitée à la détention du produit et à la gestion de l’addition / soustraction. L’arrondi final peut être effectué en convertissant le résultat en une mantisse de 53 bits, un exposant et deux bits extra-bas (le plus élevé étant le bit le plus significatif perdu lors de l’arrondi et le plus petit un OR des autres bits perdus dans l’arrondi). l’arrondi). Pour l’essentiel, les opérations sur les touches peuvent être effectuées à l’aide d’entiers (ou de registres SSE / AVX), puis la conversion finale d’une mantisse 55 bits permet de gérer l’arrondi conformément aux règles en vigueur.