contrib/gcc/config/soft-fp/op-common.h

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