soft-fp/op-common.h

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