Back to index

glibc  2.9
op-2.h
Go to the documentation of this file.
00001 /* Software floating-point emulation.
00002    Basic two-word fraction declaration and manipulation.
00003    Copyright (C) 1997,1998,1999,2006,2007 Free Software Foundation, Inc.
00004    This file is part of the GNU C Library.
00005    Contributed by Richard Henderson (rth@cygnus.com),
00006                 Jakub Jelinek (jj@ultra.linux.cz),
00007                 David S. Miller (davem@redhat.com) and
00008                 Peter Maydell (pmaydell@chiark.greenend.org.uk).
00009 
00010    The GNU C Library is free software; you can redistribute it and/or
00011    modify it under the terms of the GNU Lesser General Public
00012    License as published by the Free Software Foundation; either
00013    version 2.1 of the License, or (at your option) any later version.
00014 
00015    In addition to the permissions in the GNU Lesser General Public
00016    License, the Free Software Foundation gives you unlimited
00017    permission to link the compiled version of this file into
00018    combinations with other programs, and to distribute those
00019    combinations without any restriction coming from the use of this
00020    file.  (The Lesser General Public License restrictions do apply in
00021    other respects; for example, they cover modification of the file,
00022    and distribution when not linked into a combine executable.)
00023 
00024    The GNU C Library is distributed in the hope that it will be useful,
00025    but WITHOUT ANY WARRANTY; without even the implied warranty of
00026    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00027    Lesser General Public License for more details.
00028 
00029    You should have received a copy of the GNU Lesser General Public
00030    License along with the GNU C Library; if not, write to the Free
00031    Software Foundation, 51 Franklin Street, Fifth Floor, Boston,
00032    MA 02110-1301, USA.  */
00033 
00034 #define _FP_FRAC_DECL_2(X)  _FP_W_TYPE X##_f0, X##_f1
00035 #define _FP_FRAC_COPY_2(D,S)       (D##_f0 = S##_f0, D##_f1 = S##_f1)
00036 #define _FP_FRAC_SET_2(X,I) __FP_FRAC_SET_2(X, I)
00037 #define _FP_FRAC_HIGH_2(X)  (X##_f1)
00038 #define _FP_FRAC_LOW_2(X)   (X##_f0)
00039 #define _FP_FRAC_WORD_2(X,w)       (X##_f##w)
00040 
00041 #define _FP_FRAC_SLL_2(X,N)                                        \
00042 (void)(((N) < _FP_W_TYPE_SIZE)                                            \
00043        ? ({                                                        \
00044            if (__builtin_constant_p(N) && (N) == 1)                       \
00045              {                                                            \
00046               X##_f1 = X##_f1 + X##_f1 + (((_FP_WS_TYPE)(X##_f0)) < 0);   \
00047               X##_f0 += X##_f0;                                    \
00048              }                                                            \
00049            else                                                    \
00050              {                                                            \
00051               X##_f1 = X##_f1 << (N) | X##_f0 >> (_FP_W_TYPE_SIZE - (N)); \
00052               X##_f0 <<= (N);                                             \
00053              }                                                            \
00054            0;                                                      \
00055          })                                                        \
00056        : ({                                                        \
00057            X##_f1 = X##_f0 << ((N) - _FP_W_TYPE_SIZE);                    \
00058            X##_f0 = 0;                                                    \
00059          }))
00060 
00061 
00062 #define _FP_FRAC_SRL_2(X,N)                                    \
00063 (void)(((N) < _FP_W_TYPE_SIZE)                                        \
00064        ? ({                                                    \
00065            X##_f0 = X##_f0 >> (N) | X##_f1 << (_FP_W_TYPE_SIZE - (N));       \
00066            X##_f1 >>= (N);                                     \
00067          })                                                    \
00068        : ({                                                    \
00069            X##_f0 = X##_f1 >> ((N) - _FP_W_TYPE_SIZE);                \
00070            X##_f1 = 0;                                                \
00071          }))
00072 
00073 /* Right shift with sticky-lsb.  */
00074 #define _FP_FRAC_SRST_2(X,S, N,sz)                               \
00075 (void)(((N) < _FP_W_TYPE_SIZE)                                          \
00076        ? ({                                                      \
00077            S = (__builtin_constant_p(N) && (N) == 1                     \
00078                ? X##_f0 & 1                                      \
00079                : (X##_f0 << (_FP_W_TYPE_SIZE - (N))) != 0);             \
00080            X##_f0 = (X##_f1 << (_FP_W_TYPE_SIZE - (N)) | X##_f0 >> (N)); \
00081            X##_f1 >>= (N);                                       \
00082          })                                                      \
00083        : ({                                                      \
00084            S = ((((N) == _FP_W_TYPE_SIZE                         \
00085                  ? 0                                             \
00086                  : (X##_f1 << (2*_FP_W_TYPE_SIZE - (N))))               \
00087                 | X##_f0) != 0);                                 \
00088            X##_f0 = (X##_f1 >> ((N) - _FP_W_TYPE_SIZE));                \
00089            X##_f1 = 0;                                                  \
00090          }))
00091 
00092 #define _FP_FRAC_SRS_2(X,N,sz)                                          \
00093 (void)(((N) < _FP_W_TYPE_SIZE)                                          \
00094        ? ({                                                      \
00095            X##_f0 = (X##_f1 << (_FP_W_TYPE_SIZE - (N)) | X##_f0 >> (N) | \
00096                     (__builtin_constant_p(N) && (N) == 1                \
00097                      ? X##_f0 & 1                                \
00098                      : (X##_f0 << (_FP_W_TYPE_SIZE - (N))) != 0));      \
00099            X##_f1 >>= (N);                                       \
00100          })                                                      \
00101        : ({                                                      \
00102            X##_f0 = (X##_f1 >> ((N) - _FP_W_TYPE_SIZE) |                \
00103                     ((((N) == _FP_W_TYPE_SIZE                           \
00104                       ? 0                                        \
00105                       : (X##_f1 << (2*_FP_W_TYPE_SIZE - (N))))   \
00106                      | X##_f0) != 0));                           \
00107            X##_f1 = 0;                                                  \
00108          }))
00109 
00110 #define _FP_FRAC_ADDI_2(X,I)       \
00111   __FP_FRAC_ADDI_2(X##_f1, X##_f0, I)
00112 
00113 #define _FP_FRAC_ADD_2(R,X,Y)      \
00114   __FP_FRAC_ADD_2(R##_f1, R##_f0, X##_f1, X##_f0, Y##_f1, Y##_f0)
00115 
00116 #define _FP_FRAC_SUB_2(R,X,Y)      \
00117   __FP_FRAC_SUB_2(R##_f1, R##_f0, X##_f1, X##_f0, Y##_f1, Y##_f0)
00118 
00119 #define _FP_FRAC_DEC_2(X,Y) \
00120   __FP_FRAC_DEC_2(X##_f1, X##_f0, Y##_f1, Y##_f0)
00121 
00122 #define _FP_FRAC_CLZ_2(R,X) \
00123   do {                      \
00124     if (X##_f1)                    \
00125       __FP_CLZ(R,X##_f1);   \
00126     else                    \
00127     {                       \
00128       __FP_CLZ(R,X##_f0);   \
00129       R += _FP_W_TYPE_SIZE; \
00130     }                       \
00131   } while(0)
00132 
00133 /* Predicates */
00134 #define _FP_FRAC_NEGP_2(X)  ((_FP_WS_TYPE)X##_f1 < 0)
00135 #define _FP_FRAC_ZEROP_2(X) ((X##_f1 | X##_f0) == 0)
00136 #define _FP_FRAC_OVERP_2(fs,X)     (_FP_FRAC_HIGH_##fs(X) & _FP_OVERFLOW_##fs)
00137 #define _FP_FRAC_CLEAR_OVERP_2(fs,X)      (_FP_FRAC_HIGH_##fs(X) &= ~_FP_OVERFLOW_##fs)
00138 #define _FP_FRAC_EQ_2(X, Y) (X##_f1 == Y##_f1 && X##_f0 == Y##_f0)
00139 #define _FP_FRAC_GT_2(X, Y) \
00140   (X##_f1 > Y##_f1 || (X##_f1 == Y##_f1 && X##_f0 > Y##_f0))
00141 #define _FP_FRAC_GE_2(X, Y) \
00142   (X##_f1 > Y##_f1 || (X##_f1 == Y##_f1 && X##_f0 >= Y##_f0))
00143 
00144 #define _FP_ZEROFRAC_2             0, 0
00145 #define _FP_MINFRAC_2              0, 1
00146 #define _FP_MAXFRAC_2              (~(_FP_WS_TYPE)0), (~(_FP_WS_TYPE)0)
00147 
00148 /*
00149  * Internals 
00150  */
00151 
00152 #define __FP_FRAC_SET_2(X,I1,I0)   (X##_f0 = I0, X##_f1 = I1)
00153 
00154 #define __FP_CLZ_2(R, xh, xl)      \
00155   do {                      \
00156     if (xh)                 \
00157       __FP_CLZ(R,xh);              \
00158     else                    \
00159     {                       \
00160       __FP_CLZ(R,xl);              \
00161       R += _FP_W_TYPE_SIZE; \
00162     }                       \
00163   } while(0)
00164 
00165 #if 0
00166 
00167 #ifndef __FP_FRAC_ADDI_2
00168 #define __FP_FRAC_ADDI_2(xh, xl, i)       \
00169   (xh += ((xl += i) < i))
00170 #endif
00171 #ifndef __FP_FRAC_ADD_2
00172 #define __FP_FRAC_ADD_2(rh, rl, xh, xl, yh, yl)  \
00173   (rh = xh + yh + ((rl = xl + yl) < xl))
00174 #endif
00175 #ifndef __FP_FRAC_SUB_2
00176 #define __FP_FRAC_SUB_2(rh, rl, xh, xl, yh, yl)  \
00177   (rh = xh - yh - ((rl = xl - yl) > xl))
00178 #endif
00179 #ifndef __FP_FRAC_DEC_2
00180 #define __FP_FRAC_DEC_2(xh, xl, yh, yl)   \
00181   do {                             \
00182     UWtype _t = xl;                \
00183     xh -= yh + ((xl -= yl) > _t);  \
00184   } while (0)
00185 #endif
00186 
00187 #else
00188 
00189 #undef __FP_FRAC_ADDI_2
00190 #define __FP_FRAC_ADDI_2(xh, xl, i)       add_ssaaaa(xh, xl, xh, xl, 0, i)
00191 #undef __FP_FRAC_ADD_2
00192 #define __FP_FRAC_ADD_2                   add_ssaaaa
00193 #undef __FP_FRAC_SUB_2
00194 #define __FP_FRAC_SUB_2                   sub_ddmmss
00195 #undef __FP_FRAC_DEC_2
00196 #define __FP_FRAC_DEC_2(xh, xl, yh, yl)   sub_ddmmss(xh, xl, xh, xl, yh, yl)
00197 
00198 #endif
00199 
00200 /*
00201  * Unpack the raw bits of a native fp value.  Do not classify or
00202  * normalize the data.
00203  */
00204 
00205 #define _FP_UNPACK_RAW_2(fs, X, val)                    \
00206   do {                                           \
00207     union _FP_UNION_##fs _flo; _flo.flt = (val); \
00208                                                  \
00209     X##_f0 = _flo.bits.frac0;                           \
00210     X##_f1 = _flo.bits.frac1;                           \
00211     X##_e  = _flo.bits.exp;                      \
00212     X##_s  = _flo.bits.sign;                            \
00213   } while (0)
00214 
00215 #define _FP_UNPACK_RAW_2_P(fs, X, val)                  \
00216   do {                                           \
00217     union _FP_UNION_##fs *_flo =                 \
00218       (union _FP_UNION_##fs *)(val);                    \
00219                                                  \
00220     X##_f0 = _flo->bits.frac0;                          \
00221     X##_f1 = _flo->bits.frac1;                          \
00222     X##_e  = _flo->bits.exp;                            \
00223     X##_s  = _flo->bits.sign;                           \
00224   } while (0)
00225 
00226 
00227 /*
00228  * Repack the raw bits of a native fp value.
00229  */
00230 
00231 #define _FP_PACK_RAW_2(fs, val, X)               \
00232   do {                                           \
00233     union _FP_UNION_##fs _flo;                          \
00234                                                  \
00235     _flo.bits.frac0 = X##_f0;                           \
00236     _flo.bits.frac1 = X##_f1;                           \
00237     _flo.bits.exp   = X##_e;                            \
00238     _flo.bits.sign  = X##_s;                            \
00239                                                  \
00240     (val) = _flo.flt;                                   \
00241   } while (0)
00242 
00243 #define _FP_PACK_RAW_2_P(fs, val, X)                    \
00244   do {                                           \
00245     union _FP_UNION_##fs *_flo =                 \
00246       (union _FP_UNION_##fs *)(val);                    \
00247                                                  \
00248     _flo->bits.frac0 = X##_f0;                          \
00249     _flo->bits.frac1 = X##_f1;                          \
00250     _flo->bits.exp   = X##_e;                           \
00251     _flo->bits.sign  = X##_s;                           \
00252   } while (0)
00253 
00254 
00255 /*
00256  * Multiplication algorithms:
00257  */
00258 
00259 /* Given a 1W * 1W => 2W primitive, do the extended multiplication.  */
00260 
00261 #define _FP_MUL_MEAT_2_wide(wfracbits, R, X, Y, doit)                 \
00262   do {                                                         \
00263     _FP_FRAC_DECL_4(_z); _FP_FRAC_DECL_2(_b); _FP_FRAC_DECL_2(_c);    \
00264                                                                \
00265     doit(_FP_FRAC_WORD_4(_z,1), _FP_FRAC_WORD_4(_z,0), X##_f0, Y##_f0);      \
00266     doit(_b_f1, _b_f0, X##_f0, Y##_f1);                               \
00267     doit(_c_f1, _c_f0, X##_f1, Y##_f0);                               \
00268     doit(_FP_FRAC_WORD_4(_z,3), _FP_FRAC_WORD_4(_z,2), X##_f1, Y##_f1);      \
00269                                                                \
00270     __FP_FRAC_ADD_3(_FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2),      \
00271                   _FP_FRAC_WORD_4(_z,1), 0, _b_f1, _b_f0,             \
00272                   _FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2), \
00273                   _FP_FRAC_WORD_4(_z,1));                      \
00274     __FP_FRAC_ADD_3(_FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2),      \
00275                   _FP_FRAC_WORD_4(_z,1), 0, _c_f1, _c_f0,             \
00276                   _FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2), \
00277                   _FP_FRAC_WORD_4(_z,1));                      \
00278                                                                \
00279     /* Normalize since we know where the msb of the multiplicands     \
00280        were (bit B), we know that the msb of the of the product is    \
00281        at either 2B or 2B-1.  */                               \
00282     _FP_FRAC_SRS_4(_z, wfracbits-1, 2*wfracbits);                     \
00283     R##_f0 = _FP_FRAC_WORD_4(_z,0);                                   \
00284     R##_f1 = _FP_FRAC_WORD_4(_z,1);                                   \
00285   } while (0)
00286 
00287 /* Given a 1W * 1W => 2W primitive, do the extended multiplication.
00288    Do only 3 multiplications instead of four. This one is for machines
00289    where multiplication is much more expensive than subtraction.  */
00290 
00291 #define _FP_MUL_MEAT_2_wide_3mul(wfracbits, R, X, Y, doit)            \
00292   do {                                                         \
00293     _FP_FRAC_DECL_4(_z); _FP_FRAC_DECL_2(_b); _FP_FRAC_DECL_2(_c);    \
00294     _FP_W_TYPE _d;                                             \
00295     int _c1, _c2;                                              \
00296                                                                \
00297     _b_f0 = X##_f0 + X##_f1;                                          \
00298     _c1 = _b_f0 < X##_f0;                                      \
00299     _b_f1 = Y##_f0 + Y##_f1;                                          \
00300     _c2 = _b_f1 < Y##_f0;                                      \
00301     doit(_d, _FP_FRAC_WORD_4(_z,0), X##_f0, Y##_f0);                  \
00302     doit(_FP_FRAC_WORD_4(_z,2), _FP_FRAC_WORD_4(_z,1), _b_f0, _b_f1); \
00303     doit(_c_f1, _c_f0, X##_f1, Y##_f1);                               \
00304                                                                \
00305     _b_f0 &= -_c2;                                             \
00306     _b_f1 &= -_c1;                                             \
00307     __FP_FRAC_ADD_3(_FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2),      \
00308                   _FP_FRAC_WORD_4(_z,1), (_c1 & _c2), 0, _d,          \
00309                   0, _FP_FRAC_WORD_4(_z,2), _FP_FRAC_WORD_4(_z,1));   \
00310     __FP_FRAC_ADDI_2(_FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2),     \
00311                    _b_f0);                                     \
00312     __FP_FRAC_ADDI_2(_FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2),     \
00313                    _b_f1);                                     \
00314     __FP_FRAC_DEC_3(_FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2),      \
00315                   _FP_FRAC_WORD_4(_z,1),                       \
00316                   0, _d, _FP_FRAC_WORD_4(_z,0));               \
00317     __FP_FRAC_DEC_3(_FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2),      \
00318                   _FP_FRAC_WORD_4(_z,1), 0, _c_f1, _c_f0);            \
00319     __FP_FRAC_ADD_2(_FP_FRAC_WORD_4(_z,3), _FP_FRAC_WORD_4(_z,2),     \
00320                   _c_f1, _c_f0,                                \
00321                   _FP_FRAC_WORD_4(_z,3), _FP_FRAC_WORD_4(_z,2));      \
00322                                                                \
00323     /* Normalize since we know where the msb of the multiplicands     \
00324        were (bit B), we know that the msb of the of the product is    \
00325        at either 2B or 2B-1.  */                               \
00326     _FP_FRAC_SRS_4(_z, wfracbits-1, 2*wfracbits);                     \
00327     R##_f0 = _FP_FRAC_WORD_4(_z,0);                                   \
00328     R##_f1 = _FP_FRAC_WORD_4(_z,1);                                   \
00329   } while (0)
00330 
00331 #define _FP_MUL_MEAT_2_gmp(wfracbits, R, X, Y)                        \
00332   do {                                                         \
00333     _FP_FRAC_DECL_4(_z);                                       \
00334     _FP_W_TYPE _x[2], _y[2];                                          \
00335     _x[0] = X##_f0; _x[1] = X##_f1;                                   \
00336     _y[0] = Y##_f0; _y[1] = Y##_f1;                                   \
00337                                                                \
00338     mpn_mul_n(_z_f, _x, _y, 2);                                       \
00339                                                                \
00340     /* Normalize since we know where the msb of the multiplicands     \
00341        were (bit B), we know that the msb of the of the product is    \
00342        at either 2B or 2B-1.  */                               \
00343     _FP_FRAC_SRS_4(_z, wfracbits-1, 2*wfracbits);                     \
00344     R##_f0 = _z_f[0];                                                 \
00345     R##_f1 = _z_f[1];                                                 \
00346   } while (0)
00347 
00348 /* Do at most 120x120=240 bits multiplication using double floating
00349    point multiplication.  This is useful if floating point
00350    multiplication has much bigger throughput than integer multiply.
00351    It is supposed to work for _FP_W_TYPE_SIZE 64 and wfracbits
00352    between 106 and 120 only.  
00353    Caller guarantees that X and Y has (1LLL << (wfracbits - 1)) set.
00354    SETFETZ is a macro which will disable all FPU exceptions and set rounding
00355    towards zero,  RESETFE should optionally reset it back.  */
00356 
00357 #define _FP_MUL_MEAT_2_120_240_double(wfracbits, R, X, Y, setfetz, resetfe)  \
00358   do {                                                                \
00359     static const double _const[] = {                                         \
00360       /* 2^-24 */ 5.9604644775390625e-08,                             \
00361       /* 2^-48 */ 3.5527136788005009e-15,                             \
00362       /* 2^-72 */ 2.1175823681357508e-22,                             \
00363       /* 2^-96 */ 1.2621774483536189e-29,                             \
00364       /* 2^28 */ 2.68435456e+08,                                      \
00365       /* 2^4 */ 1.600000e+01,                                                \
00366       /* 2^-20 */ 9.5367431640625e-07,                                       \
00367       /* 2^-44 */ 5.6843418860808015e-14,                             \
00368       /* 2^-68 */ 3.3881317890172014e-21,                             \
00369       /* 2^-92 */ 2.0194839173657902e-28,                             \
00370       /* 2^-116 */ 1.2037062152420224e-35};                                  \
00371     double _a240, _b240, _c240, _d240, _e240, _f240,                         \
00372           _g240, _h240, _i240, _j240, _k240;                                 \
00373     union { double d; UDItype i; } _l240, _m240, _n240, _o240,               \
00374                                _p240, _q240, _r240, _s240;                   \
00375     UDItype _t240, _u240, _v240, _w240, _x240, _y240 = 0;                    \
00376                                                                       \
00377     if (wfracbits < 106 || wfracbits > 120)                                  \
00378       abort();                                                               \
00379                                                                       \
00380     setfetz;                                                          \
00381                                                                       \
00382     _e240 = (double)(long)(X##_f0 & 0xffffff);                               \
00383     _j240 = (double)(long)(Y##_f0 & 0xffffff);                               \
00384     _d240 = (double)(long)((X##_f0 >> 24) & 0xffffff);                       \
00385     _i240 = (double)(long)((Y##_f0 >> 24) & 0xffffff);                       \
00386     _c240 = (double)(long)(((X##_f1 << 16) & 0xffffff) | (X##_f0 >> 48));    \
00387     _h240 = (double)(long)(((Y##_f1 << 16) & 0xffffff) | (Y##_f0 >> 48));    \
00388     _b240 = (double)(long)((X##_f1 >> 8) & 0xffffff);                        \
00389     _g240 = (double)(long)((Y##_f1 >> 8) & 0xffffff);                        \
00390     _a240 = (double)(long)(X##_f1 >> 32);                             \
00391     _f240 = (double)(long)(Y##_f1 >> 32);                             \
00392     _e240 *= _const[3];                                                      \
00393     _j240 *= _const[3];                                                      \
00394     _d240 *= _const[2];                                                      \
00395     _i240 *= _const[2];                                                      \
00396     _c240 *= _const[1];                                                      \
00397     _h240 *= _const[1];                                                      \
00398     _b240 *= _const[0];                                                      \
00399     _g240 *= _const[0];                                                      \
00400     _s240.d =                                                 _e240*_j240;\
00401     _r240.d =                                    _d240*_j240 + _e240*_i240;\
00402     _q240.d =                        _c240*_j240 + _d240*_i240 + _e240*_h240;\
00403     _p240.d =            _b240*_j240 + _c240*_i240 + _d240*_h240 + _e240*_g240;\
00404     _o240.d = _a240*_j240 + _b240*_i240 + _c240*_h240 + _d240*_g240 + _e240*_f240;\
00405     _n240.d = _a240*_i240 + _b240*_h240 + _c240*_g240 + _d240*_f240;         \
00406     _m240.d = _a240*_h240 + _b240*_g240 + _c240*_f240;                       \
00407     _l240.d = _a240*_g240 + _b240*_f240;                              \
00408     _k240 =   _a240*_f240;                                            \
00409     _r240.d += _s240.d;                                                      \
00410     _q240.d += _r240.d;                                                      \
00411     _p240.d += _q240.d;                                                      \
00412     _o240.d += _p240.d;                                                      \
00413     _n240.d += _o240.d;                                                      \
00414     _m240.d += _n240.d;                                                      \
00415     _l240.d += _m240.d;                                                      \
00416     _k240 += _l240.d;                                                        \
00417     _s240.d -= ((_const[10]+_s240.d)-_const[10]);                            \
00418     _r240.d -= ((_const[9]+_r240.d)-_const[9]);                              \
00419     _q240.d -= ((_const[8]+_q240.d)-_const[8]);                              \
00420     _p240.d -= ((_const[7]+_p240.d)-_const[7]);                              \
00421     _o240.d += _const[7];                                             \
00422     _n240.d += _const[6];                                             \
00423     _m240.d += _const[5];                                             \
00424     _l240.d += _const[4];                                             \
00425     if (_s240.d != 0.0) _y240 = 1;                                    \
00426     if (_r240.d != 0.0) _y240 = 1;                                    \
00427     if (_q240.d != 0.0) _y240 = 1;                                    \
00428     if (_p240.d != 0.0) _y240 = 1;                                    \
00429     _t240 = (DItype)_k240;                                            \
00430     _u240 = _l240.i;                                                  \
00431     _v240 = _m240.i;                                                  \
00432     _w240 = _n240.i;                                                  \
00433     _x240 = _o240.i;                                                  \
00434     R##_f1 = (_t240 << (128 - (wfracbits - 1)))                              \
00435             | ((_u240 & 0xffffff) >> ((wfracbits - 1) - 104));               \
00436     R##_f0 = ((_u240 & 0xffffff) << (168 - (wfracbits - 1)))                 \
00437             | ((_v240 & 0xffffff) << (144 - (wfracbits - 1)))                \
00438             | ((_w240 & 0xffffff) << (120 - (wfracbits - 1)))                \
00439             | ((_x240 & 0xffffff) >> ((wfracbits - 1) - 96))                 \
00440             | _y240;                                                  \
00441     resetfe;                                                          \
00442   } while (0)
00443 
00444 /*
00445  * Division algorithms:
00446  */
00447 
00448 #define _FP_DIV_MEAT_2_udiv(fs, R, X, Y)                       \
00449   do {                                                         \
00450     _FP_W_TYPE _n_f2, _n_f1, _n_f0, _r_f1, _r_f0, _m_f1, _m_f0;              \
00451     if (_FP_FRAC_GT_2(X, Y))                                          \
00452       {                                                               \
00453        _n_f2 = X##_f1 >> 1;                                    \
00454        _n_f1 = X##_f1 << (_FP_W_TYPE_SIZE - 1) | X##_f0 >> 1;         \
00455        _n_f0 = X##_f0 << (_FP_W_TYPE_SIZE - 1);                \
00456       }                                                               \
00457     else                                                       \
00458       {                                                               \
00459        R##_e--;                                                \
00460        _n_f2 = X##_f1;                                                \
00461        _n_f1 = X##_f0;                                                \
00462        _n_f0 = 0;                                              \
00463       }                                                               \
00464                                                                \
00465     /* Normalize, i.e. make the most significant bit of the           \
00466        denominator set. */                                     \
00467     _FP_FRAC_SLL_2(Y, _FP_WFRACXBITS_##fs);                           \
00468                                                                \
00469     udiv_qrnnd(R##_f1, _r_f1, _n_f2, _n_f1, Y##_f1);                  \
00470     umul_ppmm(_m_f1, _m_f0, R##_f1, Y##_f0);                          \
00471     _r_f0 = _n_f0;                                             \
00472     if (_FP_FRAC_GT_2(_m, _r))                                        \
00473       {                                                               \
00474        R##_f1--;                                               \
00475        _FP_FRAC_ADD_2(_r, Y, _r);                              \
00476        if (_FP_FRAC_GE_2(_r, Y) && _FP_FRAC_GT_2(_m, _r))             \
00477          {                                                     \
00478            R##_f1--;                                           \
00479            _FP_FRAC_ADD_2(_r, Y, _r);                                 \
00480          }                                                     \
00481       }                                                               \
00482     _FP_FRAC_DEC_2(_r, _m);                                    \
00483                                                                \
00484     if (_r_f1 == Y##_f1)                                       \
00485       {                                                               \
00486        /* This is a special case, not an optimization                 \
00487           (_r/Y##_f1 would not fit into UWtype).               \
00488           As _r is guaranteed to be < Y,  R##_f0 can be either        \
00489           (UWtype)-1 or (UWtype)-2.  But as we know what kind         \
00490           of bits it is (sticky, guard, round),  we don't care.       \
00491           We also don't care what the reminder is,  because the       \
00492           guard bit will be set anyway.  -jj */                \
00493        R##_f0 = -1;                                            \
00494       }                                                               \
00495     else                                                       \
00496       {                                                               \
00497        udiv_qrnnd(R##_f0, _r_f1, _r_f1, _r_f0, Y##_f1);        \
00498        umul_ppmm(_m_f1, _m_f0, R##_f0, Y##_f0);                \
00499        _r_f0 = 0;                                              \
00500        if (_FP_FRAC_GT_2(_m, _r))                              \
00501          {                                                     \
00502            R##_f0--;                                           \
00503            _FP_FRAC_ADD_2(_r, Y, _r);                                 \
00504            if (_FP_FRAC_GE_2(_r, Y) && _FP_FRAC_GT_2(_m, _r))         \
00505              {                                                        \
00506               R##_f0--;                                        \
00507               _FP_FRAC_ADD_2(_r, Y, _r);                       \
00508              }                                                        \
00509          }                                                     \
00510        if (!_FP_FRAC_EQ_2(_r, _m))                             \
00511          R##_f0 |= _FP_WORK_STICKY;                                   \
00512       }                                                               \
00513   } while (0)
00514 
00515 
00516 #define _FP_DIV_MEAT_2_gmp(fs, R, X, Y)                               \
00517   do {                                                         \
00518     _FP_W_TYPE _x[4], _y[2], _z[4];                                   \
00519     _y[0] = Y##_f0; _y[1] = Y##_f1;                                   \
00520     _x[0] = _x[3] = 0;                                                \
00521     if (_FP_FRAC_GT_2(X, Y))                                          \
00522       {                                                               \
00523        R##_e++;                                                \
00524        _x[1] = (X##_f0 << (_FP_WFRACBITS_##fs-1 - _FP_W_TYPE_SIZE) |  \
00525                X##_f1 >> (_FP_W_TYPE_SIZE -                           \
00526                          (_FP_WFRACBITS_##fs-1 - _FP_W_TYPE_SIZE)));  \
00527        _x[2] = X##_f1 << (_FP_WFRACBITS_##fs-1 - _FP_W_TYPE_SIZE);    \
00528       }                                                               \
00529     else                                                       \
00530       {                                                               \
00531        _x[1] = (X##_f0 << (_FP_WFRACBITS_##fs - _FP_W_TYPE_SIZE) |    \
00532                X##_f1 >> (_FP_W_TYPE_SIZE -                           \
00533                          (_FP_WFRACBITS_##fs - _FP_W_TYPE_SIZE)));    \
00534        _x[2] = X##_f1 << (_FP_WFRACBITS_##fs - _FP_W_TYPE_SIZE);      \
00535       }                                                               \
00536                                                                \
00537     (void) mpn_divrem (_z, 0, _x, 4, _y, 2);                          \
00538     R##_f1 = _z[1];                                            \
00539     R##_f0 = _z[0] | ((_x[0] | _x[1]) != 0);                          \
00540   } while (0)
00541 
00542 
00543 /*
00544  * Square root algorithms:
00545  * We have just one right now, maybe Newton approximation
00546  * should be added for those machines where division is fast.
00547  */
00548  
00549 #define _FP_SQRT_MEAT_2(R, S, T, X, q)                  \
00550   do {                                           \
00551     while (q)                                    \
00552       {                                                 \
00553        T##_f1 = S##_f1 + q;                      \
00554        if (T##_f1 <= X##_f1)                            \
00555          {                                       \
00556            S##_f1 = T##_f1 + q;                  \
00557            X##_f1 -= T##_f1;                            \
00558            R##_f1 += q;                          \
00559          }                                       \
00560        _FP_FRAC_SLL_2(X, 1);                            \
00561        q >>= 1;                                  \
00562       }                                                 \
00563     q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1);         \
00564     while (q != _FP_WORK_ROUND)                         \
00565       {                                                 \
00566        T##_f0 = S##_f0 + q;                      \
00567        T##_f1 = S##_f1;                          \
00568        if (T##_f1 < X##_f1 ||                           \
00569            (T##_f1 == X##_f1 && T##_f0 <= X##_f0))      \
00570          {                                       \
00571            S##_f0 = T##_f0 + q;                  \
00572            S##_f1 += (T##_f0 > S##_f0);          \
00573            _FP_FRAC_DEC_2(X, T);                 \
00574            R##_f0 += q;                          \
00575          }                                       \
00576        _FP_FRAC_SLL_2(X, 1);                            \
00577        q >>= 1;                                  \
00578       }                                                 \
00579     if (X##_f0 | X##_f1)                         \
00580       {                                                 \
00581        if (S##_f1 < X##_f1 ||                           \
00582            (S##_f1 == X##_f1 && S##_f0 < X##_f0))       \
00583          R##_f0 |= _FP_WORK_ROUND;               \
00584        R##_f0 |= _FP_WORK_STICKY;                \
00585       }                                                 \
00586   } while (0)
00587 
00588 
00589 /*
00590  * Assembly/disassembly for converting to/from integral types.  
00591  * No shifting or overflow handled here.
00592  */
00593 
00594 #define _FP_FRAC_ASSEMBLE_2(r, X, rsize)  \
00595 (void)((rsize <= _FP_W_TYPE_SIZE)         \
00596        ? ({ r = X##_f0; })                \
00597        : ({                               \
00598            r = X##_f1;                           \
00599            r <<= _FP_W_TYPE_SIZE;         \
00600            r += X##_f0;                   \
00601          }))
00602 
00603 #define _FP_FRAC_DISASSEMBLE_2(X, r, rsize)                           \
00604   do {                                                         \
00605     X##_f0 = r;                                                       \
00606     X##_f1 = (rsize <= _FP_W_TYPE_SIZE ? 0 : r >> _FP_W_TYPE_SIZE);   \
00607   } while (0)
00608 
00609 /*
00610  * Convert FP values between word sizes
00611  */
00612 
00613 #define _FP_FRAC_COPY_1_2(D, S)           (D##_f = S##_f0)
00614 
00615 #define _FP_FRAC_COPY_2_1(D, S)           ((D##_f0 = S##_f), (D##_f1 = 0))
00616 
00617 #define _FP_FRAC_COPY_2_2(D,S)            _FP_FRAC_COPY_2(D,S)