soft-fp/op-common.h

   1 /* Software floating-point emulation. Common operations.
   2    Copyright (C) 1997-2014 Free Software Foundation, Inc.
   3    This file is part of the GNU C Library.
   4    Contributed by Richard Henderson (rth@cygnus.com),
   5                   Jakub Jelinek (jj@ultra.linux.cz),
   6                   David S. Miller (davem@redhat.com) and
   7                   Peter Maydell (pmaydell@chiark.greenend.org.uk).
   8
   9    The GNU C Library is free software; you can redistribute it and/or
  10    modify it under the terms of the GNU Lesser General Public
  11    License as published by the Free Software Foundation; either
  12    version 2.1 of the License, or (at your option) any later version.
  13
  14    In addition to the permissions in the GNU Lesser General Public
  15    License, the Free Software Foundation gives you unlimited
  16    permission to link the compiled version of this file into
  17    combinations with other programs, and to distribute those
  18    combinations without any restriction coming from the use of this
  19    file.  (The Lesser General Public License restrictions do apply in
  20    other respects; for example, they cover modification of the file,
  21    and distribution when not linked into a combine executable.)
  22
  23    The GNU C Library is distributed in the hope that it will be useful,
  24    but WITHOUT ANY WARRANTY; without even the implied warranty of
  25    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  26    Lesser General Public License for more details.
  27
  28    You should have received a copy of the GNU Lesser General Public
  29    License along with the GNU C Library; if not, see
  30    <http://www.gnu.org/licenses/>.  */
  31
  32 #define _FP_DECL(wc, X)                         \
  33   _FP_I_TYPE X##_c __attribute__ ((unused));    \
  34   _FP_I_TYPE X##_s __attribute__ ((unused));    \
  35   _FP_I_TYPE X##_e __attribute__ ((unused));    \
  36   _FP_FRAC_DECL_##wc (X)
  37
  38 /* Test whether the qNaN bit denotes a signaling NaN.  */
  39 #define _FP_FRAC_SNANP(fs, X)                           \
  40   ((_FP_QNANNEGATEDP)                                   \
  41    ? (_FP_FRAC_HIGH_RAW_##fs (X) & _FP_QNANBIT_##fs)    \
  42    : !(_FP_FRAC_HIGH_RAW_##fs (X) & _FP_QNANBIT_##fs))
  43 #define _FP_FRAC_SNANP_SEMIRAW(fs, X)                   \
  44   ((_FP_QNANNEGATEDP)                                   \
  45    ? (_FP_FRAC_HIGH_##fs (X) & _FP_QNANBIT_SH_##fs)     \
  46    : !(_FP_FRAC_HIGH_##fs (X) & _FP_QNANBIT_SH_##fs))
  47
  48 /*
  49  * Finish truly unpacking a native fp value by classifying the kind
  50  * of fp value and normalizing both the exponent and the fraction.
  51  */
  52
  53 #define _FP_UNPACK_CANONICAL(fs, wc, X)                         \
  54   do                                                            \
  55     {                                                           \
  56       switch (X##_e)                                            \
  57         {                                                       \
  58         default:                                                \
  59           _FP_FRAC_HIGH_RAW_##fs (X) |= _FP_IMPLBIT_##fs;       \
  60           _FP_FRAC_SLL_##wc (X, _FP_WORKBITS);                  \
  61           X##_e -= _FP_EXPBIAS_##fs;                            \
  62           X##_c = FP_CLS_NORMAL;                                \
  63           break;                                                \
  64                                                                 \
  65         case 0:                                                 \
  66           if (_FP_FRAC_ZEROP_##wc (X))                          \
  67             X##_c = FP_CLS_ZERO;                                \
  68           else                                                  \
  69             {                                                   \
  70               /* a denormalized number */                       \
  71               _FP_I_TYPE _shift;                                \
  72               _FP_FRAC_CLZ_##wc (_shift, X);                    \
  73               _shift -= _FP_FRACXBITS_##fs;                     \
  74               _FP_FRAC_SLL_##wc (X, (_shift+_FP_WORKBITS));     \
  75               X##_e -= _FP_EXPBIAS_##fs - 1 + _shift;           \
  76               X##_c = FP_CLS_NORMAL;                            \
  77               FP_SET_EXCEPTION (FP_EX_DENORM);                  \
  78             }                                                   \
  79           break;                                                \
  80                                                                 \
  81         case _FP_EXPMAX_##fs:                                   \
  82           if (_FP_FRAC_ZEROP_##wc (X))                          \
  83             X##_c = FP_CLS_INF;                                 \
  84           else                                                  \
  85             {                                                   \
  86               X##_c = FP_CLS_NAN;                               \
  87               /* Check for signaling NaN */                     \
  88               if (_FP_FRAC_SNANP (fs, X))                       \
  89                 FP_SET_EXCEPTION (FP_EX_INVALID);               \
  90             }                                                   \
  91           break;                                                \
  92         }                                                       \
  93     }                                                           \
  94   while (0)
  95
  96 /* Finish unpacking an fp value in semi-raw mode: the mantissa is
  97    shifted by _FP_WORKBITS but the implicit MSB is not inserted and
  98    other classification is not done.  */
  99 #define _FP_UNPACK_SEMIRAW(fs, wc, X)   _FP_FRAC_SLL_##wc (X, _FP_WORKBITS)
 100
 101 /* A semi-raw value has overflowed to infinity.  Adjust the mantissa
 102    and exponent appropriately.  */
 103 #define _FP_OVERFLOW_SEMIRAW(fs, wc, X)                 \
 104   do                                                    \
 105     {                                                   \
 106       if (FP_ROUNDMODE == FP_RND_NEAREST                \
 107           || (FP_ROUNDMODE == FP_RND_PINF && !X##_s)    \
 108           || (FP_ROUNDMODE == FP_RND_MINF && X##_s))    \
 109         {                                               \
 110           X##_e = _FP_EXPMAX_##fs;                      \
 111           _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc);     \
 112         }                                               \
 113       else                                              \
 114         {                                               \
 115           X##_e = _FP_EXPMAX_##fs - 1;                  \
 116           _FP_FRAC_SET_##wc (X, _FP_MAXFRAC_##wc);      \
 117         }                                               \
 118       FP_SET_EXCEPTION (FP_EX_INEXACT);                 \
 119       FP_SET_EXCEPTION (FP_EX_OVERFLOW);                \
 120     }                                                   \
 121   while (0)
 122
 123 /* Check for a semi-raw value being a signaling NaN and raise the
 124    invalid exception if so.  */
 125 #define _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X)     \
 126   do                                            \
 127     {                                           \
 128       if (X##_e == _FP_EXPMAX_##fs              \
 129           && !_FP_FRAC_ZEROP_##wc (X)           \
 130           && _FP_FRAC_SNANP_SEMIRAW (fs, X))    \
 131         FP_SET_EXCEPTION (FP_EX_INVALID);       \
 132     }                                           \
 133   while (0)
 134
 135 /* Choose a NaN result from an operation on two semi-raw NaN
 136    values.  */
 137 #define _FP_CHOOSENAN_SEMIRAW(fs, wc, R, X, Y, OP)                      \
 138   do                                                                    \
 139     {                                                                   \
 140       /* _FP_CHOOSENAN expects raw values, so shift as required.  */    \
 141       _FP_FRAC_SRL_##wc (X, _FP_WORKBITS);                              \
 142       _FP_FRAC_SRL_##wc (Y, _FP_WORKBITS);                              \
 143       _FP_CHOOSENAN (fs, wc, R, X, Y, OP);                              \
 144       _FP_FRAC_SLL_##wc (R, _FP_WORKBITS);                              \
 145     }                                                                   \
 146   while (0)
 147
 148 /* Make the fractional part a quiet NaN, preserving the payload
 149    if possible, otherwise make it the canonical quiet NaN and set
 150    the sign bit accordingly.  */
 151 #define _FP_SETQNAN(fs, wc, X)                                  \
 152   do                                                            \
 153     {                                                           \
 154       if (_FP_QNANNEGATEDP)                                     \
 155         {                                                       \
 156           _FP_FRAC_HIGH_RAW_##fs (X) &= _FP_QNANBIT_##fs - 1;   \
 157           if (_FP_FRAC_ZEROP_##wc (X))                          \
 158             {                                                   \
 159               X##_s = _FP_NANSIGN_##fs;                         \
 160               _FP_FRAC_SET_##wc (X, _FP_NANFRAC_##fs);          \
 161             }                                                   \
 162         }                                                       \
 163       else                                                      \
 164         _FP_FRAC_HIGH_RAW_##fs (X) |= _FP_QNANBIT_##fs;         \
 165     }                                                           \
 166   while (0)
 167 #define _FP_SETQNAN_SEMIRAW(fs, wc, X)                          \
 168   do                                                            \
 169     {                                                           \
 170       if (_FP_QNANNEGATEDP)                                     \
 171         {                                                       \
 172           _FP_FRAC_HIGH_##fs (X) &= _FP_QNANBIT_SH_##fs - 1;    \
 173           if (_FP_FRAC_ZEROP_##wc (X))                          \
 174             {                                                   \
 175               X##_s = _FP_NANSIGN_##fs;                         \
 176               _FP_FRAC_SET_##wc (X, _FP_NANFRAC_##fs);          \
 177               _FP_FRAC_SLL_##wc (X, _FP_WORKBITS);              \
 178             }                                                   \
 179         }                                                       \
 180       else                                                      \
 181         _FP_FRAC_HIGH_##fs (X) |= _FP_QNANBIT_SH_##fs;          \
 182     }                                                           \
 183   while (0)
 184
 185 /* Test whether a biased exponent is normal (not zero or maximum).  */
 186 #define _FP_EXP_NORMAL(fs, wc, X)       (((X##_e + 1) & _FP_EXPMAX_##fs) > 1)
 187
 188 /* Prepare to pack an fp value in semi-raw mode: the mantissa is
 189    rounded and shifted right, with the rounding possibly increasing
 190    the exponent (including changing a finite value to infinity).  */
 191 #define _FP_PACK_SEMIRAW(fs, wc, X)                             \
 192   do                                                            \
 193     {                                                           \
 194       int _FP_PACK_SEMIRAW_is_tiny                              \
 195         = X##_e == 0 && !_FP_FRAC_ZEROP_##wc (X);               \
 196       if (_FP_TININESS_AFTER_ROUNDING                           \
 197           && _FP_PACK_SEMIRAW_is_tiny)                          \
 198         {                                                       \
 199           FP_DECL_##fs (_FP_PACK_SEMIRAW_T);                    \
 200           _FP_FRAC_COPY_##wc (_FP_PACK_SEMIRAW_T, X);           \
 201           _FP_PACK_SEMIRAW_T##_s = X##_s;                       \
 202           _FP_PACK_SEMIRAW_T##_e = X##_e;                       \
 203           _FP_FRAC_SLL_##wc (_FP_PACK_SEMIRAW_T, 1);            \
 204           _FP_ROUND (wc, _FP_PACK_SEMIRAW_T);                   \
 205           if (_FP_FRAC_OVERP_##wc (fs, _FP_PACK_SEMIRAW_T))     \
 206             _FP_PACK_SEMIRAW_is_tiny = 0;                       \
 207         }                                                       \
 208       _FP_ROUND (wc, X);                                        \
 209       if (_FP_PACK_SEMIRAW_is_tiny)                             \
 210         {                                                       \
 211           if ((FP_CUR_EXCEPTIONS & FP_EX_INEXACT)               \
 212               || (FP_TRAPPING_EXCEPTIONS & FP_EX_UNDERFLOW))    \
 213             FP_SET_EXCEPTION (FP_EX_UNDERFLOW);                 \
 214         }                                                       \
 215       if (_FP_FRAC_HIGH_##fs (X)                                \
 216           & (_FP_OVERFLOW_##fs >> 1))                           \
 217         {                                                       \
 218           _FP_FRAC_HIGH_##fs (X) &= ~(_FP_OVERFLOW_##fs >> 1);  \
 219           X##_e++;                                              \
 220           if (X##_e == _FP_EXPMAX_##fs)                         \
 221             _FP_OVERFLOW_SEMIRAW (fs, wc, X);                   \
 222         }                                                       \
 223       _FP_FRAC_SRL_##wc (X, _FP_WORKBITS);                      \
 224       if (X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (X)) \
 225         {                                                       \
 226           if (!_FP_KEEPNANFRACP)                                \
 227             {                                                   \
 228               _FP_FRAC_SET_##wc (X, _FP_NANFRAC_##fs);          \
 229               X##_s = _FP_NANSIGN_##fs;                         \
 230             }                                                   \
 231           else                                                  \
 232             _FP_SETQNAN (fs, wc, X);                            \
 233         }                                                       \
 234     }                                                           \
 235   while (0)
 236
 237 /*
 238  * Before packing the bits back into the native fp result, take care
 239  * of such mundane things as rounding and overflow.  Also, for some
 240  * kinds of fp values, the original parts may not have been fully
 241  * extracted -- but that is ok, we can regenerate them now.
 242  */
 243
 244 #define _FP_PACK_CANONICAL(fs, wc, X)                                   \
 245   do                                                                    \
 246     {                                                                   \
 247       switch (X##_c)                                                    \
 248         {                                                               \
 249         case FP_CLS_NORMAL:                                             \
 250           X##_e += _FP_EXPBIAS_##fs;                                    \
 251           if (X##_e > 0)                                                \
 252             {                                                           \
 253               _FP_ROUND (wc, X);                                        \
 254               if (_FP_FRAC_OVERP_##wc (fs, X))                          \
 255                 {                                                       \
 256                   _FP_FRAC_CLEAR_OVERP_##wc (fs, X);                    \
 257                   X##_e++;                                              \
 258                 }                                                       \
 259               _FP_FRAC_SRL_##wc (X, _FP_WORKBITS);                      \
 260               if (X##_e >= _FP_EXPMAX_##fs)                             \
 261                 {                                                       \
 262                   /* overflow */                                        \
 263                   switch (FP_ROUNDMODE)                                 \
 264                     {                                                   \
 265                     case FP_RND_NEAREST:                                \
 266                       X##_c = FP_CLS_INF;                               \
 267                       break;                                            \
 268                     case FP_RND_PINF:                                   \
 269                       if (!X##_s)                                       \
 270                         X##_c = FP_CLS_INF;                             \
 271                       break;                                            \
 272                     case FP_RND_MINF:                                   \
 273                       if (X##_s)                                        \
 274                         X##_c = FP_CLS_INF;                             \
 275                       break;                                            \
 276                     }                                                   \
 277                   if (X##_c == FP_CLS_INF)                              \
 278                     {                                                   \
 279                       /* Overflow to infinity */                        \
 280                       X##_e = _FP_EXPMAX_##fs;                          \
 281                       _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc);         \
 282                     }                                                   \
 283                   else                                                  \
 284                     {                                                   \
 285                       /* Overflow to maximum normal */                  \
 286                       X##_e = _FP_EXPMAX_##fs - 1;                      \
 287                       _FP_FRAC_SET_##wc (X, _FP_MAXFRAC_##wc);          \
 288                     }                                                   \
 289                   FP_SET_EXCEPTION (FP_EX_OVERFLOW);                    \
 290                   FP_SET_EXCEPTION (FP_EX_INEXACT);                     \
 291                 }                                                       \
 292             }                                                           \
 293           else                                                          \
 294             {                                                           \
 295               /* we've got a denormalized number */                     \
 296               int _FP_PACK_CANONICAL_is_tiny = 1;                       \
 297               if (_FP_TININESS_AFTER_ROUNDING && X##_e == 0)            \
 298                 {                                                       \
 299                   FP_DECL_##fs (_FP_PACK_CANONICAL_T);                  \
 300                   _FP_FRAC_COPY_##wc (_FP_PACK_CANONICAL_T, X);         \
 301                   _FP_PACK_CANONICAL_T##_s = X##_s;                     \
 302                   _FP_PACK_CANONICAL_T##_e = X##_e;                     \
 303                   _FP_ROUND (wc, _FP_PACK_CANONICAL_T);                 \
 304                   if (_FP_FRAC_OVERP_##wc (fs, _FP_PACK_CANONICAL_T))   \
 305                     _FP_PACK_CANONICAL_is_tiny = 0;                     \
 306                 }                                                       \
 307               X##_e = -X##_e + 1;                                       \
 308               if (X##_e <= _FP_WFRACBITS_##fs)                          \
 309                 {                                                       \
 310                   _FP_FRAC_SRS_##wc (X, X##_e, _FP_WFRACBITS_##fs);     \
 311                   _FP_ROUND (wc, X);                                    \
 312                   if (_FP_FRAC_HIGH_##fs (X)                            \
 313                       & (_FP_OVERFLOW_##fs >> 1))                       \
 314                     {                                                   \
 315                       X##_e = 1;                                        \
 316                       _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc);         \
 317                       FP_SET_EXCEPTION (FP_EX_INEXACT);                 \
 318                     }                                                   \
 319                   else                                                  \
 320                     {                                                   \
 321                       X##_e = 0;                                        \
 322                       _FP_FRAC_SRL_##wc (X, _FP_WORKBITS);              \
 323                     }                                                   \
 324                   if (_FP_PACK_CANONICAL_is_tiny                        \
 325                       && ((FP_CUR_EXCEPTIONS & FP_EX_INEXACT)           \
 326                           || (FP_TRAPPING_EXCEPTIONS                    \
 327                               & FP_EX_UNDERFLOW)))                      \
 328                     FP_SET_EXCEPTION (FP_EX_UNDERFLOW);                 \
 329                 }                                                       \
 330               else                                                      \
 331                 {                                                       \
 332                   /* underflow to zero */                               \
 333                   X##_e = 0;                                            \
 334                   if (!_FP_FRAC_ZEROP_##wc (X))                         \
 335                     {                                                   \
 336                       _FP_FRAC_SET_##wc (X, _FP_MINFRAC_##wc);          \
 337                       _FP_ROUND (wc, X);                                \
 338                       _FP_FRAC_LOW_##wc (X) >>= (_FP_WORKBITS);         \
 339                     }                                                   \
 340                   FP_SET_EXCEPTION (FP_EX_UNDERFLOW);                   \
 341                 }                                                       \
 342             }                                                           \
 343           break;                                                        \
 344                                                                         \
 345         case FP_CLS_ZERO:                                               \
 346           X##_e = 0;                                                    \
 347           _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc);                     \
 348           break;                                                        \
 349                                                                         \
 350         case FP_CLS_INF:                                                \
 351           X##_e = _FP_EXPMAX_##fs;                                      \
 352           _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc);                     \
 353           break;                                                        \
 354                                                                         \
 355         case FP_CLS_NAN:                                                \
 356           X##_e = _FP_EXPMAX_##fs;                                      \
 357           if (!_FP_KEEPNANFRACP)                                        \
 358             {                                                           \
 359               _FP_FRAC_SET_##wc (X, _FP_NANFRAC_##fs);                  \
 360               X##_s = _FP_NANSIGN_##fs;                                 \
 361             }                                                           \
 362           else                                                          \
 363             _FP_SETQNAN (fs, wc, X);                                    \
 364           break;                                                        \
 365         }                                                               \
 366     }                                                                   \
 367   while (0)
 368
 369 /* This one accepts raw argument and not cooked,  returns
 370  * 1 if X is a signaling NaN.
 371  */
 372 #define _FP_ISSIGNAN(fs, wc, X)                 \
 373   ({                                            \
 374     int __ret = 0;                              \
 375     if (X##_e == _FP_EXPMAX_##fs)               \
 376       {                                         \
 377         if (!_FP_FRAC_ZEROP_##wc (X)            \
 378             && _FP_FRAC_SNANP (fs, X))          \
 379           __ret = 1;                            \
 380       }                                         \
 381     __ret;                                      \
 382   })
 383
 384
 385
 386
 387
 388 /* Addition on semi-raw values.  */
 389 #define _FP_ADD_INTERNAL(fs, wc, R, X, Y, OP)                           \
 390   do                                                                    \
 391     {                                                                   \
 392       if (X##_s == Y##_s)                                               \
 393         {                                                               \
 394           /* Addition.  */                                              \
 395           R##_s = X##_s;                                                \
 396           int ediff = X##_e - Y##_e;                                    \
 397           if (ediff > 0)                                                \
 398             {                                                           \
 399               R##_e = X##_e;                                            \
 400               if (Y##_e == 0)                                           \
 401                 {                                                       \
 402                   /* Y is zero or denormalized.  */                     \
 403                   if (_FP_FRAC_ZEROP_##wc (Y))                          \
 404                     {                                                   \
 405                       _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, X);             \
 406                       _FP_FRAC_COPY_##wc (R, X);                        \
 407                       goto add_done;                                    \
 408                     }                                                   \
 409                   else                                                  \
 410                     {                                                   \
 411                       FP_SET_EXCEPTION (FP_EX_DENORM);                  \
 412                       ediff--;                                          \
 413                       if (ediff == 0)                                   \
 414                         {                                               \
 415                           _FP_FRAC_ADD_##wc (R, X, Y);                  \
 416                           goto add3;                                    \
 417                         }                                               \
 418                       if (X##_e == _FP_EXPMAX_##fs)                     \
 419                         {                                               \
 420                           _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, X);         \
 421                           _FP_FRAC_COPY_##wc (R, X);                    \
 422                           goto add_done;                                \
 423                         }                                               \
 424                       goto add1;                                        \
 425                     }                                                   \
 426                 }                                                       \
 427               else if (X##_e == _FP_EXPMAX_##fs)                        \
 428                 {                                                       \
 429                   /* X is NaN or Inf, Y is normal.  */                  \
 430                   _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, X);                 \
 431                   _FP_FRAC_COPY_##wc (R, X);                            \
 432                   goto add_done;                                        \
 433                 }                                                       \
 434                                                                         \
 435               /* Insert implicit MSB of Y.  */                          \
 436               _FP_FRAC_HIGH_##fs (Y) |= _FP_IMPLBIT_SH_##fs;            \
 437                                                                         \
 438             add1:                                                       \
 439               /* Shift the mantissa of Y to the right EDIFF steps;      \
 440                  remember to account later for the implicit MSB of X.  */ \
 441               if (ediff <= _FP_WFRACBITS_##fs)                          \
 442                 _FP_FRAC_SRS_##wc (Y, ediff, _FP_WFRACBITS_##fs);       \
 443               else if (!_FP_FRAC_ZEROP_##wc (Y))                        \
 444                 _FP_FRAC_SET_##wc (Y, _FP_MINFRAC_##wc);                \
 445               _FP_FRAC_ADD_##wc (R, X, Y);                              \
 446             }                                                           \
 447           else if (ediff < 0)                                           \
 448             {                                                           \
 449               ediff = -ediff;                                           \
 450               R##_e = Y##_e;                                            \
 451               if (X##_e == 0)                                           \
 452                 {                                                       \
 453                   /* X is zero or denormalized.  */                     \
 454                   if (_FP_FRAC_ZEROP_##wc (X))                          \
 455                     {                                                   \
 456                       _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, Y);             \
 457                       _FP_FRAC_COPY_##wc (R, Y);                        \
 458                       goto add_done;                                    \
 459                     }                                                   \
 460                   else                                                  \
 461                     {                                                   \
 462                       FP_SET_EXCEPTION (FP_EX_DENORM);                  \
 463                       ediff--;                                          \
 464                       if (ediff == 0)                                   \
 465                         {                                               \
 466                           _FP_FRAC_ADD_##wc (R, Y, X);                  \
 467                           goto add3;                                    \
 468                         }                                               \
 469                       if (Y##_e == _FP_EXPMAX_##fs)                     \
 470                         {                                               \
 471                           _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, Y);         \
 472                           _FP_FRAC_COPY_##wc (R, Y);                    \
 473                           goto add_done;                                \
 474                         }                                               \
 475                       goto add2;                                        \
 476                     }                                                   \
 477                 }                                                       \
 478               else if (Y##_e == _FP_EXPMAX_##fs)                        \
 479                 {                                                       \
 480                   /* Y is NaN or Inf, X is normal.  */                  \
 481                   _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, Y);                 \
 482                   _FP_FRAC_COPY_##wc (R, Y);                            \
 483                   goto add_done;                                        \
 484                 }                                                       \
 485                                                                         \
 486               /* Insert implicit MSB of X.  */                          \
 487               _FP_FRAC_HIGH_##fs (X) |= _FP_IMPLBIT_SH_##fs;            \
 488                                                                         \
 489             add2:                                                       \
 490               /* Shift the mantissa of X to the right EDIFF steps;      \
 491                  remember to account later for the implicit MSB of Y.  */ \
 492               if (ediff <= _FP_WFRACBITS_##fs)                          \
 493                 _FP_FRAC_SRS_##wc (X, ediff, _FP_WFRACBITS_##fs);       \
 494               else if (!_FP_FRAC_ZEROP_##wc (X))                        \
 495                 _FP_FRAC_SET_##wc (X, _FP_MINFRAC_##wc);                \
 496               _FP_FRAC_ADD_##wc (R, Y, X);                              \
 497             }                                                           \
 498           else                                                          \
 499             {                                                           \
 500               /* ediff == 0.  */                                        \
 501               if (!_FP_EXP_NORMAL (fs, wc, X))                          \
 502                 {                                                       \
 503                   if (X##_e == 0)                                       \
 504                     {                                                   \
 505                       /* X and Y are zero or denormalized.  */          \
 506                       R##_e = 0;                                        \
 507                       if (_FP_FRAC_ZEROP_##wc (X))                      \
 508                         {                                               \
 509                           if (!_FP_FRAC_ZEROP_##wc (Y))                 \
 510                             FP_SET_EXCEPTION (FP_EX_DENORM);            \
 511                           _FP_FRAC_COPY_##wc (R, Y);                    \
 512                           goto add_done;                                \
 513                         }                                               \
 514                       else if (_FP_FRAC_ZEROP_##wc (Y))                 \
 515                         {                                               \
 516                           FP_SET_EXCEPTION (FP_EX_DENORM);              \
 517                           _FP_FRAC_COPY_##wc (R, X);                    \
 518                           goto add_done;                                \
 519                         }                                               \
 520                       else                                              \
 521                         {                                               \
 522                           FP_SET_EXCEPTION (FP_EX_DENORM);              \
 523                           _FP_FRAC_ADD_##wc (R, X, Y);                  \
 524                           if (_FP_FRAC_HIGH_##fs (R) & _FP_IMPLBIT_SH_##fs) \
 525                             {                                           \
 526                               /* Normalized result.  */                 \
 527                               _FP_FRAC_HIGH_##fs (R)                    \
 528                                 &= ~(_FP_W_TYPE) _FP_IMPLBIT_SH_##fs;   \
 529                               R##_e = 1;                                \
 530                             }                                           \
 531                           goto add_done;                                \
 532                         }                                               \
 533                     }                                                   \
 534                   else                                                  \
 535                     {                                                   \
 536                       /* X and Y are NaN or Inf.  */                    \
 537                       _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, X);             \
 538                       _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, Y);             \
 539                       R##_e = _FP_EXPMAX_##fs;                          \
 540                       if (_FP_FRAC_ZEROP_##wc (X))                      \
 541                         _FP_FRAC_COPY_##wc (R, Y);                      \
 542                       else if (_FP_FRAC_ZEROP_##wc (Y))                 \
 543                         _FP_FRAC_COPY_##wc (R, X);                      \
 544                       else                                              \
 545                         _FP_CHOOSENAN_SEMIRAW (fs, wc, R, X, Y, OP);    \
 546                       goto add_done;                                    \
 547                     }                                                   \
 548                 }                                                       \
 549               /* The exponents of X and Y, both normal, are equal.  The \
 550                  implicit MSBs will always add to increase the          \
 551                  exponent.  */                                          \
 552               _FP_FRAC_ADD_##wc (R, X, Y);                              \
 553               R##_e = X##_e + 1;                                        \
 554               _FP_FRAC_SRS_##wc (R, 1, _FP_WFRACBITS_##fs);             \
 555               if (R##_e == _FP_EXPMAX_##fs)                             \
 556                 /* Overflow to infinity (depending on rounding mode).  */ \
 557                 _FP_OVERFLOW_SEMIRAW (fs, wc, R);                       \
 558               goto add_done;                                            \
 559             }                                                           \
 560         add3:                                                           \
 561           if (_FP_FRAC_HIGH_##fs (R) & _FP_IMPLBIT_SH_##fs)             \
 562             {                                                           \
 563               /* Overflow.  */                                          \
 564               _FP_FRAC_HIGH_##fs (R) &= ~(_FP_W_TYPE) _FP_IMPLBIT_SH_##fs; \
 565               R##_e++;                                                  \
 566               _FP_FRAC_SRS_##wc (R, 1, _FP_WFRACBITS_##fs);             \
 567               if (R##_e == _FP_EXPMAX_##fs)                             \
 568                 /* Overflow to infinity (depending on rounding mode).  */ \
 569                 _FP_OVERFLOW_SEMIRAW (fs, wc, R);                       \
 570             }                                                           \
 571         add_done: ;                                                     \
 572         }                                                               \
 573       else                                                              \
 574         {                                                               \
 575           /* Subtraction.  */                                           \
 576           int ediff = X##_e - Y##_e;                                    \
 577           if (ediff > 0)                                                \
 578             {                                                           \
 579               R##_e = X##_e;                                            \
 580               R##_s = X##_s;                                            \
 581               if (Y##_e == 0)                                           \
 582                 {                                                       \
 583                   /* Y is zero or denormalized.  */                     \
 584                   if (_FP_FRAC_ZEROP_##wc (Y))                          \
 585                     {                                                   \
 586                       _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, X);             \
 587                       _FP_FRAC_COPY_##wc (R, X);                        \
 588                       goto sub_done;                                    \
 589                     }                                                   \
 590                   else                                                  \
 591                     {                                                   \
 592                       FP_SET_EXCEPTION (FP_EX_DENORM);                  \
 593                       ediff--;                                          \
 594                       if (ediff == 0)                                   \
 595                         {                                               \
 596                           _FP_FRAC_SUB_##wc (R, X, Y);                  \
 597                           goto sub3;                                    \
 598                         }                                               \
 599                       if (X##_e == _FP_EXPMAX_##fs)                     \
 600                         {                                               \
 601                           _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, X);         \
 602                           _FP_FRAC_COPY_##wc (R, X);                    \
 603                           goto sub_done;                                \
 604                         }                                               \
 605                       goto sub1;                                        \
 606                     }                                                   \
 607                 }                                                       \
 608               else if (X##_e == _FP_EXPMAX_##fs)                        \
 609                 {                                                       \
 610                   /* X is NaN or Inf, Y is normal.  */                  \
 611                   _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, X);                 \
 612                   _FP_FRAC_COPY_##wc (R, X);                            \
 613                   goto sub_done;                                        \
 614                 }                                                       \
 615                                                                         \
 616               /* Insert implicit MSB of Y.  */                          \
 617               _FP_FRAC_HIGH_##fs (Y) |= _FP_IMPLBIT_SH_##fs;            \
 618                                                                         \
 619             sub1:                                                       \
 620               /* Shift the mantissa of Y to the right EDIFF steps;      \
 621                  remember to account later for the implicit MSB of X.  */ \
 622               if (ediff <= _FP_WFRACBITS_##fs)                          \
 623                 _FP_FRAC_SRS_##wc (Y, ediff, _FP_WFRACBITS_##fs);       \
 624               else if (!_FP_FRAC_ZEROP_##wc (Y))                        \
 625                 _FP_FRAC_SET_##wc (Y, _FP_MINFRAC_##wc);                \
 626               _FP_FRAC_SUB_##wc (R, X, Y);                              \
 627             }                                                           \
 628           else if (ediff < 0)                                           \
 629             {                                                           \
 630               ediff = -ediff;                                           \
 631               R##_e = Y##_e;                                            \
 632               R##_s = Y##_s;                                            \
 633               if (X##_e == 0)                                           \
 634                 {                                                       \
 635                   /* X is zero or denormalized.  */                     \
 636                   if (_FP_FRAC_ZEROP_##wc (X))                          \
 637                     {                                                   \
 638                       _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, Y);             \
 639                       _FP_FRAC_COPY_##wc (R, Y);                        \
 640                       goto sub_done;                                    \
 641                     }                                                   \
 642                   else                                                  \
 643                     {                                                   \
 644                       FP_SET_EXCEPTION (FP_EX_DENORM);                  \
 645                       ediff--;                                          \
 646                       if (ediff == 0)                                   \
 647                         {                                               \
 648                           _FP_FRAC_SUB_##wc (R, Y, X);                  \
 649                           goto sub3;                                    \
 650                         }                                               \
 651                       if (Y##_e == _FP_EXPMAX_##fs)                     \
 652                         {                                               \
 653                           _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, Y);         \
 654                           _FP_FRAC_COPY_##wc (R, Y);                    \
 655                           goto sub_done;                                \
 656                         }                                               \
 657                       goto sub2;                                        \
 658                     }                                                   \
 659                 }                                                       \
 660               else if (Y##_e == _FP_EXPMAX_##fs)                        \
 661                 {                                                       \
 662                   /* Y is NaN or Inf, X is normal.  */                  \
 663                   _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, Y);                 \
 664                   _FP_FRAC_COPY_##wc (R, Y);                            \
 665                   goto sub_done;                                        \
 666                 }                                                       \
 667                                                                         \
 668               /* Insert implicit MSB of X.  */                          \
 669               _FP_FRAC_HIGH_##fs (X) |= _FP_IMPLBIT_SH_##fs;            \
 670                                                                         \
 671             sub2:                                                       \
 672               /* Shift the mantissa of X to the right EDIFF steps;      \
 673                  remember to account later for the implicit MSB of Y.  */ \
 674               if (ediff <= _FP_WFRACBITS_##fs)                          \
 675                 _FP_FRAC_SRS_##wc (X, ediff, _FP_WFRACBITS_##fs);       \
 676               else if (!_FP_FRAC_ZEROP_##wc (X))                        \
 677                 _FP_FRAC_SET_##wc (X, _FP_MINFRAC_##wc);                \
 678               _FP_FRAC_SUB_##wc (R, Y, X);                              \
 679             }                                                           \
 680           else                                                          \
 681             {                                                           \
 682               /* ediff == 0.  */                                        \
 683               if (!_FP_EXP_NORMAL (fs, wc, X))                          \
 684                 {                                                       \
 685                   if (X##_e == 0)                                       \
 686                     {                                                   \
 687                       /* X and Y are zero or denormalized.  */          \
 688                       R##_e = 0;                                        \
 689                       if (_FP_FRAC_ZEROP_##wc (X))                      \
 690                         {                                               \
 691                           _FP_FRAC_COPY_##wc (R, Y);                    \
 692                           if (_FP_FRAC_ZEROP_##wc (Y))                  \
 693                             R##_s = (FP_ROUNDMODE == FP_RND_MINF);      \
 694                           else                                          \
 695                             {                                           \
 696                               FP_SET_EXCEPTION (FP_EX_DENORM);          \
 697                               R##_s = Y##_s;                            \
 698                             }                                           \
 699                           goto sub_done;                                \
 700                         }                                               \
 701                       else if (_FP_FRAC_ZEROP_##wc (Y))                 \
 702                         {                                               \
 703                           FP_SET_EXCEPTION (FP_EX_DENORM);              \
 704                           _FP_FRAC_COPY_##wc (R, X);                    \
 705                           R##_s = X##_s;                                \
 706                           goto sub_done;                                \
 707                         }                                               \
 708                       else                                              \
 709                         {                                               \
 710                           FP_SET_EXCEPTION (FP_EX_DENORM);              \
 711                           _FP_FRAC_SUB_##wc (R, X, Y);                  \
 712                           R##_s = X##_s;                                \
 713                           if (_FP_FRAC_HIGH_##fs (R) & _FP_IMPLBIT_SH_##fs) \
 714                             {                                           \
 715                               /* |X| < |Y|, negate result.  */          \
 716                               _FP_FRAC_SUB_##wc (R, Y, X);              \
 717                               R##_s = Y##_s;                            \
 718                             }                                           \
 719                           else if (_FP_FRAC_ZEROP_##wc (R))             \
 720                             R##_s = (FP_ROUNDMODE == FP_RND_MINF);      \
 721                           goto sub_done;                                \
 722                         }                                               \
 723                     }                                                   \
 724                   else                                                  \
 725                     {                                                   \
 726                       /* X and Y are NaN or Inf, of opposite signs.  */ \
 727                       _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, X);             \
 728                       _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, Y);             \
 729                       R##_e = _FP_EXPMAX_##fs;                          \
 730                       if (_FP_FRAC_ZEROP_##wc (X))                      \
 731                         {                                               \
 732                           if (_FP_FRAC_ZEROP_##wc (Y))                  \
 733                             {                                           \
 734                               /* Inf - Inf.  */                         \
 735                               R##_s = _FP_NANSIGN_##fs;                 \
 736                               _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs);  \
 737                               _FP_FRAC_SLL_##wc (R, _FP_WORKBITS);      \
 738                               FP_SET_EXCEPTION (FP_EX_INVALID);         \
 739                             }                                           \
 740                           else                                          \
 741                             {                                           \
 742                               /* Inf - NaN.  */                         \
 743                               R##_s = Y##_s;                            \
 744                               _FP_FRAC_COPY_##wc (R, Y);                \
 745                             }                                           \
 746                         }                                               \
 747                       else                                              \
 748                         {                                               \
 749                           if (_FP_FRAC_ZEROP_##wc (Y))                  \
 750                             {                                           \
 751                               /* NaN - Inf.  */                         \
 752                               R##_s = X##_s;                            \
 753                               _FP_FRAC_COPY_##wc (R, X);                \
 754                             }                                           \
 755                           else                                          \
 756                             {                                           \
 757                               /* NaN - NaN.  */                         \
 758                               _FP_CHOOSENAN_SEMIRAW (fs, wc, R, X, Y, OP); \
 759                             }                                           \
 760                         }                                               \
 761                       goto sub_done;                                    \
 762                     }                                                   \
 763                 }                                                       \
 764               /* The exponents of X and Y, both normal, are equal.  The \
 765                  implicit MSBs cancel.  */                              \
 766               R##_e = X##_e;                                            \
 767               _FP_FRAC_SUB_##wc (R, X, Y);                              \
 768               R##_s = X##_s;                                            \
 769               if (_FP_FRAC_HIGH_##fs (R) & _FP_IMPLBIT_SH_##fs)         \
 770                 {                                                       \
 771                   /* |X| < |Y|, negate result.  */                      \
 772                   _FP_FRAC_SUB_##wc (R, Y, X);                          \
 773                   R##_s = Y##_s;                                        \
 774                 }                                                       \
 775               else if (_FP_FRAC_ZEROP_##wc (R))                         \
 776                 {                                                       \
 777                   R##_e = 0;                                            \
 778                   R##_s = (FP_ROUNDMODE == FP_RND_MINF);                \
 779                   goto sub_done;                                        \
 780                 }                                                       \
 781               goto norm;                                                \
 782             }                                                           \
 783         sub3:                                                           \
 784           if (_FP_FRAC_HIGH_##fs (R) & _FP_IMPLBIT_SH_##fs)             \
 785             {                                                           \
 786               int diff;                                                 \
 787               /* Carry into most significant bit of larger one of X and Y, \
 788                  canceling it; renormalize.  */                         \
 789               _FP_FRAC_HIGH_##fs (R) &= _FP_IMPLBIT_SH_##fs - 1;        \
 790             norm:                                                       \
 791               _FP_FRAC_CLZ_##wc (diff, R);                              \
 792               diff -= _FP_WFRACXBITS_##fs;                              \
 793               _FP_FRAC_SLL_##wc (R, diff);                              \
 794               if (R##_e <= diff)                                        \
 795                 {                                                       \
 796                   /* R is denormalized.  */                             \
 797                   diff = diff - R##_e + 1;                              \
 798                   _FP_FRAC_SRS_##wc (R, diff, _FP_WFRACBITS_##fs);      \
 799                   R##_e = 0;                                            \
 800                 }                                                       \
 801               else                                                      \
 802                 {                                                       \
 803                   R##_e -= diff;                                        \
 804                   _FP_FRAC_HIGH_##fs (R) &= ~(_FP_W_TYPE) _FP_IMPLBIT_SH_##fs; \
 805                 }                                                       \
 806             }                                                           \
 807         sub_done: ;                                                     \
 808         }                                                               \
 809     }                                                                   \
 810   while (0)
 811
 812 #define _FP_ADD(fs, wc, R, X, Y) _FP_ADD_INTERNAL (fs, wc, R, X, Y, '+')
 813 #define _FP_SUB(fs, wc, R, X, Y)                                        \
 814   do                                                                    \
 815     {                                                                   \
 816       if (!(Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (Y)))      \
 817         Y##_s ^= 1;                                                     \
 818       _FP_ADD_INTERNAL (fs, wc, R, X, Y, '-');                          \
 819     }                                                                   \
 820   while (0)
 821
 822
 823 /*
 824  * Main negation routine.  The input value is raw.
 825  */
 826
 827 #define _FP_NEG(fs, wc, R, X)                   \
 828   do                                            \
 829     {                                           \
 830       _FP_FRAC_COPY_##wc (R, X);                \
 831       R##_e = X##_e;                            \
 832       R##_s = 1 ^ X##_s;                        \
 833     }                                           \
 834   while (0)
 835
 836
 837 /*
 838  * Main multiplication routine.  The input values should be cooked.
 839  */
 840
 841 #define _FP_MUL(fs, wc, R, X, Y)                                \
 842   do                                                            \
 843     {                                                           \
 844       R##_s = X##_s ^ Y##_s;                                    \
 845       R##_e = X##_e + Y##_e + 1;                                \
 846       switch (_FP_CLS_COMBINE (X##_c, Y##_c))                   \
 847         {                                                       \
 848         case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NORMAL):    \
 849           R##_c = FP_CLS_NORMAL;                                \
 850                                                                 \
 851           _FP_MUL_MEAT_##fs (R, X, Y);                          \
 852                                                                 \
 853           if (_FP_FRAC_OVERP_##wc (fs, R))                      \
 854             _FP_FRAC_SRS_##wc (R, 1, _FP_WFRACBITS_##fs);       \
 855           else                                                  \
 856             R##_e--;                                            \
 857           break;                                                \
 858                                                                 \
 859         case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NAN):          \
 860           _FP_CHOOSENAN (fs, wc, R, X, Y, '*');                 \
 861           break;                                                \
 862                                                                 \
 863         case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NORMAL):       \
 864         case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_INF):          \
 865         case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_ZERO):         \
 866           R##_s = X##_s;                                        \
 867                                                                 \
 868         case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_INF):          \
 869         case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NORMAL):       \
 870         case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NORMAL):      \
 871         case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_ZERO):        \
 872           _FP_FRAC_COPY_##wc (R, X);                            \
 873           R##_c = X##_c;                                        \
 874           break;                                                \
 875                                                                 \
 876         case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NAN):       \
 877         case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NAN):          \
 878         case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NAN):         \
 879           R##_s = Y##_s;                                        \
 880                                                                 \
 881         case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_INF):       \
 882         case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_ZERO):      \
 883           _FP_FRAC_COPY_##wc (R, Y);                            \
 884           R##_c = Y##_c;                                        \
 885           break;                                                \
 886                                                                 \
 887         case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_ZERO):         \
 888         case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_INF):         \
 889           R##_s = _FP_NANSIGN_##fs;                             \
 890           R##_c = FP_CLS_NAN;                                   \
 891           _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs);              \
 892           FP_SET_EXCEPTION (FP_EX_INVALID);                     \
 893           break;                                                \
 894                                                                 \
 895         default:                                                \
 896           abort ();                                             \
 897         }                                                       \
 898     }                                                           \
 899   while (0)
 900
 901
 902 /* Fused multiply-add.  The input values should be cooked.  */
 903
 904 #define _FP_FMA(fs, wc, dwc, R, X, Y, Z)                                \
 905   do                                                                    \
 906     {                                                                   \
 907       FP_DECL_##fs (T);                                                 \
 908       T##_s = X##_s ^ Y##_s;                                            \
 909       T##_e = X##_e + Y##_e + 1;                                        \
 910       switch (_FP_CLS_COMBINE (X##_c, Y##_c))                           \
 911         {                                                               \
 912         case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NORMAL):            \
 913           switch (Z##_c)                                                \
 914             {                                                           \
 915             case FP_CLS_INF:                                            \
 916             case FP_CLS_NAN:                                            \
 917               R##_s = Z##_s;                                            \
 918               _FP_FRAC_COPY_##wc (R, Z);                                \
 919               R##_c = Z##_c;                                            \
 920               break;                                                    \
 921                                                                         \
 922             case FP_CLS_ZERO:                                           \
 923               R##_c = FP_CLS_NORMAL;                                    \
 924               R##_s = T##_s;                                            \
 925               R##_e = T##_e;                                            \
 926                                                                         \
 927               _FP_MUL_MEAT_##fs (R, X, Y);                              \
 928                                                                         \
 929               if (_FP_FRAC_OVERP_##wc (fs, R))                          \
 930                 _FP_FRAC_SRS_##wc (R, 1, _FP_WFRACBITS_##fs);           \
 931               else                                                      \
 932                 R##_e--;                                                \
 933               break;                                                    \
 934                                                                         \
 935             case FP_CLS_NORMAL:;                                        \
 936               _FP_FRAC_DECL_##dwc (TD);                                 \
 937               _FP_FRAC_DECL_##dwc (ZD);                                 \
 938               _FP_FRAC_DECL_##dwc (RD);                                 \
 939               _FP_MUL_MEAT_DW_##fs (TD, X, Y);                          \
 940               R##_e = T##_e;                                            \
 941               int tsh = _FP_FRAC_HIGHBIT_DW_##dwc (fs, TD) == 0;        \
 942               T##_e -= tsh;                                             \
 943               int ediff = T##_e - Z##_e;                                \
 944               if (ediff >= 0)                                           \
 945                 {                                                       \
 946                   int shift = _FP_WFRACBITS_##fs - tsh - ediff;         \
 947                   if (shift <= -_FP_WFRACBITS_##fs)                     \
 948                     _FP_FRAC_SET_##dwc (ZD, _FP_MINFRAC_##dwc);         \
 949                   else                                                  \
 950                     {                                                   \
 951                       _FP_FRAC_COPY_##dwc##_##wc (ZD, Z);               \
 952                       if (shift < 0)                                    \
 953                         _FP_FRAC_SRS_##dwc (ZD, -shift,                 \
 954                                             _FP_WFRACBITS_DW_##fs);     \
 955                       else if (shift > 0)                               \
 956                         _FP_FRAC_SLL_##dwc (ZD, shift);                 \
 957                     }                                                   \
 958                   R##_s = T##_s;                                        \
 959                   if (T##_s == Z##_s)                                   \
 960                     _FP_FRAC_ADD_##dwc (RD, TD, ZD);                    \
 961                   else                                                  \
 962                     {                                                   \
 963                       _FP_FRAC_SUB_##dwc (RD, TD, ZD);                  \
 964                       if (_FP_FRAC_NEGP_##dwc (RD))                     \
 965                         {                                               \
 966                           R##_s = Z##_s;                                \
 967                           _FP_FRAC_SUB_##dwc (RD, ZD, TD);              \
 968                         }                                               \
 969                     }                                                   \
 970                 }                                                       \
 971               else                                                      \
 972                 {                                                       \
 973                   R##_e = Z##_e;                                        \
 974                   R##_s = Z##_s;                                        \
 975                   _FP_FRAC_COPY_##dwc##_##wc (ZD, Z);                   \
 976                   _FP_FRAC_SLL_##dwc (ZD, _FP_WFRACBITS_##fs);          \
 977                   int shift = -ediff - tsh;                             \
 978                   if (shift >= _FP_WFRACBITS_DW_##fs)                   \
 979                     _FP_FRAC_SET_##dwc (TD, _FP_MINFRAC_##dwc);         \
 980                   else if (shift > 0)                                   \
 981                     _FP_FRAC_SRS_##dwc (TD, shift,                      \
 982                                         _FP_WFRACBITS_DW_##fs);         \
 983                   if (Z##_s == T##_s)                                   \
 984                     _FP_FRAC_ADD_##dwc (RD, ZD, TD);                    \
 985                   else                                                  \
 986                     _FP_FRAC_SUB_##dwc (RD, ZD, TD);                    \
 987                 }                                                       \
 988               if (_FP_FRAC_ZEROP_##dwc (RD))                            \
 989                 {                                                       \
 990                   if (T##_s == Z##_s)                                   \
 991                     R##_s = Z##_s;                                      \
 992                   else                                                  \
 993                     R##_s = (FP_ROUNDMODE == FP_RND_MINF);              \
 994                   _FP_FRAC_SET_##wc (R, _FP_ZEROFRAC_##wc);             \
 995                   R##_c = FP_CLS_ZERO;                                  \
 996                 }                                                       \
 997               else                                                      \
 998                 {                                                       \
 999                   int rlz;                                              \
1000                   _FP_FRAC_CLZ_##dwc (rlz, RD);                         \
1001                   rlz -= _FP_WFRACXBITS_DW_##fs;                        \
1002                   R##_e -= rlz;                                         \
1003                   int shift = _FP_WFRACBITS_##fs - rlz;                 \
1004                   if (shift > 0)                                        \
1005                     _FP_FRAC_SRS_##dwc (RD, shift,                      \
1006                                         _FP_WFRACBITS_DW_##fs);         \
1007                   else if (shift < 0)                                   \
1008                     _FP_FRAC_SLL_##dwc (RD, -shift);                    \
1009                   _FP_FRAC_COPY_##wc##_##dwc (R, RD);                   \
1010                   R##_c = FP_CLS_NORMAL;                                \
1011                 }                                                       \
1012               break;                                                    \
1013             }                                                           \
1014           goto done_fma;                                                \
1015                                                                         \
1016         case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NAN):                  \
1017           _FP_CHOOSENAN (fs, wc, T, X, Y, '*');                         \
1018           break;                                                        \
1019                                                                         \
1020         case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NORMAL):               \
1021         case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_INF):                  \
1022         case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_ZERO):                 \
1023           T##_s = X##_s;                                                \
1024                                                                         \
1025         case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_INF):                  \
1026         case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NORMAL):               \
1027         case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NORMAL):              \
1028         case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_ZERO):                \
1029           _FP_FRAC_COPY_##wc (T, X);                                    \
1030           T##_c = X##_c;                                                \
1031           break;                                                        \
1032                                                                         \
1033         case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NAN):               \
1034         case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NAN):                  \
1035         case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NAN):                 \
1036           T##_s = Y##_s;                                                \
1037                                                                         \
1038         case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_INF):               \
1039         case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_ZERO):              \
1040           _FP_FRAC_COPY_##wc (T, Y);                                    \
1041           T##_c = Y##_c;                                                \
1042           break;                                                        \
1043                                                                         \
1044         case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_ZERO):                 \
1045         case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_INF):                 \
1046           T##_s = _FP_NANSIGN_##fs;                                     \
1047           T##_c = FP_CLS_NAN;                                           \
1048           _FP_FRAC_SET_##wc (T, _FP_NANFRAC_##fs);                      \
1049           FP_SET_EXCEPTION (FP_EX_INVALID);                             \
1050           break;                                                        \
1051                                                                         \
1052         default:                                                        \
1053           abort ();                                                     \
1054         }                                                               \
1055                                                                         \
1056       /* T = X * Y is zero, infinity or NaN.  */                        \
1057       switch (_FP_CLS_COMBINE (T##_c, Z##_c))                           \
1058         {                                                               \
1059         case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NAN):                  \
1060           _FP_CHOOSENAN (fs, wc, R, T, Z, '+');                         \
1061           break;                                                        \
1062                                                                         \
1063         case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NORMAL):               \
1064         case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_INF):                  \
1065         case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_ZERO):                 \
1066         case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NORMAL):               \
1067         case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_ZERO):                 \
1068           R##_s = T##_s;                                                \
1069           _FP_FRAC_COPY_##wc (R, T);                                    \
1070           R##_c = T##_c;                                                \
1071           break;                                                        \
1072                                                                         \
1073         case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NAN):                  \
1074         case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NAN):                 \
1075         case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NORMAL):              \
1076         case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_INF):                 \
1077           R##_s = Z##_s;                                                \
1078           _FP_FRAC_COPY_##wc (R, Z);                                    \
1079           R##_c = Z##_c;                                                \
1080           break;                                                        \
1081                                                                         \
1082         case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_INF):                  \
1083           if (T##_s == Z##_s)                                           \
1084             {                                                           \
1085               R##_s = Z##_s;                                            \
1086               _FP_FRAC_COPY_##wc (R, Z);                                \
1087               R##_c = Z##_c;                                            \
1088             }                                                           \
1089           else                                                          \
1090             {                                                           \
1091               R##_s = _FP_NANSIGN_##fs;                                 \
1092               R##_c = FP_CLS_NAN;                                       \
1093               _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs);                  \
1094               FP_SET_EXCEPTION (FP_EX_INVALID);                         \
1095             }                                                           \
1096           break;                                                        \
1097                                                                         \
1098         case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_ZERO):                \
1099           if (T##_s == Z##_s)                                           \
1100             R##_s = Z##_s;                                              \
1101           else                                                          \
1102             R##_s = (FP_ROUNDMODE == FP_RND_MINF);                      \
1103           _FP_FRAC_COPY_##wc (R, Z);                                    \
1104           R##_c = Z##_c;                                                \
1105           break;                                                        \
1106                                                                         \
1107         default:                                                        \
1108           abort ();                                                     \
1109         }                                                               \
1110     done_fma: ;                                                         \
1111     }                                                                   \
1112   while (0)
1113
1114
1115 /*
1116  * Main division routine.  The input values should be cooked.
1117  */
1118
1119 #define _FP_DIV(fs, wc, R, X, Y)                                \
1120   do                                                            \
1121     {                                                           \
1122       R##_s = X##_s ^ Y##_s;                                    \
1123       R##_e = X##_e - Y##_e;                                    \
1124       switch (_FP_CLS_COMBINE (X##_c, Y##_c))                   \
1125         {                                                       \
1126         case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NORMAL):    \
1127           R##_c = FP_CLS_NORMAL;                                \
1128                                                                 \
1129           _FP_DIV_MEAT_##fs (R, X, Y);                          \
1130           break;                                                \
1131                                                                 \
1132         case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NAN):          \
1133           _FP_CHOOSENAN (fs, wc, R, X, Y, '/');                 \
1134           break;                                                \
1135                                                                 \
1136         case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NORMAL):       \
1137         case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_INF):          \
1138         case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_ZERO):         \
1139           R##_s = X##_s;                                        \
1140           _FP_FRAC_COPY_##wc (R, X);                            \
1141           R##_c = X##_c;                                        \
1142           break;                                                \
1143                                                                 \
1144         case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NAN):       \
1145         case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NAN):          \
1146         case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NAN):         \
1147           R##_s = Y##_s;                                        \
1148           _FP_FRAC_COPY_##wc (R, Y);                            \
1149           R##_c = Y##_c;                                        \
1150           break;                                                \
1151                                                                 \
1152         case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_INF):       \
1153         case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_INF):         \
1154         case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NORMAL):      \
1155           R##_c = FP_CLS_ZERO;                                  \
1156           break;                                                \
1157                                                                 \
1158         case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_ZERO):      \
1159           FP_SET_EXCEPTION (FP_EX_DIVZERO);                     \
1160         case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_ZERO):         \
1161         case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NORMAL):       \
1162           R##_c = FP_CLS_INF;                                   \
1163           break;                                                \
1164                                                                 \
1165         case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_INF):          \
1166         case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_ZERO):        \
1167           R##_s = _FP_NANSIGN_##fs;                             \
1168           R##_c = FP_CLS_NAN;                                   \
1169           _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs);              \
1170           FP_SET_EXCEPTION (FP_EX_INVALID);                     \
1171           break;                                                \
1172                                                                 \
1173         default:                                                \
1174           abort ();                                             \
1175         }                                                       \
1176     }                                                           \
1177   while (0)
1178
1179
1180 /*
1181  * Main differential comparison routine.  The inputs should be raw not
1182  * cooked.  The return is -1,0,1 for normal values, 2 otherwise.
1183  */
1184
1185 #define _FP_CMP(fs, wc, ret, X, Y, un)                                  \
1186   do                                                                    \
1187     {                                                                   \
1188       /* NANs are unordered */                                          \
1189       if ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (X))        \
1190           || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (Y)))    \
1191         {                                                               \
1192           ret = un;                                                     \
1193         }                                                               \
1194       else                                                              \
1195         {                                                               \
1196           int __is_zero_x;                                              \
1197           int __is_zero_y;                                              \
1198                                                                         \
1199           __is_zero_x = (!X##_e && _FP_FRAC_ZEROP_##wc (X)) ? 1 : 0;    \
1200           __is_zero_y = (!Y##_e && _FP_FRAC_ZEROP_##wc (Y)) ? 1 : 0;    \
1201                                                                         \
1202           if (__is_zero_x && __is_zero_y)                               \
1203             ret = 0;                                                    \
1204           else if (__is_zero_x)                                         \
1205             ret = Y##_s ? 1 : -1;                                       \
1206           else if (__is_zero_y)                                         \
1207             ret = X##_s ? -1 : 1;                                       \
1208           else if (X##_s != Y##_s)                                      \
1209             ret = X##_s ? -1 : 1;                                       \
1210           else if (X##_e > Y##_e)                                       \
1211             ret = X##_s ? -1 : 1;                                       \
1212           else if (X##_e < Y##_e)                                       \
1213             ret = X##_s ? 1 : -1;                                       \
1214           else if (_FP_FRAC_GT_##wc (X, Y))                             \
1215             ret = X##_s ? -1 : 1;                                       \
1216           else if (_FP_FRAC_GT_##wc (Y, X))                             \
1217             ret = X##_s ? 1 : -1;                                       \
1218           else                                                          \
1219             ret = 0;                                                    \
1220         }                                                               \
1221     }                                                                   \
1222   while (0)
1223
1224
1225 /* Simplification for strict equality.  */
1226
1227 #define _FP_CMP_EQ(fs, wc, ret, X, Y)                                   \
1228   do                                                                    \
1229     {                                                                   \
1230       /* NANs are unordered */                                          \
1231       if ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (X))        \
1232           || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (Y)))    \
1233         {                                                               \
1234           ret = 1;                                                      \
1235         }                                                               \
1236       else                                                              \
1237         {                                                               \
1238           ret = !(X##_e == Y##_e                                        \
1239                   && _FP_FRAC_EQ_##wc (X, Y)                            \
1240                   && (X##_s == Y##_s || (!X##_e && _FP_FRAC_ZEROP_##wc (X)))); \
1241         }                                                               \
1242     }                                                                   \
1243   while (0)
1244
1245 /* Version to test unordered.  */
1246
1247 #define _FP_CMP_UNORD(fs, wc, ret, X, Y)                                \
1248   do                                                                    \
1249     {                                                                   \
1250       ret = ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (X))     \
1251              || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (Y))); \
1252     }                                                                   \
1253   while (0)
1254
1255 /*
1256  * Main square root routine.  The input value should be cooked.
1257  */
1258
1259 #define _FP_SQRT(fs, wc, R, X)                                  \
1260   do                                                            \
1261     {                                                           \
1262       _FP_FRAC_DECL_##wc (T);                                   \
1263       _FP_FRAC_DECL_##wc (S);                                   \
1264       _FP_W_TYPE q;                                             \
1265       switch (X##_c)                                            \
1266         {                                                       \
1267         case FP_CLS_NAN:                                        \
1268           _FP_FRAC_COPY_##wc (R, X);                            \
1269           R##_s = X##_s;                                        \
1270           R##_c = FP_CLS_NAN;                                   \
1271           break;                                                \
1272         case FP_CLS_INF:                                        \
1273           if (X##_s)                                            \
1274             {                                                   \
1275               R##_s = _FP_NANSIGN_##fs;                         \
1276               R##_c = FP_CLS_NAN; /* NAN */                     \
1277               _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs);          \
1278               FP_SET_EXCEPTION (FP_EX_INVALID);                 \
1279             }                                                   \
1280           else                                                  \
1281             {                                                   \
1282               R##_s = 0;                                        \
1283               R##_c = FP_CLS_INF; /* sqrt(+inf) = +inf */       \
1284             }                                                   \
1285           break;                                                \
1286         case FP_CLS_ZERO:                                       \
1287           R##_s = X##_s;                                        \
1288           R##_c = FP_CLS_ZERO; /* sqrt(+-0) = +-0 */            \
1289           break;                                                \
1290         case FP_CLS_NORMAL:                                     \
1291           R##_s = 0;                                            \
1292           if (X##_s)                                            \
1293             {                                                   \
1294               R##_c = FP_CLS_NAN; /* NAN */                     \
1295               R##_s = _FP_NANSIGN_##fs;                         \
1296               _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs);          \
1297               FP_SET_EXCEPTION (FP_EX_INVALID);                 \
1298               break;                                            \
1299             }                                                   \
1300           R##_c = FP_CLS_NORMAL;                                \
1301           if (X##_e & 1)                                        \
1302             _FP_FRAC_SLL_##wc (X, 1);                           \
1303           R##_e = X##_e >> 1;                                   \
1304           _FP_FRAC_SET_##wc (S, _FP_ZEROFRAC_##wc);             \
1305           _FP_FRAC_SET_##wc (R, _FP_ZEROFRAC_##wc);             \
1306           q = _FP_OVERFLOW_##fs >> 1;                           \
1307           _FP_SQRT_MEAT_##wc (R, S, T, X, q);                   \
1308         }                                                       \
1309     }                                                           \
1310   while (0)
1311
1312 /*
1313  * Convert from FP to integer.  Input is raw.
1314  */
1315
1316 /* RSIGNED can have following values:
1317  * 0:  the number is required to be 0..(2^rsize)-1, if not, NV is set plus
1318  *     the result is either 0 or (2^rsize)-1 depending on the sign in such
1319  *     case.
1320  * 1:  the number is required to be -(2^(rsize-1))..(2^(rsize-1))-1, if not,
1321  *     NV is set plus the result is either -(2^(rsize-1)) or (2^(rsize-1))-1
1322  *     depending on the sign in such case.
1323  * -1: the number is required to be -(2^(rsize-1))..(2^rsize)-1, if not, NV is
1324  *     set plus the result is either -(2^(rsize-1)) or (2^(rsize-1))-1
1325  *     depending on the sign in such case.
1326  */
1327 #define _FP_TO_INT(fs, wc, r, X, rsize, rsigned)                        \
1328   do                                                                    \
1329     {                                                                   \
1330       if (X##_e < _FP_EXPBIAS_##fs)                                     \
1331         {                                                               \
1332           r = 0;                                                        \
1333           if (X##_e == 0)                                               \
1334             {                                                           \
1335               if (!_FP_FRAC_ZEROP_##wc (X))                             \
1336                 {                                                       \
1337                   FP_SET_EXCEPTION (FP_EX_INEXACT);                     \
1338                   FP_SET_EXCEPTION (FP_EX_DENORM);                      \
1339                 }                                                       \
1340             }                                                           \
1341           else                                                          \
1342             FP_SET_EXCEPTION (FP_EX_INEXACT);                           \
1343         }                                                               \
1344       else if (X##_e >= _FP_EXPBIAS_##fs + rsize - (rsigned > 0 || X##_s) \
1345                || (!rsigned && X##_s))                                  \
1346         {                                                               \
1347           /* Overflow or converting to the most negative integer.  */   \
1348           if (rsigned)                                                  \
1349             {                                                           \
1350               r = 1;                                                    \
1351               r <<= rsize - 1;                                          \
1352               r -= 1 - X##_s;                                           \
1353             } else {                                                    \
1354             r = 0;                                                      \
1355             if (!X##_s)                                                 \
1356               r = ~r;                                                   \
1357           }                                                             \
1358                                                                         \
1359           if (rsigned && X##_s && X##_e == _FP_EXPBIAS_##fs + rsize - 1) \
1360             {                                                           \
1361               /* Possibly converting to most negative integer; check the \
1362                  mantissa.  */                                          \
1363               int inexact = 0;                                          \
1364               (void) ((_FP_FRACBITS_##fs > rsize)                       \
1365                       ? ({                                              \
1366                           _FP_FRAC_SRST_##wc (X, inexact,               \
1367                                               _FP_FRACBITS_##fs - rsize, \
1368                                               _FP_FRACBITS_##fs);       \
1369                           0;                                            \
1370                         })                                              \
1371                       : 0);                                             \
1372               if (!_FP_FRAC_ZEROP_##wc (X))                             \
1373                 FP_SET_EXCEPTION (FP_EX_INVALID);                       \
1374               else if (inexact)                                         \
1375                 FP_SET_EXCEPTION (FP_EX_INEXACT);                       \
1376             }                                                           \
1377           else                                                          \
1378             FP_SET_EXCEPTION (FP_EX_INVALID);                           \
1379         }                                                               \
1380       else                                                              \
1381         {                                                               \
1382           _FP_FRAC_HIGH_RAW_##fs (X) |= _FP_IMPLBIT_##fs;               \
1383           if (X##_e >= _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs - 1)        \
1384             {                                                           \
1385               _FP_FRAC_ASSEMBLE_##wc (r, X, rsize);                     \
1386               r <<= X##_e - _FP_EXPBIAS_##fs - _FP_FRACBITS_##fs + 1;   \
1387             }                                                           \
1388           else                                                          \
1389             {                                                           \
1390               int inexact;                                              \
1391               _FP_FRAC_SRST_##wc (X, inexact,                           \
1392                                   (_FP_FRACBITS_##fs + _FP_EXPBIAS_##fs - 1 \
1393                                    - X##_e),                            \
1394                                   _FP_FRACBITS_##fs);                   \
1395               if (inexact)                                              \
1396                 FP_SET_EXCEPTION (FP_EX_INEXACT);                       \
1397               _FP_FRAC_ASSEMBLE_##wc (r, X, rsize);                     \
1398             }                                                           \
1399           if (rsigned && X##_s)                                         \
1400             r = -r;                                                     \
1401         }                                                               \
1402     }                                                                   \
1403   while (0)
1404
1405 /* Convert integer to fp.  Output is raw.  RTYPE is unsigned even if
1406    input is signed.  */
1407 #define _FP_FROM_INT(fs, wc, X, r, rsize, rtype)                        \
1408   do                                                                    \
1409     {                                                                   \
1410       if (r)                                                            \
1411         {                                                               \
1412           rtype ur_;                                                    \
1413                                                                         \
1414           if ((X##_s = (r < 0)))                                        \
1415             r = -(rtype) r;                                             \
1416                                                                         \
1417           ur_ = (rtype) r;                                              \
1418           (void) ((rsize <= _FP_W_TYPE_SIZE)                            \
1419                   ? ({                                                  \
1420                       int lz_;                                          \
1421                       __FP_CLZ (lz_, (_FP_W_TYPE) ur_);                 \
1422                       X##_e = _FP_EXPBIAS_##fs + _FP_W_TYPE_SIZE - 1 - lz_; \
1423                     })                                                  \
1424                   : ((rsize <= 2 * _FP_W_TYPE_SIZE)                     \
1425                      ? ({                                               \
1426                          int lz_;                                       \
1427                          __FP_CLZ_2 (lz_,                               \
1428                                      (_FP_W_TYPE) (ur_ >> _FP_W_TYPE_SIZE), \
1429                                      (_FP_W_TYPE) ur_);                 \
1430                          X##_e = (_FP_EXPBIAS_##fs + 2 * _FP_W_TYPE_SIZE - 1 \
1431                                   - lz_);                               \
1432                        })                                               \
1433                      : (abort (), 0)));                                 \
1434                                                                         \
1435           if (rsize - 1 + _FP_EXPBIAS_##fs >= _FP_EXPMAX_##fs           \
1436               && X##_e >= _FP_EXPMAX_##fs)                              \
1437             {                                                           \
1438               /* Exponent too big; overflow to infinity.  (May also     \
1439                  happen after rounding below.)  */                      \
1440               _FP_OVERFLOW_SEMIRAW (fs, wc, X);                         \
1441               goto pack_semiraw;                                        \
1442             }                                                           \
1443                                                                         \
1444           if (rsize <= _FP_FRACBITS_##fs                                \
1445               || X##_e < _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs)          \
1446             {                                                           \
1447               /* Exactly representable; shift left.  */                 \
1448               _FP_FRAC_DISASSEMBLE_##wc (X, ur_, rsize);                \
1449               if (_FP_EXPBIAS_##fs + _FP_FRACBITS_##fs - 1 - X##_e > 0) \
1450                 _FP_FRAC_SLL_##wc (X, (_FP_EXPBIAS_##fs                 \
1451                                        + _FP_FRACBITS_##fs - 1 - X##_e)); \
1452             }                                                           \
1453           else                                                          \
1454             {                                                           \
1455               /* More bits in integer than in floating type; need to    \
1456                  round.  */                                             \
1457               if (_FP_EXPBIAS_##fs + _FP_WFRACBITS_##fs - 1 < X##_e)    \
1458                 ur_ = ((ur_ >> (X##_e - _FP_EXPBIAS_##fs                \
1459                                 - _FP_WFRACBITS_##fs + 1))              \
1460                        | ((ur_ << (rsize - (X##_e - _FP_EXPBIAS_##fs    \
1461                                             - _FP_WFRACBITS_##fs + 1))) \
1462                           != 0));                                       \
1463               _FP_FRAC_DISASSEMBLE_##wc (X, ur_, rsize);                \
1464               if ((_FP_EXPBIAS_##fs + _FP_WFRACBITS_##fs - 1 - X##_e) > 0) \
1465                 _FP_FRAC_SLL_##wc (X, (_FP_EXPBIAS_##fs                 \
1466                                        + _FP_WFRACBITS_##fs - 1 - X##_e)); \
1467               _FP_FRAC_HIGH_##fs (X) &= ~(_FP_W_TYPE) _FP_IMPLBIT_SH_##fs; \
1468             pack_semiraw:                                               \
1469               _FP_PACK_SEMIRAW (fs, wc, X);                             \
1470             }                                                           \
1471         }                                                               \
1472       else                                                              \
1473         {                                                               \
1474           X##_s = 0;                                                    \
1475           X##_e = 0;                                                    \
1476           _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc);                     \
1477         }                                                               \
1478     }                                                                   \
1479   while (0)
1480
1481
1482 /* Extend from a narrower floating-point format to a wider one.  Input
1483    and output are raw.  */
1484 #define FP_EXTEND(dfs, sfs, dwc, swc, D, S)                             \
1485   do                                                                    \
1486     {                                                                   \
1487       if (_FP_FRACBITS_##dfs < _FP_FRACBITS_##sfs                       \
1488           || (_FP_EXPMAX_##dfs - _FP_EXPBIAS_##dfs                      \
1489               < _FP_EXPMAX_##sfs - _FP_EXPBIAS_##sfs)                   \
1490           || (_FP_EXPBIAS_##dfs < _FP_EXPBIAS_##sfs + _FP_FRACBITS_##sfs - 1 \
1491               && _FP_EXPBIAS_##dfs != _FP_EXPBIAS_##sfs))               \
1492         abort ();                                                       \
1493       D##_s = S##_s;                                                    \
1494       _FP_FRAC_COPY_##dwc##_##swc (D, S);                               \
1495       if (_FP_EXP_NORMAL (sfs, swc, S))                                 \
1496         {                                                               \
1497           D##_e = S##_e + _FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs;        \
1498           _FP_FRAC_SLL_##dwc (D, (_FP_FRACBITS_##dfs - _FP_FRACBITS_##sfs)); \
1499         }                                                               \
1500       else                                                              \
1501         {                                                               \
1502           if (S##_e == 0)                                               \
1503             {                                                           \
1504               if (_FP_FRAC_ZEROP_##swc (S))                             \
1505                 D##_e = 0;                                              \
1506               else if (_FP_EXPBIAS_##dfs                                \
1507                        < _FP_EXPBIAS_##sfs + _FP_FRACBITS_##sfs - 1)    \
1508                 {                                                       \
1509                   FP_SET_EXCEPTION (FP_EX_DENORM);                      \
1510                   _FP_FRAC_SLL_##dwc (D, (_FP_FRACBITS_##dfs            \
1511                                           - _FP_FRACBITS_##sfs));       \
1512                   D##_e = 0;                                            \
1513                 }                                                       \
1514               else                                                      \
1515                 {                                                       \
1516                   int _lz;                                              \
1517                   FP_SET_EXCEPTION (FP_EX_DENORM);                      \
1518                   _FP_FRAC_CLZ_##swc (_lz, S);                          \
1519                   _FP_FRAC_SLL_##dwc (D,                                \
1520                                       _lz + _FP_FRACBITS_##dfs          \
1521                                       - _FP_FRACTBITS_##sfs);           \
1522                   D##_e = (_FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs + 1    \
1523                            + _FP_FRACXBITS_##sfs - _lz);                \
1524                 }                                                       \
1525             }                                                           \
1526           else                                                          \
1527             {                                                           \
1528               D##_e = _FP_EXPMAX_##dfs;                                 \
1529               if (!_FP_FRAC_ZEROP_##swc (S))                            \
1530                 {                                                       \
1531                   if (_FP_FRAC_SNANP (sfs, S))                          \
1532                     FP_SET_EXCEPTION (FP_EX_INVALID);                   \
1533                   _FP_FRAC_SLL_##dwc (D, (_FP_FRACBITS_##dfs            \
1534                                           - _FP_FRACBITS_##sfs));       \
1535                   _FP_SETQNAN (dfs, dwc, D);                            \
1536                 }                                                       \
1537             }                                                           \
1538         }                                                               \
1539     }                                                                   \
1540   while (0)
1541
1542 /* Truncate from a wider floating-point format to a narrower one.
1543    Input and output are semi-raw.  */
1544 #define FP_TRUNC(dfs, sfs, dwc, swc, D, S)                              \
1545   do                                                                    \
1546     {                                                                   \
1547       if (_FP_FRACBITS_##sfs < _FP_FRACBITS_##dfs                       \
1548           || (_FP_EXPBIAS_##sfs < _FP_EXPBIAS_##dfs + _FP_FRACBITS_##dfs - 1 \
1549               && _FP_EXPBIAS_##sfs != _FP_EXPBIAS_##dfs))               \
1550         abort ();                                                       \
1551       D##_s = S##_s;                                                    \
1552       if (_FP_EXP_NORMAL (sfs, swc, S))                                 \
1553         {                                                               \
1554           D##_e = S##_e + _FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs;        \
1555           if (D##_e >= _FP_EXPMAX_##dfs)                                \
1556             _FP_OVERFLOW_SEMIRAW (dfs, dwc, D);                         \
1557           else                                                          \
1558             {                                                           \
1559               if (D##_e <= 0)                                           \
1560                 {                                                       \
1561                   if (D##_e < 1 - _FP_FRACBITS_##dfs)                   \
1562                     {                                                   \
1563                       _FP_FRAC_SET_##swc (S, _FP_ZEROFRAC_##swc);       \
1564                       _FP_FRAC_LOW_##swc (S) |= 1;                      \
1565                     }                                                   \
1566                   else                                                  \
1567                     {                                                   \
1568                       _FP_FRAC_HIGH_##sfs (S) |= _FP_IMPLBIT_SH_##sfs;  \
1569                       _FP_FRAC_SRS_##swc (S, (_FP_WFRACBITS_##sfs       \
1570                                               - _FP_WFRACBITS_##dfs     \
1571                                               + 1 - D##_e),             \
1572                                           _FP_WFRACBITS_##sfs);         \
1573                     }                                                   \
1574                   D##_e = 0;                                            \
1575                 }                                                       \
1576               else                                                      \
1577                 _FP_FRAC_SRS_##swc (S, (_FP_WFRACBITS_##sfs             \
1578                                         - _FP_WFRACBITS_##dfs),         \
1579                                     _FP_WFRACBITS_##sfs);               \
1580               _FP_FRAC_COPY_##dwc##_##swc (D, S);                       \
1581             }                                                           \
1582         }                                                               \
1583       else                                                              \
1584         {                                                               \
1585           if (S##_e == 0)                                               \
1586             {                                                           \
1587               D##_e = 0;                                                \
1588               if (_FP_FRAC_ZEROP_##swc (S))                             \
1589                 _FP_FRAC_SET_##dwc (D, _FP_ZEROFRAC_##dwc);             \
1590               else                                                      \
1591                 {                                                       \
1592                   FP_SET_EXCEPTION (FP_EX_DENORM);                      \
1593                   if (_FP_EXPBIAS_##sfs                                 \
1594                       < _FP_EXPBIAS_##dfs + _FP_FRACBITS_##dfs - 1)     \
1595                     {                                                   \
1596                       _FP_FRAC_SRS_##swc (S, (_FP_WFRACBITS_##sfs       \
1597                                               - _FP_WFRACBITS_##dfs),   \
1598                                           _FP_WFRACBITS_##sfs);         \
1599                       _FP_FRAC_COPY_##dwc##_##swc (D, S);               \
1600                     }                                                   \
1601                   else                                                  \
1602                     {                                                   \
1603                       _FP_FRAC_SET_##dwc (D, _FP_ZEROFRAC_##dwc);       \
1604                       _FP_FRAC_LOW_##dwc (D) |= 1;                      \
1605                     }                                                   \
1606                 }                                                       \
1607             }                                                           \
1608           else                                                          \
1609             {                                                           \
1610               D##_e = _FP_EXPMAX_##dfs;                                 \
1611               if (_FP_FRAC_ZEROP_##swc (S))                             \
1612                 _FP_FRAC_SET_##dwc (D, _FP_ZEROFRAC_##dwc);             \
1613               else                                                      \
1614                 {                                                       \
1615                   _FP_CHECK_SIGNAN_SEMIRAW (sfs, swc, S);               \
1616                   _FP_FRAC_SRL_##swc (S, (_FP_WFRACBITS_##sfs           \
1617                                           - _FP_WFRACBITS_##dfs));      \
1618                   _FP_FRAC_COPY_##dwc##_##swc (D, S);                   \
1619                   /* Semi-raw NaN must have all workbits cleared.  */   \
1620                   _FP_FRAC_LOW_##dwc (D)                                \
1621                     &= ~(_FP_W_TYPE) ((1 << _FP_WORKBITS) - 1);         \
1622                   _FP_SETQNAN_SEMIRAW (dfs, dwc, D);                    \
1623                 }                                                       \
1624             }                                                           \
1625         }                                                               \
1626     }                                                                   \
1627   while (0)
1628
1629 /*
1630  * Helper primitives.
1631  */
1632
1633 /* Count leading zeros in a word.  */
1634
1635 #ifndef __FP_CLZ
1636 /* GCC 3.4 and later provide the builtins for us.  */
1637 # define __FP_CLZ(r, x)                                                 \
1638   do                                                                    \
1639     {                                                                   \
1640       if (sizeof (_FP_W_TYPE) == sizeof (unsigned int))                 \
1641         r = __builtin_clz (x);                                          \
1642       else if (sizeof (_FP_W_TYPE) == sizeof (unsigned long))           \
1643         r = __builtin_clzl (x);                                         \
1644       else if (sizeof (_FP_W_TYPE) == sizeof (unsigned long long))      \
1645         r = __builtin_clzll (x);                                        \
1646       else                                                              \
1647         abort ();                                                       \
1648     }                                                                   \
1649   while (0)
1650 #endif /* ndef __FP_CLZ */
1651
1652 #define _FP_DIV_HELP_imm(q, r, n, d)            \
1653   do                                            \
1654     {                                           \
1655       q = n / d, r = n % d;                     \
1656     }                                           \
1657   while (0)
1658
1659
1660 /* A restoring bit-by-bit division primitive.  */
1661
1662 #define _FP_DIV_MEAT_N_loop(fs, wc, R, X, Y)                            \
1663   do                                                                    \
1664     {                                                                   \
1665       int count = _FP_WFRACBITS_##fs;                                   \
1666       _FP_FRAC_DECL_##wc (u);                                           \
1667       _FP_FRAC_DECL_##wc (v);                                           \
1668       _FP_FRAC_COPY_##wc (u, X);                                        \
1669       _FP_FRAC_COPY_##wc (v, Y);                                        \
1670       _FP_FRAC_SET_##wc (R, _FP_ZEROFRAC_##wc);                         \
1671       /* Normalize U and V.  */                                         \
1672       _FP_FRAC_SLL_##wc (u, _FP_WFRACXBITS_##fs);                       \
1673       _FP_FRAC_SLL_##wc (v, _FP_WFRACXBITS_##fs);                       \
1674       /* First round.  Since the operands are normalized, either the    \
1675          first or second bit will be set in the fraction.  Produce a    \
1676          normalized result by checking which and adjusting the loop     \
1677          count and exponent accordingly.  */                            \
1678       if (_FP_FRAC_GE_1 (u, v))                                         \
1679         {                                                               \
1680           _FP_FRAC_SUB_##wc (u, u, v);                                  \
1681           _FP_FRAC_LOW_##wc (R) |= 1;                                   \
1682           count--;                                                      \
1683         }                                                               \
1684       else                                                              \
1685         R##_e--;                                                        \
1686       /* Subsequent rounds.  */                                         \
1687       do                                                                \
1688         {                                                               \
1689           int msb = (_FP_WS_TYPE) _FP_FRAC_HIGH_##wc (u) < 0;           \
1690           _FP_FRAC_SLL_##wc (u, 1);                                     \
1691           _FP_FRAC_SLL_##wc (R, 1);                                     \
1692           if (msb || _FP_FRAC_GE_1 (u, v))                              \
1693             {                                                           \
1694               _FP_FRAC_SUB_##wc (u, u, v);                              \
1695               _FP_FRAC_LOW_##wc (R) |= 1;                               \
1696             }                                                           \
1697         }                                                               \
1698       while (--count > 0);                                              \
1699       /* If there's anything left in U, the result is inexact.  */      \
1700       _FP_FRAC_LOW_##wc (R) |= !_FP_FRAC_ZEROP_##wc (u);                \
1701     }                                                                   \
1702   while (0)
1703
1704 #define _FP_DIV_MEAT_1_loop(fs, R, X, Y)  _FP_DIV_MEAT_N_loop (fs, 1, R, X, Y)
1705 #define _FP_DIV_MEAT_2_loop(fs, R, X, Y)  _FP_DIV_MEAT_N_loop (fs, 2, R, X, Y)
1706 #define _FP_DIV_MEAT_4_loop(fs, R, X, Y)  _FP_DIV_MEAT_N_loop (fs, 4, R, X, Y)