soft-fp/op-common.h

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