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