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