Back to index

glibc  2.9
op-common.h
Go to the documentation of this file.
00001 /* Software floating-point emulation. Common operations.
00002    Copyright (C) 1997,1998,1999,2006,2007 Free Software Foundation, Inc.
00003    This file is part of the GNU C Library.
00004    Contributed by Richard Henderson (rth@cygnus.com),
00005                 Jakub Jelinek (jj@ultra.linux.cz),
00006                 David S. Miller (davem@redhat.com) and
00007                 Peter Maydell (pmaydell@chiark.greenend.org.uk).
00008 
00009    The GNU C Library is free software; you can redistribute it and/or
00010    modify it under the terms of the GNU Lesser General Public
00011    License as published by the Free Software Foundation; either
00012    version 2.1 of the License, or (at your option) any later version.
00013 
00014    In addition to the permissions in the GNU Lesser General Public
00015    License, the Free Software Foundation gives you unlimited
00016    permission to link the compiled version of this file into
00017    combinations with other programs, and to distribute those
00018    combinations without any restriction coming from the use of this
00019    file.  (The Lesser General Public License restrictions do apply in
00020    other respects; for example, they cover modification of the file,
00021    and distribution when not linked into a combine executable.)
00022 
00023    The GNU C Library is distributed in the hope that it will be useful,
00024    but WITHOUT ANY WARRANTY; without even the implied warranty of
00025    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00026    Lesser General Public License for more details.
00027 
00028    You should have received a copy of the GNU Lesser General Public
00029    License along with the GNU C Library; if not, write to the Free
00030    Software Foundation, 51 Franklin Street, Fifth Floor, Boston,
00031    MA 02110-1301, USA.  */
00032 
00033 #define _FP_DECL(wc, X)                                        \
00034   _FP_I_TYPE X##_c __attribute__((unused)), X##_s, X##_e;      \
00035   _FP_FRAC_DECL_##wc(X)
00036 
00037 /*
00038  * Finish truely unpacking a native fp value by classifying the kind
00039  * of fp value and normalizing both the exponent and the fraction.
00040  */
00041 
00042 #define _FP_UNPACK_CANONICAL(fs, wc, X)                               \
00043 do {                                                           \
00044   switch (X##_e)                                               \
00045   {                                                            \
00046   default:                                                     \
00047     _FP_FRAC_HIGH_RAW_##fs(X) |= _FP_IMPLBIT_##fs;                    \
00048     _FP_FRAC_SLL_##wc(X, _FP_WORKBITS);                               \
00049     X##_e -= _FP_EXPBIAS_##fs;                                        \
00050     X##_c = FP_CLS_NORMAL;                                     \
00051     break;                                                     \
00052                                                                \
00053   case 0:                                                      \
00054     if (_FP_FRAC_ZEROP_##wc(X))                                       \
00055       X##_c = FP_CLS_ZERO;                                     \
00056     else                                                       \
00057       {                                                               \
00058        /* a denormalized number */                             \
00059        _FP_I_TYPE _shift;                                      \
00060        _FP_FRAC_CLZ_##wc(_shift, X);                                  \
00061        _shift -= _FP_FRACXBITS_##fs;                                  \
00062        _FP_FRAC_SLL_##wc(X, (_shift+_FP_WORKBITS));                   \
00063        X##_e -= _FP_EXPBIAS_##fs - 1 + _shift;                        \
00064        X##_c = FP_CLS_NORMAL;                                         \
00065        FP_SET_EXCEPTION(FP_EX_DENORM);                                \
00066       }                                                               \
00067     break;                                                     \
00068                                                                \
00069   case _FP_EXPMAX_##fs:                                               \
00070     if (_FP_FRAC_ZEROP_##wc(X))                                       \
00071       X##_c = FP_CLS_INF;                                      \
00072     else                                                       \
00073       {                                                               \
00074        X##_c = FP_CLS_NAN;                                     \
00075        /* Check for signaling NaN */                                  \
00076        if (!(_FP_FRAC_HIGH_RAW_##fs(X) & _FP_QNANBIT_##fs))           \
00077          FP_SET_EXCEPTION(FP_EX_INVALID);                      \
00078       }                                                               \
00079     break;                                                     \
00080   }                                                            \
00081 } while (0)
00082 
00083 /* Finish unpacking an fp value in semi-raw mode: the mantissa is
00084    shifted by _FP_WORKBITS but the implicit MSB is not inserted and
00085    other classification is not done.  */
00086 #define _FP_UNPACK_SEMIRAW(fs, wc, X)     _FP_FRAC_SLL_##wc(X, _FP_WORKBITS)
00087 
00088 /* A semi-raw value has overflowed to infinity.  Adjust the mantissa
00089    and exponent appropriately.  */
00090 #define _FP_OVERFLOW_SEMIRAW(fs, wc, X)                 \
00091 do {                                             \
00092   if (FP_ROUNDMODE == FP_RND_NEAREST                    \
00093       || (FP_ROUNDMODE == FP_RND_PINF && !X##_s) \
00094       || (FP_ROUNDMODE == FP_RND_MINF && X##_s)) \
00095     {                                            \
00096       X##_e = _FP_EXPMAX_##fs;                          \
00097       _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc);          \
00098     }                                            \
00099   else                                           \
00100     {                                            \
00101       X##_e = _FP_EXPMAX_##fs - 1;               \
00102       _FP_FRAC_SET_##wc(X, _FP_MAXFRAC_##wc);           \
00103     }                                            \
00104     FP_SET_EXCEPTION(FP_EX_INEXACT);                    \
00105     FP_SET_EXCEPTION(FP_EX_OVERFLOW);                   \
00106 } while (0)
00107 
00108 /* Check for a semi-raw value being a signaling NaN and raise the
00109    invalid exception if so.  */
00110 #define _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X)                    \
00111 do {                                                    \
00112   if (X##_e == _FP_EXPMAX_##fs                                 \
00113       && !_FP_FRAC_ZEROP_##wc(X)                        \
00114       && !(_FP_FRAC_HIGH_##fs(X) & _FP_QNANBIT_SH_##fs))       \
00115     FP_SET_EXCEPTION(FP_EX_INVALID);                           \
00116 } while (0)
00117 
00118 /* Choose a NaN result from an operation on two semi-raw NaN
00119    values.  */
00120 #define _FP_CHOOSENAN_SEMIRAW(fs, wc, R, X, Y, OP)                    \
00121 do {                                                           \
00122   /* _FP_CHOOSENAN expects raw values, so shift as required.  */      \
00123   _FP_FRAC_SRL_##wc(X, _FP_WORKBITS);                                 \
00124   _FP_FRAC_SRL_##wc(Y, _FP_WORKBITS);                                 \
00125   _FP_CHOOSENAN(fs, wc, R, X, Y, OP);                                 \
00126   _FP_FRAC_SLL_##wc(R, _FP_WORKBITS);                                 \
00127 } while (0)
00128 
00129 /* Test whether a biased exponent is normal (not zero or maximum).  */
00130 #define _FP_EXP_NORMAL(fs, wc, X)  (((X##_e + 1) & _FP_EXPMAX_##fs) > 1)
00131 
00132 /* Prepare to pack an fp value in semi-raw mode: the mantissa is
00133    rounded and shifted right, with the rounding possibly increasing
00134    the exponent (including changing a finite value to infinity).  */
00135 #define _FP_PACK_SEMIRAW(fs, wc, X)                            \
00136 do {                                                    \
00137   _FP_ROUND(wc, X);                                     \
00138   if (_FP_FRAC_HIGH_##fs(X)                             \
00139       & (_FP_OVERFLOW_##fs >> 1))                       \
00140     {                                                   \
00141       _FP_FRAC_HIGH_##fs(X) &= ~(_FP_OVERFLOW_##fs >> 1);      \
00142       X##_e++;                                                 \
00143       if (X##_e == _FP_EXPMAX_##fs)                            \
00144        _FP_OVERFLOW_SEMIRAW(fs, wc, X);                 \
00145     }                                                   \
00146   _FP_FRAC_SRL_##wc(X, _FP_WORKBITS);                          \
00147   if (!_FP_EXP_NORMAL(fs, wc, X) && !_FP_FRAC_ZEROP_##wc(X))   \
00148     {                                                   \
00149       if (X##_e == 0)                                          \
00150        FP_SET_EXCEPTION(FP_EX_UNDERFLOW);               \
00151       else                                              \
00152        {                                                \
00153          if (!_FP_KEEPNANFRACP)                         \
00154            {                                            \
00155              _FP_FRAC_SET_##wc(X, _FP_NANFRAC_##fs);           \
00156              X##_s = _FP_NANSIGN_##fs;                         \
00157            }                                            \
00158          else                                           \
00159            _FP_FRAC_HIGH_RAW_##fs(X) |= _FP_QNANBIT_##fs;      \
00160        }                                                \
00161     }                                                   \
00162 } while (0)
00163 
00164 /*
00165  * Before packing the bits back into the native fp result, take care
00166  * of such mundane things as rounding and overflow.  Also, for some
00167  * kinds of fp values, the original parts may not have been fully
00168  * extracted -- but that is ok, we can regenerate them now.
00169  */
00170 
00171 #define _FP_PACK_CANONICAL(fs, wc, X)                          \
00172 do {                                                    \
00173   switch (X##_c)                                        \
00174   {                                                     \
00175   case FP_CLS_NORMAL:                                          \
00176     X##_e += _FP_EXPBIAS_##fs;                                 \
00177     if (X##_e > 0)                                      \
00178       {                                                        \
00179        _FP_ROUND(wc, X);                                \
00180        if (_FP_FRAC_OVERP_##wc(fs, X))                         \
00181          {                                              \
00182            _FP_FRAC_CLEAR_OVERP_##wc(fs, X);                   \
00183            X##_e++;                                     \
00184          }                                              \
00185        _FP_FRAC_SRL_##wc(X, _FP_WORKBITS);                     \
00186        if (X##_e >= _FP_EXPMAX_##fs)                           \
00187          {                                              \
00188            /* overflow */                               \
00189            switch (FP_ROUNDMODE)                        \
00190              {                                                 \
00191              case FP_RND_NEAREST:                       \
00192               X##_c = FP_CLS_INF;                       \
00193               break;                                    \
00194              case FP_RND_PINF:                                 \
00195               if (!X##_s) X##_c = FP_CLS_INF;                  \
00196               break;                                    \
00197              case FP_RND_MINF:                                 \
00198               if (X##_s) X##_c = FP_CLS_INF;                   \
00199               break;                                    \
00200              }                                                 \
00201            if (X##_c == FP_CLS_INF)                            \
00202              {                                                 \
00203               /* Overflow to infinity */                \
00204               X##_e = _FP_EXPMAX_##fs;                  \
00205               _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc);  \
00206              }                                                 \
00207            else                                         \
00208              {                                                 \
00209               /* Overflow to maximum normal */          \
00210               X##_e = _FP_EXPMAX_##fs - 1;                     \
00211               _FP_FRAC_SET_##wc(X, _FP_MAXFRAC_##wc);          \
00212              }                                                 \
00213            FP_SET_EXCEPTION(FP_EX_OVERFLOW);                   \
00214             FP_SET_EXCEPTION(FP_EX_INEXACT);                   \
00215          }                                              \
00216       }                                                        \
00217     else                                                \
00218       {                                                        \
00219        /* we've got a denormalized number */                   \
00220        X##_e = -X##_e + 1;                              \
00221        if (X##_e <= _FP_WFRACBITS_##fs)                 \
00222          {                                              \
00223            _FP_FRAC_SRS_##wc(X, X##_e, _FP_WFRACBITS_##fs);    \
00224            _FP_ROUND(wc, X);                                   \
00225            if (_FP_FRAC_HIGH_##fs(X)                           \
00226               & (_FP_OVERFLOW_##fs >> 1))               \
00227              {                                                 \
00228                X##_e = 1;                               \
00229                _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \
00230              }                                                 \
00231            else                                         \
00232              {                                                 \
00233               X##_e = 0;                                \
00234               _FP_FRAC_SRL_##wc(X, _FP_WORKBITS);              \
00235               FP_SET_EXCEPTION(FP_EX_UNDERFLOW);        \
00236              }                                                 \
00237          }                                              \
00238        else                                             \
00239          {                                              \
00240            /* underflow to zero */                      \
00241            X##_e = 0;                                          \
00242            if (!_FP_FRAC_ZEROP_##wc(X))                 \
00243              {                                                 \
00244                _FP_FRAC_SET_##wc(X, _FP_MINFRAC_##wc);         \
00245                _FP_ROUND(wc, X);                        \
00246                _FP_FRAC_LOW_##wc(X) >>= (_FP_WORKBITS); \
00247              }                                                 \
00248            FP_SET_EXCEPTION(FP_EX_UNDERFLOW);                  \
00249          }                                              \
00250       }                                                        \
00251     break;                                              \
00252                                                         \
00253   case FP_CLS_ZERO:                                     \
00254     X##_e = 0;                                                 \
00255     _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc);                   \
00256     break;                                              \
00257                                                         \
00258   case FP_CLS_INF:                                      \
00259     X##_e = _FP_EXPMAX_##fs;                                   \
00260     _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc);                   \
00261     break;                                              \
00262                                                         \
00263   case FP_CLS_NAN:                                      \
00264     X##_e = _FP_EXPMAX_##fs;                                   \
00265     if (!_FP_KEEPNANFRACP)                              \
00266       {                                                        \
00267        _FP_FRAC_SET_##wc(X, _FP_NANFRAC_##fs);                 \
00268        X##_s = _FP_NANSIGN_##fs;                        \
00269       }                                                        \
00270     else                                                \
00271       _FP_FRAC_HIGH_RAW_##fs(X) |= _FP_QNANBIT_##fs;           \
00272     break;                                              \
00273   }                                                     \
00274 } while (0)
00275 
00276 /* This one accepts raw argument and not cooked,  returns
00277  * 1 if X is a signaling NaN.
00278  */
00279 #define _FP_ISSIGNAN(fs, wc, X)                                \
00280 ({                                                      \
00281   int __ret = 0;                                        \
00282   if (X##_e == _FP_EXPMAX_##fs)                                \
00283     {                                                   \
00284       if (!_FP_FRAC_ZEROP_##wc(X)                       \
00285          && !(_FP_FRAC_HIGH_RAW_##fs(X) & _FP_QNANBIT_##fs))   \
00286        __ret = 1;                                       \
00287     }                                                   \
00288   __ret;                                                \
00289 })
00290 
00291 
00292 
00293 
00294 
00295 /* Addition on semi-raw values.  */
00296 #define _FP_ADD_INTERNAL(fs, wc, R, X, Y, OP)                          \
00297 do {                                                            \
00298   if (X##_s == Y##_s)                                                  \
00299     {                                                           \
00300       /* Addition.  */                                                 \
00301       R##_s = X##_s;                                            \
00302       int ediff = X##_e - Y##_e;                                \
00303       if (ediff > 0)                                            \
00304        {                                                        \
00305          R##_e = X##_e;                                         \
00306          if (Y##_e == 0)                                        \
00307            {                                                    \
00308              /* Y is zero or denormalized.  */                         \
00309              if (_FP_FRAC_ZEROP_##wc(Y))                        \
00310               {                                                 \
00311                 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X);                   \
00312                 _FP_FRAC_COPY_##wc(R, X);                       \
00313                 goto add_done;                                  \
00314               }                                                 \
00315              else                                               \
00316               {                                                 \
00317                 FP_SET_EXCEPTION(FP_EX_DENORM);                 \
00318                 ediff--;                                        \
00319                 if (ediff == 0)                                 \
00320                   {                                             \
00321                     _FP_FRAC_ADD_##wc(R, X, Y);                 \
00322                     goto add3;                                  \
00323                   }                                             \
00324                 if (X##_e == _FP_EXPMAX_##fs)                          \
00325                   {                                             \
00326                     _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X);               \
00327                     _FP_FRAC_COPY_##wc(R, X);                          \
00328                     goto add_done;                              \
00329                   }                                             \
00330                 goto add1;                                      \
00331               }                                                 \
00332            }                                                    \
00333          else if (X##_e == _FP_EXPMAX_##fs)                            \
00334            {                                                    \
00335              /* X is NaN or Inf, Y is normal.  */                      \
00336              _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X);                      \
00337              _FP_FRAC_COPY_##wc(R, X);                                 \
00338              goto add_done;                                     \
00339            }                                                    \
00340                                                                 \
00341          /* Insert implicit MSB of Y.  */                       \
00342          _FP_FRAC_HIGH_##fs(Y) |= _FP_IMPLBIT_SH_##fs;                 \
00343                                                                 \
00344        add1:                                                    \
00345          /* Shift the mantissa of Y to the right EDIFF steps;          \
00346             remember to account later for the implicit MSB of X.  */   \
00347          if (ediff <= _FP_WFRACBITS_##fs)                       \
00348            _FP_FRAC_SRS_##wc(Y, ediff, _FP_WFRACBITS_##fs);            \
00349          else if (!_FP_FRAC_ZEROP_##wc(Y))                             \
00350            _FP_FRAC_SET_##wc(Y, _FP_MINFRAC_##wc);                     \
00351          _FP_FRAC_ADD_##wc(R, X, Y);                                   \
00352        }                                                        \
00353       else if (ediff < 0)                                       \
00354        {                                                        \
00355          ediff = -ediff;                                        \
00356          R##_e = Y##_e;                                         \
00357          if (X##_e == 0)                                        \
00358            {                                                    \
00359              /* X is zero or denormalized.  */                         \
00360              if (_FP_FRAC_ZEROP_##wc(X))                        \
00361               {                                                 \
00362                 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y);                   \
00363                 _FP_FRAC_COPY_##wc(R, Y);                       \
00364                 goto add_done;                                  \
00365               }                                                 \
00366              else                                               \
00367               {                                                 \
00368                 FP_SET_EXCEPTION(FP_EX_DENORM);                 \
00369                 ediff--;                                        \
00370                 if (ediff == 0)                                 \
00371                   {                                             \
00372                     _FP_FRAC_ADD_##wc(R, Y, X);                 \
00373                     goto add3;                                  \
00374                   }                                             \
00375                 if (Y##_e == _FP_EXPMAX_##fs)                          \
00376                   {                                             \
00377                     _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y);               \
00378                     _FP_FRAC_COPY_##wc(R, Y);                          \
00379                     goto add_done;                              \
00380                   }                                             \
00381                 goto add2;                                      \
00382               }                                                 \
00383            }                                                    \
00384          else if (Y##_e == _FP_EXPMAX_##fs)                            \
00385            {                                                    \
00386              /* Y is NaN or Inf, X is normal.  */                      \
00387              _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y);                      \
00388              _FP_FRAC_COPY_##wc(R, Y);                                 \
00389              goto add_done;                                     \
00390            }                                                    \
00391                                                                 \
00392          /* Insert implicit MSB of X.  */                       \
00393          _FP_FRAC_HIGH_##fs(X) |= _FP_IMPLBIT_SH_##fs;                 \
00394                                                                 \
00395        add2:                                                    \
00396          /* Shift the mantissa of X to the right EDIFF steps;          \
00397             remember to account later for the implicit MSB of Y.  */   \
00398          if (ediff <= _FP_WFRACBITS_##fs)                       \
00399            _FP_FRAC_SRS_##wc(X, ediff, _FP_WFRACBITS_##fs);            \
00400          else if (!_FP_FRAC_ZEROP_##wc(X))                             \
00401            _FP_FRAC_SET_##wc(X, _FP_MINFRAC_##wc);                     \
00402          _FP_FRAC_ADD_##wc(R, Y, X);                                   \
00403        }                                                        \
00404       else                                                      \
00405        {                                                        \
00406          /* ediff == 0.  */                                     \
00407          if (!_FP_EXP_NORMAL(fs, wc, X))                        \
00408            {                                                    \
00409              if (X##_e == 0)                                           \
00410               {                                                 \
00411                 /* X and Y are zero or denormalized.  */               \
00412                 R##_e = 0;                                      \
00413                 if (_FP_FRAC_ZEROP_##wc(X))                            \
00414                   {                                             \
00415                     if (!_FP_FRAC_ZEROP_##wc(Y))                \
00416                      FP_SET_EXCEPTION(FP_EX_DENORM);                   \
00417                     _FP_FRAC_COPY_##wc(R, Y);                          \
00418                     goto add_done;                              \
00419                   }                                             \
00420                 else if (_FP_FRAC_ZEROP_##wc(Y))                \
00421                   {                                             \
00422                     FP_SET_EXCEPTION(FP_EX_DENORM);                    \
00423                     _FP_FRAC_COPY_##wc(R, X);                          \
00424                     goto add_done;                              \
00425                   }                                             \
00426                 else                                            \
00427                   {                                             \
00428                     FP_SET_EXCEPTION(FP_EX_DENORM);                    \
00429                     _FP_FRAC_ADD_##wc(R, X, Y);                 \
00430                     if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs)   \
00431                      {                                          \
00432                        /* Normalized result.  */                \
00433                        _FP_FRAC_HIGH_##fs(R)                           \
00434                          &= ~(_FP_W_TYPE)_FP_IMPLBIT_SH_##fs;   \
00435                        R##_e = 1;                               \
00436                      }                                          \
00437                     goto add_done;                              \
00438                   }                                             \
00439               }                                                 \
00440              else                                               \
00441               {                                                 \
00442                 /* X and Y are NaN or Inf.  */                  \
00443                 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X);                   \
00444                 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y);                   \
00445                 R##_e = _FP_EXPMAX_##fs;                        \
00446                 if (_FP_FRAC_ZEROP_##wc(X))                            \
00447                   _FP_FRAC_COPY_##wc(R, Y);                            \
00448                 else if (_FP_FRAC_ZEROP_##wc(Y))                \
00449                   _FP_FRAC_COPY_##wc(R, X);                            \
00450                 else                                            \
00451                   _FP_CHOOSENAN_SEMIRAW(fs, wc, R, X, Y, OP);          \
00452                 goto add_done;                                  \
00453               }                                                 \
00454            }                                                    \
00455          /* The exponents of X and Y, both normal, are equal.  The     \
00456             implicit MSBs will always add to increase the              \
00457             exponent.  */                                       \
00458          _FP_FRAC_ADD_##wc(R, X, Y);                                   \
00459          R##_e = X##_e + 1;                                     \
00460          _FP_FRAC_SRS_##wc(R, 1, _FP_WFRACBITS_##fs);                  \
00461          if (R##_e == _FP_EXPMAX_##fs)                                 \
00462            /* Overflow to infinity (depending on rounding mode).  */   \
00463            _FP_OVERFLOW_SEMIRAW(fs, wc, R);                            \
00464          goto add_done;                                         \
00465        }                                                        \
00466     add3:                                                       \
00467       if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs)                 \
00468        {                                                        \
00469          /* Overflow.  */                                       \
00470          _FP_FRAC_HIGH_##fs(R) &= ~(_FP_W_TYPE)_FP_IMPLBIT_SH_##fs;    \
00471          R##_e++;                                               \
00472          _FP_FRAC_SRS_##wc(R, 1, _FP_WFRACBITS_##fs);                  \
00473          if (R##_e == _FP_EXPMAX_##fs)                                 \
00474            /* Overflow to infinity (depending on rounding mode).  */   \
00475            _FP_OVERFLOW_SEMIRAW(fs, wc, R);                            \
00476        }                                                        \
00477     add_done: ;                                                        \
00478     }                                                           \
00479   else                                                          \
00480     {                                                           \
00481       /* Subtraction.  */                                       \
00482       int ediff = X##_e - Y##_e;                                \
00483       if (ediff > 0)                                            \
00484        {                                                        \
00485          R##_e = X##_e;                                         \
00486          R##_s = X##_s;                                         \
00487          if (Y##_e == 0)                                        \
00488            {                                                    \
00489              /* Y is zero or denormalized.  */                         \
00490              if (_FP_FRAC_ZEROP_##wc(Y))                        \
00491               {                                                 \
00492                 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X);                   \
00493                 _FP_FRAC_COPY_##wc(R, X);                       \
00494                 goto sub_done;                                  \
00495               }                                                 \
00496              else                                               \
00497               {                                                 \
00498                 FP_SET_EXCEPTION(FP_EX_DENORM);                 \
00499                 ediff--;                                        \
00500                 if (ediff == 0)                                 \
00501                   {                                             \
00502                     _FP_FRAC_SUB_##wc(R, X, Y);                 \
00503                     goto sub3;                                  \
00504                   }                                             \
00505                 if (X##_e == _FP_EXPMAX_##fs)                          \
00506                   {                                             \
00507                     _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X);               \
00508                     _FP_FRAC_COPY_##wc(R, X);                          \
00509                     goto sub_done;                              \
00510                   }                                             \
00511                 goto sub1;                                      \
00512               }                                                 \
00513            }                                                    \
00514          else if (X##_e == _FP_EXPMAX_##fs)                            \
00515            {                                                    \
00516              /* X is NaN or Inf, Y is normal.  */                      \
00517              _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X);                      \
00518              _FP_FRAC_COPY_##wc(R, X);                                 \
00519              goto sub_done;                                     \
00520            }                                                    \
00521                                                                 \
00522          /* Insert implicit MSB of Y.  */                       \
00523          _FP_FRAC_HIGH_##fs(Y) |= _FP_IMPLBIT_SH_##fs;                 \
00524                                                                 \
00525        sub1:                                                    \
00526          /* Shift the mantissa of Y to the right EDIFF steps;          \
00527             remember to account later for the implicit MSB of X.  */   \
00528          if (ediff <= _FP_WFRACBITS_##fs)                       \
00529            _FP_FRAC_SRS_##wc(Y, ediff, _FP_WFRACBITS_##fs);            \
00530          else if (!_FP_FRAC_ZEROP_##wc(Y))                             \
00531            _FP_FRAC_SET_##wc(Y, _FP_MINFRAC_##wc);                     \
00532          _FP_FRAC_SUB_##wc(R, X, Y);                                   \
00533        }                                                        \
00534       else if (ediff < 0)                                       \
00535        {                                                        \
00536          ediff = -ediff;                                        \
00537          R##_e = Y##_e;                                         \
00538          R##_s = Y##_s;                                         \
00539          if (X##_e == 0)                                        \
00540            {                                                    \
00541              /* X is zero or denormalized.  */                         \
00542              if (_FP_FRAC_ZEROP_##wc(X))                        \
00543               {                                                 \
00544                 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y);                   \
00545                 _FP_FRAC_COPY_##wc(R, Y);                       \
00546                 goto sub_done;                                  \
00547               }                                                 \
00548              else                                               \
00549               {                                                 \
00550                 FP_SET_EXCEPTION(FP_EX_DENORM);                 \
00551                 ediff--;                                        \
00552                 if (ediff == 0)                                 \
00553                   {                                             \
00554                     _FP_FRAC_SUB_##wc(R, Y, X);                 \
00555                     goto sub3;                                  \
00556                   }                                             \
00557                 if (Y##_e == _FP_EXPMAX_##fs)                          \
00558                   {                                             \
00559                     _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y);               \
00560                     _FP_FRAC_COPY_##wc(R, Y);                          \
00561                     goto sub_done;                              \
00562                   }                                             \
00563                 goto sub2;                                      \
00564               }                                                 \
00565            }                                                    \
00566          else if (Y##_e == _FP_EXPMAX_##fs)                            \
00567            {                                                    \
00568              /* Y is NaN or Inf, X is normal.  */                      \
00569              _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y);                      \
00570              _FP_FRAC_COPY_##wc(R, Y);                                 \
00571              goto sub_done;                                     \
00572            }                                                    \
00573                                                                 \
00574          /* Insert implicit MSB of X.  */                       \
00575          _FP_FRAC_HIGH_##fs(X) |= _FP_IMPLBIT_SH_##fs;                 \
00576                                                                 \
00577        sub2:                                                    \
00578          /* Shift the mantissa of X to the right EDIFF steps;          \
00579             remember to account later for the implicit MSB of Y.  */   \
00580          if (ediff <= _FP_WFRACBITS_##fs)                       \
00581            _FP_FRAC_SRS_##wc(X, ediff, _FP_WFRACBITS_##fs);            \
00582          else if (!_FP_FRAC_ZEROP_##wc(X))                             \
00583            _FP_FRAC_SET_##wc(X, _FP_MINFRAC_##wc);                     \
00584          _FP_FRAC_SUB_##wc(R, Y, X);                                   \
00585        }                                                        \
00586       else                                                      \
00587        {                                                        \
00588          /* ediff == 0.  */                                     \
00589          if (!_FP_EXP_NORMAL(fs, wc, X))                        \
00590            {                                                    \
00591              if (X##_e == 0)                                           \
00592               {                                                 \
00593                 /* X and Y are zero or denormalized.  */               \
00594                 R##_e = 0;                                      \
00595                 if (_FP_FRAC_ZEROP_##wc(X))                            \
00596                   {                                             \
00597                     _FP_FRAC_COPY_##wc(R, Y);                          \
00598                     if (_FP_FRAC_ZEROP_##wc(Y))                 \
00599                      R##_s = (FP_ROUNDMODE == FP_RND_MINF);            \
00600                     else                                        \
00601                      {                                          \
00602                        FP_SET_EXCEPTION(FP_EX_DENORM);          \
00603                        R##_s = Y##_s;                           \
00604                      }                                          \
00605                     goto sub_done;                              \
00606                   }                                             \
00607                 else if (_FP_FRAC_ZEROP_##wc(Y))                \
00608                   {                                             \
00609                     FP_SET_EXCEPTION(FP_EX_DENORM);                    \
00610                     _FP_FRAC_COPY_##wc(R, X);                          \
00611                     R##_s = X##_s;                              \
00612                     goto sub_done;                              \
00613                   }                                             \
00614                 else                                            \
00615                   {                                             \
00616                     FP_SET_EXCEPTION(FP_EX_DENORM);                    \
00617                     _FP_FRAC_SUB_##wc(R, X, Y);                 \
00618                     R##_s = X##_s;                              \
00619                     if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs)   \
00620                      {                                          \
00621                        /* |X| < |Y|, negate result.  */         \
00622                        _FP_FRAC_SUB_##wc(R, Y, X);                     \
00623                        R##_s = Y##_s;                           \
00624                      }                                          \
00625                     else if (_FP_FRAC_ZEROP_##wc(R))                   \
00626                      R##_s = (FP_ROUNDMODE == FP_RND_MINF);            \
00627                     goto sub_done;                              \
00628                   }                                             \
00629               }                                                 \
00630              else                                               \
00631               {                                                 \
00632                 /* X and Y are NaN or Inf, of opposite signs.  */      \
00633                 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X);                   \
00634                 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y);                   \
00635                 R##_e = _FP_EXPMAX_##fs;                        \
00636                 if (_FP_FRAC_ZEROP_##wc(X))                            \
00637                   {                                             \
00638                     if (_FP_FRAC_ZEROP_##wc(Y))                 \
00639                      {                                          \
00640                        /* Inf - Inf.  */                        \
00641                        R##_s = _FP_NANSIGN_##fs;                \
00642                        _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs);  \
00643                        _FP_FRAC_SLL_##wc(R, _FP_WORKBITS);             \
00644                        FP_SET_EXCEPTION(FP_EX_INVALID);         \
00645                      }                                          \
00646                     else                                        \
00647                      {                                          \
00648                        /* Inf - NaN.  */                        \
00649                        R##_s = Y##_s;                           \
00650                        _FP_FRAC_COPY_##wc(R, Y);                \
00651                      }                                          \
00652                   }                                             \
00653                 else                                            \
00654                   {                                             \
00655                     if (_FP_FRAC_ZEROP_##wc(Y))                 \
00656                      {                                          \
00657                        /* NaN - Inf.  */                        \
00658                        R##_s = X##_s;                           \
00659                        _FP_FRAC_COPY_##wc(R, X);                \
00660                      }                                          \
00661                     else                                        \
00662                      {                                          \
00663                        /* NaN - NaN.  */                        \
00664                        _FP_CHOOSENAN_SEMIRAW(fs, wc, R, X, Y, OP);     \
00665                      }                                          \
00666                   }                                             \
00667                 goto sub_done;                                  \
00668               }                                                 \
00669            }                                                    \
00670          /* The exponents of X and Y, both normal, are equal.  The     \
00671             implicit MSBs cancel.  */                                  \
00672          R##_e = X##_e;                                         \
00673          _FP_FRAC_SUB_##wc(R, X, Y);                                   \
00674          R##_s = X##_s;                                         \
00675          if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs)              \
00676            {                                                    \
00677              /* |X| < |Y|, negate result.  */                          \
00678              _FP_FRAC_SUB_##wc(R, Y, X);                        \
00679              R##_s = Y##_s;                                     \
00680            }                                                    \
00681          else if (_FP_FRAC_ZEROP_##wc(R))                       \
00682            {                                                    \
00683              R##_e = 0;                                         \
00684              R##_s = (FP_ROUNDMODE == FP_RND_MINF);                    \
00685              goto sub_done;                                     \
00686            }                                                    \
00687          goto norm;                                             \
00688        }                                                        \
00689     sub3:                                                       \
00690       if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs)                 \
00691        {                                                        \
00692          int diff;                                              \
00693          /* Carry into most significant bit of larger one of X and Y,  \
00694             canceling it; renormalize.  */                             \
00695          _FP_FRAC_HIGH_##fs(R) &= _FP_IMPLBIT_SH_##fs - 1;             \
00696        norm:                                                    \
00697          _FP_FRAC_CLZ_##wc(diff, R);                                   \
00698          diff -= _FP_WFRACXBITS_##fs;                                  \
00699          _FP_FRAC_SLL_##wc(R, diff);                                   \
00700          if (R##_e <= diff)                                     \
00701            {                                                    \
00702              /* R is denormalized.  */                                 \
00703              diff = diff - R##_e + 1;                                  \
00704              _FP_FRAC_SRS_##wc(R, diff, _FP_WFRACBITS_##fs);           \
00705              R##_e = 0;                                         \
00706            }                                                    \
00707          else                                                   \
00708            {                                                    \
00709              R##_e -= diff;                                     \
00710              _FP_FRAC_HIGH_##fs(R) &= ~(_FP_W_TYPE)_FP_IMPLBIT_SH_##fs; \
00711            }                                                    \
00712        }                                                        \
00713     sub_done: ;                                                        \
00714     }                                                           \
00715 } while (0)
00716 
00717 #define _FP_ADD(fs, wc, R, X, Y) _FP_ADD_INTERNAL(fs, wc, R, X, Y, '+')
00718 #define _FP_SUB(fs, wc, R, X, Y)                                   \
00719   do {                                                             \
00720     if (!(Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y))) Y##_s ^= 1; \
00721     _FP_ADD_INTERNAL(fs, wc, R, X, Y, '-');                               \
00722   } while (0)
00723 
00724 
00725 /*
00726  * Main negation routine.  FIXME -- when we care about setting exception
00727  * bits reliably, this will not do.  We should examine all of the fp classes.
00728  */
00729 
00730 #define _FP_NEG(fs, wc, R, X)             \
00731   do {                             \
00732     _FP_FRAC_COPY_##wc(R, X);             \
00733     R##_c = X##_c;                 \
00734     R##_e = X##_e;                 \
00735     R##_s = 1 ^ X##_s;                    \
00736   } while (0)
00737 
00738 
00739 /*
00740  * Main multiplication routine.  The input values should be cooked.
00741  */
00742 
00743 #define _FP_MUL(fs, wc, R, X, Y)                 \
00744 do {                                             \
00745   R##_s = X##_s ^ Y##_s;                         \
00746   switch (_FP_CLS_COMBINE(X##_c, Y##_c))         \
00747   {                                              \
00748   case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NORMAL):    \
00749     R##_c = FP_CLS_NORMAL;                       \
00750     R##_e = X##_e + Y##_e + 1;                          \
00751                                                  \
00752     _FP_MUL_MEAT_##fs(R,X,Y);                           \
00753                                                  \
00754     if (_FP_FRAC_OVERP_##wc(fs, R))                     \
00755       _FP_FRAC_SRS_##wc(R, 1, _FP_WFRACBITS_##fs);      \
00756     else                                         \
00757       R##_e--;                                          \
00758     break;                                       \
00759                                                  \
00760   case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NAN):          \
00761     _FP_CHOOSENAN(fs, wc, R, X, Y, '*');         \
00762     break;                                       \
00763                                                  \
00764   case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NORMAL):       \
00765   case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_INF):          \
00766   case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_ZERO):         \
00767     R##_s = X##_s;                               \
00768                                                  \
00769   case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_INF):          \
00770   case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NORMAL):       \
00771   case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NORMAL):      \
00772   case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_ZERO): \
00773     _FP_FRAC_COPY_##wc(R, X);                           \
00774     R##_c = X##_c;                               \
00775     break;                                       \
00776                                                  \
00777   case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NAN):       \
00778   case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NAN):          \
00779   case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NAN):         \
00780     R##_s = Y##_s;                               \
00781                                                  \
00782   case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_INF):       \
00783   case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_ZERO):      \
00784     _FP_FRAC_COPY_##wc(R, Y);                           \
00785     R##_c = Y##_c;                               \
00786     break;                                       \
00787                                                  \
00788   case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_ZERO):         \
00789   case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_INF):         \
00790     R##_s = _FP_NANSIGN_##fs;                           \
00791     R##_c = FP_CLS_NAN;                                 \
00792     _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs);             \
00793     FP_SET_EXCEPTION(FP_EX_INVALID);                    \
00794     break;                                       \
00795                                                  \
00796   default:                                       \
00797     abort();                                     \
00798   }                                              \
00799 } while (0)
00800 
00801 
00802 /*
00803  * Main division routine.  The input values should be cooked.
00804  */
00805 
00806 #define _FP_DIV(fs, wc, R, X, Y)                 \
00807 do {                                             \
00808   R##_s = X##_s ^ Y##_s;                         \
00809   switch (_FP_CLS_COMBINE(X##_c, Y##_c))         \
00810   {                                              \
00811   case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NORMAL):    \
00812     R##_c = FP_CLS_NORMAL;                       \
00813     R##_e = X##_e - Y##_e;                       \
00814                                                  \
00815     _FP_DIV_MEAT_##fs(R,X,Y);                           \
00816     break;                                       \
00817                                                  \
00818   case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NAN):          \
00819     _FP_CHOOSENAN(fs, wc, R, X, Y, '/');         \
00820     break;                                       \
00821                                                  \
00822   case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NORMAL):       \
00823   case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_INF):          \
00824   case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_ZERO):         \
00825     R##_s = X##_s;                               \
00826     _FP_FRAC_COPY_##wc(R, X);                           \
00827     R##_c = X##_c;                               \
00828     break;                                       \
00829                                                  \
00830   case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NAN):       \
00831   case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NAN):          \
00832   case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NAN):         \
00833     R##_s = Y##_s;                               \
00834     _FP_FRAC_COPY_##wc(R, Y);                           \
00835     R##_c = Y##_c;                               \
00836     break;                                       \
00837                                                  \
00838   case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_INF):       \
00839   case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_INF):         \
00840   case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NORMAL):      \
00841     R##_c = FP_CLS_ZERO;                         \
00842     break;                                       \
00843                                                  \
00844   case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_ZERO):      \
00845     FP_SET_EXCEPTION(FP_EX_DIVZERO);                    \
00846   case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_ZERO):         \
00847   case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NORMAL):       \
00848     R##_c = FP_CLS_INF;                                 \
00849     break;                                       \
00850                                                  \
00851   case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_INF):          \
00852   case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_ZERO): \
00853     R##_s = _FP_NANSIGN_##fs;                           \
00854     R##_c = FP_CLS_NAN;                                 \
00855     _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs);             \
00856     FP_SET_EXCEPTION(FP_EX_INVALID);                    \
00857     break;                                       \
00858                                                  \
00859   default:                                       \
00860     abort();                                     \
00861   }                                              \
00862 } while (0)
00863 
00864 
00865 /*
00866  * Main differential comparison routine.  The inputs should be raw not
00867  * cooked.  The return is -1,0,1 for normal values, 2 otherwise.
00868  */
00869 
00870 #define _FP_CMP(fs, wc, ret, X, Y, un)                                \
00871   do {                                                         \
00872     /* NANs are unordered */                                          \
00873     if ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(X))         \
00874        || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y)))      \
00875       {                                                               \
00876        ret = un;                                               \
00877       }                                                               \
00878     else                                                       \
00879       {                                                               \
00880        int __is_zero_x;                                        \
00881        int __is_zero_y;                                        \
00882                                                                \
00883        __is_zero_x = (!X##_e && _FP_FRAC_ZEROP_##wc(X)) ? 1 : 0;      \
00884        __is_zero_y = (!Y##_e && _FP_FRAC_ZEROP_##wc(Y)) ? 1 : 0;      \
00885                                                                \
00886        if (__is_zero_x && __is_zero_y)                                \
00887               ret = 0;                                         \
00888        else if (__is_zero_x)                                          \
00889               ret = Y##_s ? 1 : -1;                                   \
00890        else if (__is_zero_y)                                          \
00891               ret = X##_s ? -1 : 1;                                   \
00892        else if (X##_s != Y##_s)                                \
00893          ret = X##_s ? -1 : 1;                                        \
00894        else if (X##_e > Y##_e)                                        \
00895          ret = X##_s ? -1 : 1;                                        \
00896        else if (X##_e < Y##_e)                                        \
00897          ret = X##_s ? 1 : -1;                                        \
00898        else if (_FP_FRAC_GT_##wc(X, Y))                        \
00899          ret = X##_s ? -1 : 1;                                        \
00900        else if (_FP_FRAC_GT_##wc(Y, X))                        \
00901          ret = X##_s ? 1 : -1;                                        \
00902        else                                                    \
00903          ret = 0;                                              \
00904       }                                                               \
00905   } while (0)
00906 
00907 
00908 /* Simplification for strict equality.  */
00909 
00910 #define _FP_CMP_EQ(fs, wc, ret, X, Y)                                     \
00911   do {                                                             \
00912     /* NANs are unordered */                                              \
00913     if ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(X))             \
00914        || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y)))          \
00915       {                                                                   \
00916        ret = 1;                                                    \
00917       }                                                                   \
00918     else                                                           \
00919       {                                                                   \
00920        ret = !(X##_e == Y##_e                                             \
00921               && _FP_FRAC_EQ_##wc(X, Y)                            \
00922               && (X##_s == Y##_s || (!X##_e && _FP_FRAC_ZEROP_##wc(X)))); \
00923       }                                                                   \
00924   } while (0)
00925 
00926 /* Version to test unordered.  */
00927 
00928 #define _FP_CMP_UNORD(fs, wc, ret, X, Y)                       \
00929   do {                                                         \
00930     ret = ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(X))      \
00931           || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y)));  \
00932   } while (0)
00933 
00934 /*
00935  * Main square root routine.  The input value should be cooked.
00936  */
00937 
00938 #define _FP_SQRT(fs, wc, R, X)                                        \
00939 do {                                                           \
00940     _FP_FRAC_DECL_##wc(T); _FP_FRAC_DECL_##wc(S);                     \
00941     _FP_W_TYPE q;                                              \
00942     switch (X##_c)                                             \
00943     {                                                          \
00944     case FP_CLS_NAN:                                           \
00945        _FP_FRAC_COPY_##wc(R, X);                               \
00946        R##_s = X##_s;                                                 \
00947        R##_c = FP_CLS_NAN;                                     \
00948        break;                                                  \
00949     case FP_CLS_INF:                                           \
00950        if (X##_s)                                              \
00951          {                                                     \
00952            R##_s = _FP_NANSIGN_##fs;                                  \
00953            R##_c = FP_CLS_NAN; /* NAN */                       \
00954            _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs);                    \
00955            FP_SET_EXCEPTION(FP_EX_INVALID);                           \
00956          }                                                     \
00957        else                                                    \
00958          {                                                     \
00959            R##_s = 0;                                                 \
00960            R##_c = FP_CLS_INF; /* sqrt(+inf) = +inf */                \
00961          }                                                     \
00962        break;                                                  \
00963     case FP_CLS_ZERO:                                                 \
00964        R##_s = X##_s;                                                 \
00965        R##_c = FP_CLS_ZERO; /* sqrt(+-0) = +-0 */                     \
00966        break;                                                  \
00967     case FP_CLS_NORMAL:                                               \
00968        R##_s = 0;                                              \
00969         if (X##_s)                                             \
00970           {                                                    \
00971            R##_c = FP_CLS_NAN; /* sNAN */                      \
00972            R##_s = _FP_NANSIGN_##fs;                                  \
00973            _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs);                    \
00974            FP_SET_EXCEPTION(FP_EX_INVALID);                           \
00975            break;                                              \
00976           }                                                    \
00977        R##_c = FP_CLS_NORMAL;                                         \
00978         if (X##_e & 1)                                                \
00979           _FP_FRAC_SLL_##wc(X, 1);                             \
00980         R##_e = X##_e >> 1;                                    \
00981         _FP_FRAC_SET_##wc(S, _FP_ZEROFRAC_##wc);               \
00982         _FP_FRAC_SET_##wc(R, _FP_ZEROFRAC_##wc);               \
00983         q = _FP_OVERFLOW_##fs >> 1;                                   \
00984         _FP_SQRT_MEAT_##wc(R, S, T, X, q);                            \
00985     }                                                          \
00986   } while (0)
00987 
00988 /*
00989  * Convert from FP to integer.  Input is raw.
00990  */
00991 
00992 /* RSIGNED can have following values:
00993  * 0:  the number is required to be 0..(2^rsize)-1, if not, NV is set plus
00994  *     the result is either 0 or (2^rsize)-1 depending on the sign in such
00995  *     case.
00996  * 1:  the number is required to be -(2^(rsize-1))..(2^(rsize-1))-1, if not,
00997  *     NV is set plus the result is either -(2^(rsize-1)) or (2^(rsize-1))-1
00998  *     depending on the sign in such case.
00999  * -1: the number is required to be -(2^(rsize-1))..(2^rsize)-1, if not, NV is
01000  *     set plus the result is either -(2^(rsize-1)) or (2^(rsize-1))-1
01001  *     depending on the sign in such case.
01002  */
01003 #define _FP_TO_INT(fs, wc, r, X, rsize, rsigned)               \
01004 do {                                                           \
01005   if (X##_e < _FP_EXPBIAS_##fs)                                       \
01006     {                                                          \
01007       r = 0;                                                   \
01008       if (X##_e == 0)                                                 \
01009        {                                                       \
01010          if (!_FP_FRAC_ZEROP_##wc(X))                                 \
01011            {                                                   \
01012              FP_SET_EXCEPTION(FP_EX_INEXACT);                         \
01013              FP_SET_EXCEPTION(FP_EX_DENORM);                          \
01014            }                                                   \
01015        }                                                       \
01016       else                                                     \
01017        FP_SET_EXCEPTION(FP_EX_INEXACT);                        \
01018     }                                                          \
01019   else if (X##_e >= _FP_EXPBIAS_##fs + rsize - (rsigned > 0 || X##_s) \
01020           || (!rsigned && X##_s))                              \
01021     {                                                          \
01022       /* Overflow or converting to the most negative integer.  */     \
01023       if (rsigned)                                             \
01024        {                                                       \
01025          r = 1;                                                \
01026          r <<= rsize - 1;                                      \
01027          r -= 1 - X##_s;                                       \
01028        } else {                                                \
01029          r = 0;                                                \
01030          if (X##_s)                                            \
01031            r = ~r;                                             \
01032        }                                                       \
01033                                                                \
01034       if (rsigned && X##_s && X##_e == _FP_EXPBIAS_##fs + rsize - 1)  \
01035        {                                                       \
01036          /* Possibly converting to most negative integer; check the   \
01037             mantissa.  */                                      \
01038          int inexact = 0;                                      \
01039          (void)((_FP_FRACBITS_##fs > rsize)                           \
01040                ? ({ _FP_FRAC_SRST_##wc(X, inexact,                    \
01041                                     _FP_FRACBITS_##fs - rsize, \
01042                                     _FP_FRACBITS_##fs); 0; })  \
01043                : 0);                                           \
01044          if (!_FP_FRAC_ZEROP_##wc(X))                                 \
01045            FP_SET_EXCEPTION(FP_EX_INVALID);                           \
01046          else if (inexact)                                     \
01047            FP_SET_EXCEPTION(FP_EX_INEXACT);                           \
01048        }                                                       \
01049       else                                                     \
01050        FP_SET_EXCEPTION(FP_EX_INVALID);                        \
01051     }                                                          \
01052   else                                                         \
01053     {                                                          \
01054       _FP_FRAC_HIGH_RAW_##fs(X) |= _FP_IMPLBIT_##fs;                  \
01055       if (X##_e >= _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs - 1)          \
01056        {                                                       \
01057          _FP_FRAC_ASSEMBLE_##wc(r, X, rsize);                         \
01058          r <<= X##_e - _FP_EXPBIAS_##fs - _FP_FRACBITS_##fs + 1;      \
01059        }                                                       \
01060       else                                                     \
01061        {                                                       \
01062          int inexact;                                                 \
01063          _FP_FRAC_SRST_##wc(X, inexact,                        \
01064                          (_FP_FRACBITS_##fs + _FP_EXPBIAS_##fs - 1    \
01065                           - X##_e),                                   \
01066                          _FP_FRACBITS_##fs);                          \
01067          if (inexact)                                                 \
01068            FP_SET_EXCEPTION(FP_EX_INEXACT);                           \
01069          _FP_FRAC_ASSEMBLE_##wc(r, X, rsize);                         \
01070        }                                                       \
01071       if (rsigned && X##_s)                                    \
01072        r = -r;                                                        \
01073     }                                                          \
01074 } while (0)
01075 
01076 /* Convert integer to fp.  Output is raw.  RTYPE is unsigned even if
01077    input is signed.  */
01078 #define _FP_FROM_INT(fs, wc, X, r, rsize, rtype)                    \
01079   do {                                                              \
01080     if (r)                                                          \
01081       {                                                                    \
01082        rtype ur_;                                                   \
01083                                                                     \
01084        if ((X##_s = (r < 0)))                                              \
01085          r = -(rtype)r;                                             \
01086                                                                     \
01087        ur_ = (rtype) r;                                             \
01088        (void)((rsize <= _FP_W_TYPE_SIZE)                            \
01089               ? ({                                                  \
01090                   int lz_;                                          \
01091                   __FP_CLZ(lz_, (_FP_W_TYPE)ur_);                          \
01092                   X##_e = _FP_EXPBIAS_##fs + _FP_W_TYPE_SIZE - 1 - lz_;    \
01093                 })                                                  \
01094               : ((rsize <= 2 * _FP_W_TYPE_SIZE)                     \
01095                 ? ({                                                \
01096                      int lz_;                                              \
01097                      __FP_CLZ_2(lz_, (_FP_W_TYPE)(ur_ >> _FP_W_TYPE_SIZE), \
01098                               (_FP_W_TYPE)ur_);                     \
01099                      X##_e = (_FP_EXPBIAS_##fs + 2 * _FP_W_TYPE_SIZE - 1   \
01100                             - lz_);                                        \
01101                    })                                                      \
01102                 : (abort(), 0)));                                   \
01103                                                                     \
01104        if (rsize - 1 + _FP_EXPBIAS_##fs >= _FP_EXPMAX_##fs                 \
01105            && X##_e >= _FP_EXPMAX_##fs)                             \
01106          {                                                          \
01107            /* Exponent too big; overflow to infinity.  (May also           \
01108               happen after rounding below.)  */                     \
01109            _FP_OVERFLOW_SEMIRAW(fs, wc, X);                                \
01110            goto pack_semiraw;                                              \
01111          }                                                          \
01112                                                                     \
01113        if (rsize <= _FP_FRACBITS_##fs                                      \
01114            || X##_e < _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs)                \
01115          {                                                          \
01116            /* Exactly representable; shift left.  */                       \
01117            _FP_FRAC_DISASSEMBLE_##wc(X, ur_, rsize);                       \
01118            _FP_FRAC_SLL_##wc(X, (_FP_EXPBIAS_##fs                          \
01119                               + _FP_FRACBITS_##fs - 1 - X##_e));           \
01120          }                                                          \
01121        else                                                         \
01122          {                                                          \
01123            /* More bits in integer than in floating type; need to          \
01124               round.  */                                            \
01125            if (_FP_EXPBIAS_##fs + _FP_WFRACBITS_##fs - 1 < X##_e)          \
01126              ur_ = ((ur_ >> (X##_e - _FP_EXPBIAS_##fs                      \
01127                            - _FP_WFRACBITS_##fs + 1))               \
01128                    | ((ur_ << (rsize - (X##_e - _FP_EXPBIAS_##fs           \
01129                                      - _FP_WFRACBITS_##fs + 1)))           \
01130                      != 0));                                               \
01131            _FP_FRAC_DISASSEMBLE_##wc(X, ur_, rsize);                       \
01132            if ((_FP_EXPBIAS_##fs + _FP_WFRACBITS_##fs - 1 - X##_e) > 0)     \
01133              _FP_FRAC_SLL_##wc(X, (_FP_EXPBIAS_##fs                        \
01134                                 + _FP_WFRACBITS_##fs - 1 - X##_e));        \
01135            _FP_FRAC_HIGH_##fs(X) &= ~(_FP_W_TYPE)_FP_IMPLBIT_SH_##fs;      \
01136          pack_semiraw:                                                     \
01137            _FP_PACK_SEMIRAW(fs, wc, X);                             \
01138          }                                                          \
01139       }                                                                    \
01140     else                                                            \
01141       {                                                                    \
01142        X##_s = 0;                                                   \
01143        X##_e = 0;                                                   \
01144        _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc);                     \
01145       }                                                                    \
01146   } while (0)
01147 
01148 
01149 /* Extend from a narrower floating-point format to a wider one.  Input
01150    and output are raw.  */
01151 #define FP_EXTEND(dfs,sfs,dwc,swc,D,S)                                 \
01152 do {                                                            \
01153   if (_FP_FRACBITS_##dfs < _FP_FRACBITS_##sfs                          \
01154       || (_FP_EXPMAX_##dfs - _FP_EXPBIAS_##dfs                         \
01155          < _FP_EXPMAX_##sfs - _FP_EXPBIAS_##sfs)                \
01156       || (_FP_EXPBIAS_##dfs < _FP_EXPBIAS_##sfs + _FP_FRACBITS_##sfs - 1 \
01157          && _FP_EXPBIAS_##dfs != _FP_EXPBIAS_##sfs))                   \
01158     abort();                                                    \
01159   D##_s = S##_s;                                                \
01160   _FP_FRAC_COPY_##dwc##_##swc(D, S);                                   \
01161   if (_FP_EXP_NORMAL(sfs, swc, S))                              \
01162     {                                                           \
01163       D##_e = S##_e + _FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs;           \
01164       _FP_FRAC_SLL_##dwc(D, (_FP_FRACBITS_##dfs - _FP_FRACBITS_##sfs));       \
01165     }                                                           \
01166   else                                                          \
01167     {                                                           \
01168       if (S##_e == 0)                                                  \
01169        {                                                        \
01170          if (_FP_FRAC_ZEROP_##swc(S))                                  \
01171            D##_e = 0;                                                  \
01172          else if (_FP_EXPBIAS_##dfs                                    \
01173                  < _FP_EXPBIAS_##sfs + _FP_FRACBITS_##sfs - 1)  \
01174            {                                                    \
01175              FP_SET_EXCEPTION(FP_EX_DENORM);                           \
01176              _FP_FRAC_SLL_##dwc(D, (_FP_FRACBITS_##dfs                 \
01177                                  - _FP_FRACBITS_##sfs));               \
01178              D##_e = 0;                                         \
01179            }                                                    \
01180          else                                                   \
01181            {                                                    \
01182              int _lz;                                                  \
01183              FP_SET_EXCEPTION(FP_EX_DENORM);                           \
01184              _FP_FRAC_CLZ_##swc(_lz, S);                        \
01185              _FP_FRAC_SLL_##dwc(D,                              \
01186                              _lz + _FP_FRACBITS_##dfs           \
01187                              - _FP_FRACTBITS_##sfs);            \
01188              D##_e = (_FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs + 1        \
01189                      + _FP_FRACXBITS_##sfs - _lz);                     \
01190            }                                                    \
01191        }                                                        \
01192       else                                                      \
01193        {                                                        \
01194          D##_e = _FP_EXPMAX_##dfs;                              \
01195          if (!_FP_FRAC_ZEROP_##swc(S))                                 \
01196            {                                                    \
01197              if (!(_FP_FRAC_HIGH_RAW_##sfs(S) & _FP_QNANBIT_##sfs))    \
01198               FP_SET_EXCEPTION(FP_EX_INVALID);                  \
01199              _FP_FRAC_SLL_##dwc(D, (_FP_FRACBITS_##dfs                 \
01200                                  - _FP_FRACBITS_##sfs));               \
01201            }                                                    \
01202        }                                                        \
01203     }                                                           \
01204 } while (0)
01205 
01206 /* Truncate from a wider floating-point format to a narrower one.
01207    Input and output are semi-raw.  */
01208 #define FP_TRUNC(dfs,sfs,dwc,swc,D,S)                                      \
01209 do {                                                                \
01210   if (_FP_FRACBITS_##sfs < _FP_FRACBITS_##dfs                              \
01211       || (_FP_EXPBIAS_##sfs < _FP_EXPBIAS_##dfs + _FP_FRACBITS_##dfs - 1     \
01212          && _FP_EXPBIAS_##sfs != _FP_EXPBIAS_##dfs))                       \
01213     abort();                                                        \
01214   D##_s = S##_s;                                                    \
01215   if (_FP_EXP_NORMAL(sfs, swc, S))                                  \
01216     {                                                               \
01217       D##_e = S##_e + _FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs;               \
01218       if (D##_e >= _FP_EXPMAX_##dfs)                                       \
01219        _FP_OVERFLOW_SEMIRAW(dfs, dwc, D);                           \
01220       else                                                          \
01221        {                                                            \
01222          if (D##_e <= 0)                                            \
01223            {                                                        \
01224              if (D##_e < 1 - _FP_FRACBITS_##dfs)                    \
01225               {                                                     \
01226                 _FP_FRAC_SET_##swc(S, _FP_ZEROFRAC_##swc);                 \
01227                 _FP_FRAC_LOW_##swc(S) |= 1;                                \
01228               }                                                     \
01229              else                                                   \
01230               {                                                     \
01231                 _FP_FRAC_HIGH_##sfs(S) |= _FP_IMPLBIT_SH_##sfs;            \
01232                 _FP_FRAC_SRS_##swc(S, (_FP_WFRACBITS_##sfs                 \
01233                                     - _FP_WFRACBITS_##dfs + 1 - D##_e), \
01234                                  _FP_WFRACBITS_##sfs);              \
01235               }                                                     \
01236              D##_e = 0;                                             \
01237            }                                                        \
01238          else                                                       \
01239            _FP_FRAC_SRS_##swc(S, (_FP_WFRACBITS_##sfs                      \
01240                                - _FP_WFRACBITS_##dfs),              \
01241                             _FP_WFRACBITS_##sfs);                          \
01242          _FP_FRAC_COPY_##dwc##_##swc(D, S);                                \
01243        }                                                            \
01244     }                                                               \
01245   else                                                              \
01246     {                                                               \
01247       if (S##_e == 0)                                                      \
01248        {                                                            \
01249          D##_e = 0;                                                 \
01250          if (_FP_FRAC_ZEROP_##swc(S))                                      \
01251            _FP_FRAC_SET_##dwc(D, _FP_ZEROFRAC_##dwc);                      \
01252          else                                                       \
01253            {                                                        \
01254              FP_SET_EXCEPTION(FP_EX_DENORM);                               \
01255              if (_FP_EXPBIAS_##sfs                                  \
01256                 < _FP_EXPBIAS_##dfs + _FP_FRACBITS_##dfs - 1)              \
01257               {                                                     \
01258                 _FP_FRAC_SRS_##swc(S, (_FP_WFRACBITS_##sfs                 \
01259                                     - _FP_WFRACBITS_##dfs),         \
01260                                  _FP_WFRACBITS_##sfs);              \
01261                 _FP_FRAC_COPY_##dwc##_##swc(D, S);                         \
01262               }                                                     \
01263              else                                                   \
01264               {                                                     \
01265                 _FP_FRAC_SET_##dwc(D, _FP_ZEROFRAC_##dwc);                 \
01266                 _FP_FRAC_LOW_##dwc(D) |= 1;                                \
01267               }                                                     \
01268            }                                                        \
01269        }                                                            \
01270       else                                                          \
01271        {                                                            \
01272          D##_e = _FP_EXPMAX_##dfs;                                  \
01273          if (_FP_FRAC_ZEROP_##swc(S))                                      \
01274            _FP_FRAC_SET_##dwc(D, _FP_ZEROFRAC_##dwc);                      \
01275          else                                                       \
01276            {                                                        \
01277              _FP_CHECK_SIGNAN_SEMIRAW(sfs, swc, S);                        \
01278              _FP_FRAC_SRL_##swc(S, (_FP_WFRACBITS_##sfs             \
01279                                  - _FP_WFRACBITS_##dfs));                  \
01280              _FP_FRAC_COPY_##dwc##_##swc(D, S);                     \
01281              /* Semi-raw NaN must have all workbits cleared.  */           \
01282              _FP_FRAC_LOW_##dwc(D)                                  \
01283               &= ~(_FP_W_TYPE) ((1 << _FP_WORKBITS) - 1);                  \
01284              _FP_FRAC_HIGH_##dfs(D) |= _FP_QNANBIT_SH_##dfs;               \
01285            }                                                        \
01286        }                                                            \
01287     }                                                               \
01288 } while (0)
01289 
01290 /*
01291  * Helper primitives.
01292  */
01293 
01294 /* Count leading zeros in a word.  */
01295 
01296 #ifndef __FP_CLZ
01297 /* GCC 3.4 and later provide the builtins for us.  */
01298 #define __FP_CLZ(r, x)                                                      \
01299   do {                                                               \
01300     if (sizeof (_FP_W_TYPE) == sizeof (unsigned int))                       \
01301       r = __builtin_clz (x);                                                \
01302     else if (sizeof (_FP_W_TYPE) == sizeof (unsigned long))                 \
01303       r = __builtin_clzl (x);                                               \
01304     else if (sizeof (_FP_W_TYPE) == sizeof (unsigned long long))            \
01305       r = __builtin_clzll (x);                                              \
01306     else                                                             \
01307       abort ();                                                             \
01308   } while (0)
01309 #endif /* ndef __FP_CLZ */
01310 
01311 #define _FP_DIV_HELP_imm(q, r, n, d)             \
01312   do {                                    \
01313     q = n / d, r = n % d;                 \
01314   } while (0)
01315 
01316 
01317 /* A restoring bit-by-bit division primitive.  */
01318 
01319 #define _FP_DIV_MEAT_N_loop(fs, wc, R, X, Y)                          \
01320   do {                                                         \
01321     int count = _FP_WFRACBITS_##fs;                                   \
01322     _FP_FRAC_DECL_##wc (u);                                    \
01323     _FP_FRAC_DECL_##wc (v);                                    \
01324     _FP_FRAC_COPY_##wc (u, X);                                        \
01325     _FP_FRAC_COPY_##wc (v, Y);                                        \
01326     _FP_FRAC_SET_##wc (R, _FP_ZEROFRAC_##wc);                         \
01327     /* Normalize U and V.  */                                         \
01328     _FP_FRAC_SLL_##wc (u, _FP_WFRACXBITS_##fs);                       \
01329     _FP_FRAC_SLL_##wc (v, _FP_WFRACXBITS_##fs);                       \
01330     /* First round.  Since the operands are normalized, either the    \
01331        first or second bit will be set in the fraction.  Produce a    \
01332        normalized result by checking which and adjusting the loop     \
01333        count and exponent accordingly.  */                            \
01334     if (_FP_FRAC_GE_1 (u, v))                                         \
01335       {                                                               \
01336        _FP_FRAC_SUB_##wc (u, u, v);                                   \
01337        _FP_FRAC_LOW_##wc (R) |= 1;                             \
01338        count--;                                                \
01339       }                                                               \
01340     else                                                       \
01341       R##_e--;                                                        \
01342     /* Subsequent rounds.  */                                         \
01343     do {                                                       \
01344       int msb = (_FP_WS_TYPE) _FP_FRAC_HIGH_##wc (u) < 0;             \
01345       _FP_FRAC_SLL_##wc (u, 1);                                       \
01346       _FP_FRAC_SLL_##wc (R, 1);                                       \
01347       if (msb || _FP_FRAC_GE_1 (u, v))                                \
01348        {                                                       \
01349          _FP_FRAC_SUB_##wc (u, u, v);                                 \
01350          _FP_FRAC_LOW_##wc (R) |= 1;                                  \
01351        }                                                       \
01352     } while (--count > 0);                                     \
01353     /* If there's anything left in U, the result is inexact.  */      \
01354     _FP_FRAC_LOW_##wc (R) |= !_FP_FRAC_ZEROP_##wc (u);                \
01355   } while (0)
01356 
01357 #define _FP_DIV_MEAT_1_loop(fs, R, X, Y)  _FP_DIV_MEAT_N_loop (fs, 1, R, X, Y)
01358 #define _FP_DIV_MEAT_2_loop(fs, R, X, Y)  _FP_DIV_MEAT_N_loop (fs, 2, R, X, Y)
01359 #define _FP_DIV_MEAT_4_loop(fs, R, X, Y)  _FP_DIV_MEAT_N_loop (fs, 4, R, X, Y)