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