Back to home page

OSCL-LXR

 
 

    


0001 /* Software floating-point emulation.
0002    Basic four-word fraction declaration and manipulation.
0003    Copyright (C) 1997,1998,1999 Free Software Foundation, Inc.
0004    This file is part of the GNU C Library.
0005    Contributed by Richard Henderson (rth@cygnus.com),
0006           Jakub Jelinek (jj@ultra.linux.cz),
0007           David S. Miller (davem@redhat.com) and
0008           Peter Maydell (pmaydell@chiark.greenend.org.uk).
0009 
0010    The GNU C Library is free software; you can redistribute it and/or
0011    modify it under the terms of the GNU Library General Public License as
0012    published by the Free Software Foundation; either version 2 of the
0013    License, or (at your option) any later version.
0014 
0015    The GNU C Library is distributed in the hope that it will be useful,
0016    but WITHOUT ANY WARRANTY; without even the implied warranty of
0017    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
0018    Library General Public License for more details.
0019 
0020    You should have received a copy of the GNU Library General Public
0021    License along with the GNU C Library; see the file COPYING.LIB.  If
0022    not, write to the Free Software Foundation, Inc.,
0023    59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  */
0024 
0025 #ifndef __MATH_EMU_OP_4_H__
0026 #define __MATH_EMU_OP_4_H__
0027 
0028 #define _FP_FRAC_DECL_4(X)  _FP_W_TYPE X##_f[4]
0029 #define _FP_FRAC_COPY_4(D,S)            \
0030   (D##_f[0] = S##_f[0], D##_f[1] = S##_f[1],    \
0031    D##_f[2] = S##_f[2], D##_f[3] = S##_f[3])
0032 #define _FP_FRAC_SET_4(X,I) __FP_FRAC_SET_4(X, I)
0033 #define _FP_FRAC_HIGH_4(X)  (X##_f[3])
0034 #define _FP_FRAC_LOW_4(X)   (X##_f[0])
0035 #define _FP_FRAC_WORD_4(X,w)    (X##_f[w])
0036 
0037 #define _FP_FRAC_SLL_4(X,N)                     \
0038   do {                                  \
0039     _FP_I_TYPE _up, _down, _skip, _i;                   \
0040     _skip = (N) / _FP_W_TYPE_SIZE;                  \
0041     _up = (N) % _FP_W_TYPE_SIZE;                    \
0042     _down = _FP_W_TYPE_SIZE - _up;                  \
0043     if (!_up)                               \
0044       for (_i = 3; _i >= _skip; --_i)                   \
0045     X##_f[_i] = X##_f[_i-_skip];                    \
0046     else                                \
0047       {                                 \
0048     for (_i = 3; _i > _skip; --_i)                  \
0049       X##_f[_i] = X##_f[_i-_skip] << _up                \
0050               | X##_f[_i-_skip-1] >> _down;         \
0051     X##_f[_i--] = X##_f[0] << _up;                  \
0052       }                                 \
0053     for (; _i >= 0; --_i)                       \
0054       X##_f[_i] = 0;                            \
0055   } while (0)
0056 
0057 /* This one was broken too */
0058 #define _FP_FRAC_SRL_4(X,N)                     \
0059   do {                                  \
0060     _FP_I_TYPE _up, _down, _skip, _i;                   \
0061     _skip = (N) / _FP_W_TYPE_SIZE;                  \
0062     _down = (N) % _FP_W_TYPE_SIZE;                  \
0063     _up = _FP_W_TYPE_SIZE - _down;                  \
0064     if (!_down)                             \
0065       for (_i = 0; _i <= 3-_skip; ++_i)                 \
0066     X##_f[_i] = X##_f[_i+_skip];                    \
0067     else                                \
0068       {                                 \
0069     for (_i = 0; _i < 3-_skip; ++_i)                \
0070       X##_f[_i] = X##_f[_i+_skip] >> _down              \
0071               | X##_f[_i+_skip+1] << _up;           \
0072     X##_f[_i++] = X##_f[3] >> _down;                \
0073       }                                 \
0074     for (; _i < 4; ++_i)                        \
0075       X##_f[_i] = 0;                            \
0076   } while (0)
0077 
0078 
0079 /* Right shift with sticky-lsb. 
0080  * What this actually means is that we do a standard right-shift,
0081  * but that if any of the bits that fall off the right hand side
0082  * were one then we always set the LSbit.
0083  */
0084 #define _FP_FRAC_SRS_4(X,N,size)                    \
0085   do {                                  \
0086     _FP_I_TYPE _up, _down, _skip, _i;                   \
0087     _FP_W_TYPE _s;                          \
0088     _skip = (N) / _FP_W_TYPE_SIZE;                  \
0089     _down = (N) % _FP_W_TYPE_SIZE;                  \
0090     _up = _FP_W_TYPE_SIZE - _down;                  \
0091     for (_s = _i = 0; _i < _skip; ++_i)                 \
0092       _s |= X##_f[_i];                          \
0093     _s |= X##_f[_i] << _up;                     \
0094 /* s is now != 0 if we want to set the LSbit */             \
0095     if (!_down)                             \
0096       for (_i = 0; _i <= 3-_skip; ++_i)                 \
0097     X##_f[_i] = X##_f[_i+_skip];                    \
0098     else                                \
0099       {                                 \
0100     for (_i = 0; _i < 3-_skip; ++_i)                \
0101       X##_f[_i] = X##_f[_i+_skip] >> _down              \
0102               | X##_f[_i+_skip+1] << _up;           \
0103     X##_f[_i++] = X##_f[3] >> _down;                \
0104       }                                 \
0105     for (; _i < 4; ++_i)                        \
0106       X##_f[_i] = 0;                            \
0107     /* don't fix the LSB until the very end when we're sure f[0] is stable */   \
0108     X##_f[0] |= (_s != 0);                      \
0109   } while (0)
0110 
0111 #define _FP_FRAC_ADD_4(R,X,Y)                       \
0112   __FP_FRAC_ADD_4(R##_f[3], R##_f[2], R##_f[1], R##_f[0],       \
0113           X##_f[3], X##_f[2], X##_f[1], X##_f[0],       \
0114           Y##_f[3], Y##_f[2], Y##_f[1], Y##_f[0])
0115 
0116 #define _FP_FRAC_SUB_4(R,X,Y)                       \
0117   __FP_FRAC_SUB_4(R##_f[3], R##_f[2], R##_f[1], R##_f[0],       \
0118           X##_f[3], X##_f[2], X##_f[1], X##_f[0],       \
0119           Y##_f[3], Y##_f[2], Y##_f[1], Y##_f[0])
0120 
0121 #define _FP_FRAC_DEC_4(X,Y)                     \
0122   __FP_FRAC_DEC_4(X##_f[3], X##_f[2], X##_f[1], X##_f[0],       \
0123           Y##_f[3], Y##_f[2], Y##_f[1], Y##_f[0])
0124 
0125 #define _FP_FRAC_ADDI_4(X,I)                        \
0126   __FP_FRAC_ADDI_4(X##_f[3], X##_f[2], X##_f[1], X##_f[0], I)
0127 
0128 #define _FP_ZEROFRAC_4  0,0,0,0
0129 #define _FP_MINFRAC_4   0,0,0,1
0130 #define _FP_MAXFRAC_4   (~(_FP_WS_TYPE)0), (~(_FP_WS_TYPE)0), (~(_FP_WS_TYPE)0), (~(_FP_WS_TYPE)0)
0131 
0132 #define _FP_FRAC_ZEROP_4(X)     ((X##_f[0] | X##_f[1] | X##_f[2] | X##_f[3]) == 0)
0133 #define _FP_FRAC_NEGP_4(X)      ((_FP_WS_TYPE)X##_f[3] < 0)
0134 #define _FP_FRAC_OVERP_4(fs,X)  (_FP_FRAC_HIGH_##fs(X) & _FP_OVERFLOW_##fs)
0135 #define _FP_FRAC_CLEAR_OVERP_4(fs,X)  (_FP_FRAC_HIGH_##fs(X) &= ~_FP_OVERFLOW_##fs)
0136 
0137 #define _FP_FRAC_EQ_4(X,Y)              \
0138  (X##_f[0] == Y##_f[0] && X##_f[1] == Y##_f[1]      \
0139   && X##_f[2] == Y##_f[2] && X##_f[3] == Y##_f[3])
0140 
0141 #define _FP_FRAC_GT_4(X,Y)              \
0142  (X##_f[3] > Y##_f[3] ||                \
0143   (X##_f[3] == Y##_f[3] && (X##_f[2] > Y##_f[2] ||  \
0144    (X##_f[2] == Y##_f[2] && (X##_f[1] > Y##_f[1] || \
0145     (X##_f[1] == Y##_f[1] && X##_f[0] > Y##_f[0])   \
0146    ))                           \
0147   ))                            \
0148  )
0149 
0150 #define _FP_FRAC_GE_4(X,Y)              \
0151  (X##_f[3] > Y##_f[3] ||                \
0152   (X##_f[3] == Y##_f[3] && (X##_f[2] > Y##_f[2] ||  \
0153    (X##_f[2] == Y##_f[2] && (X##_f[1] > Y##_f[1] || \
0154     (X##_f[1] == Y##_f[1] && X##_f[0] >= Y##_f[0])  \
0155    ))                           \
0156   ))                            \
0157  )
0158 
0159 
0160 #define _FP_FRAC_CLZ_4(R,X)     \
0161   do {                  \
0162     if (X##_f[3])           \
0163     {                   \
0164     __FP_CLZ(R,X##_f[3]);       \
0165     }                   \
0166     else if (X##_f[2])          \
0167     {                   \
0168     __FP_CLZ(R,X##_f[2]);       \
0169     R += _FP_W_TYPE_SIZE;       \
0170     }                   \
0171     else if (X##_f[1])          \
0172     {                   \
0173     __FP_CLZ(R,X##_f[2]);       \
0174     R += _FP_W_TYPE_SIZE*2;     \
0175     }                   \
0176     else                \
0177     {                   \
0178     __FP_CLZ(R,X##_f[0]);       \
0179     R += _FP_W_TYPE_SIZE*3;     \
0180     }                   \
0181   } while(0)
0182 
0183 
0184 #define _FP_UNPACK_RAW_4(fs, X, val)                \
0185   do {                              \
0186     union _FP_UNION_##fs _flo; _flo.flt = (val);        \
0187     X##_f[0] = _flo.bits.frac0;                 \
0188     X##_f[1] = _flo.bits.frac1;                 \
0189     X##_f[2] = _flo.bits.frac2;                 \
0190     X##_f[3] = _flo.bits.frac3;                 \
0191     X##_e  = _flo.bits.exp;                 \
0192     X##_s  = _flo.bits.sign;                    \
0193   } while (0)
0194 
0195 #define _FP_UNPACK_RAW_4_P(fs, X, val)              \
0196   do {                              \
0197     union _FP_UNION_##fs *_flo =                \
0198       (union _FP_UNION_##fs *)(val);                \
0199                                 \
0200     X##_f[0] = _flo->bits.frac0;                \
0201     X##_f[1] = _flo->bits.frac1;                \
0202     X##_f[2] = _flo->bits.frac2;                \
0203     X##_f[3] = _flo->bits.frac3;                \
0204     X##_e  = _flo->bits.exp;                    \
0205     X##_s  = _flo->bits.sign;                   \
0206   } while (0)
0207 
0208 #define _FP_PACK_RAW_4(fs, val, X)              \
0209   do {                              \
0210     union _FP_UNION_##fs _flo;                  \
0211     _flo.bits.frac0 = X##_f[0];                 \
0212     _flo.bits.frac1 = X##_f[1];                 \
0213     _flo.bits.frac2 = X##_f[2];                 \
0214     _flo.bits.frac3 = X##_f[3];                 \
0215     _flo.bits.exp   = X##_e;                    \
0216     _flo.bits.sign  = X##_s;                    \
0217     (val) = _flo.flt;                       \
0218   } while (0)
0219 
0220 #define _FP_PACK_RAW_4_P(fs, val, X)                \
0221   do {                              \
0222     union _FP_UNION_##fs *_flo =                \
0223       (union _FP_UNION_##fs *)(val);                \
0224                                 \
0225     _flo->bits.frac0 = X##_f[0];                \
0226     _flo->bits.frac1 = X##_f[1];                \
0227     _flo->bits.frac2 = X##_f[2];                \
0228     _flo->bits.frac3 = X##_f[3];                \
0229     _flo->bits.exp   = X##_e;                   \
0230     _flo->bits.sign  = X##_s;                   \
0231   } while (0)
0232 
0233 /*
0234  * Multiplication algorithms:
0235  */
0236 
0237 /* Given a 1W * 1W => 2W primitive, do the extended multiplication.  */
0238 
0239 #define _FP_MUL_MEAT_4_wide(wfracbits, R, X, Y, doit)               \
0240   do {                                      \
0241     _FP_FRAC_DECL_8(_z); _FP_FRAC_DECL_2(_b); _FP_FRAC_DECL_2(_c);      \
0242     _FP_FRAC_DECL_2(_d); _FP_FRAC_DECL_2(_e); _FP_FRAC_DECL_2(_f);      \
0243                                         \
0244     doit(_FP_FRAC_WORD_8(_z,1), _FP_FRAC_WORD_8(_z,0), X##_f[0], Y##_f[0]); \
0245     doit(_b_f1, _b_f0, X##_f[0], Y##_f[1]);                 \
0246     doit(_c_f1, _c_f0, X##_f[1], Y##_f[0]);                 \
0247     doit(_d_f1, _d_f0, X##_f[1], Y##_f[1]);                 \
0248     doit(_e_f1, _e_f0, X##_f[0], Y##_f[2]);                 \
0249     doit(_f_f1, _f_f0, X##_f[2], Y##_f[0]);                 \
0250     __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,3),_FP_FRAC_WORD_8(_z,2),        \
0251             _FP_FRAC_WORD_8(_z,1), 0,_b_f1,_b_f0,           \
0252             0,0,_FP_FRAC_WORD_8(_z,1));                 \
0253     __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,3),_FP_FRAC_WORD_8(_z,2),        \
0254             _FP_FRAC_WORD_8(_z,1), 0,_c_f1,_c_f0,           \
0255             _FP_FRAC_WORD_8(_z,3),_FP_FRAC_WORD_8(_z,2),        \
0256             _FP_FRAC_WORD_8(_z,1));                 \
0257     __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,4),_FP_FRAC_WORD_8(_z,3),        \
0258             _FP_FRAC_WORD_8(_z,2), 0,_d_f1,_d_f0,           \
0259             0,_FP_FRAC_WORD_8(_z,3),_FP_FRAC_WORD_8(_z,2));     \
0260     __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,4),_FP_FRAC_WORD_8(_z,3),        \
0261             _FP_FRAC_WORD_8(_z,2), 0,_e_f1,_e_f0,           \
0262             _FP_FRAC_WORD_8(_z,4),_FP_FRAC_WORD_8(_z,3),        \
0263             _FP_FRAC_WORD_8(_z,2));                 \
0264     __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,4),_FP_FRAC_WORD_8(_z,3),        \
0265             _FP_FRAC_WORD_8(_z,2), 0,_f_f1,_f_f0,           \
0266             _FP_FRAC_WORD_8(_z,4),_FP_FRAC_WORD_8(_z,3),        \
0267             _FP_FRAC_WORD_8(_z,2));                 \
0268     doit(_b_f1, _b_f0, X##_f[0], Y##_f[3]);                 \
0269     doit(_c_f1, _c_f0, X##_f[3], Y##_f[0]);                 \
0270     doit(_d_f1, _d_f0, X##_f[1], Y##_f[2]);                 \
0271     doit(_e_f1, _e_f0, X##_f[2], Y##_f[1]);                 \
0272     __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,5),_FP_FRAC_WORD_8(_z,4),        \
0273             _FP_FRAC_WORD_8(_z,3), 0,_b_f1,_b_f0,           \
0274             0,_FP_FRAC_WORD_8(_z,4),_FP_FRAC_WORD_8(_z,3));     \
0275     __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,5),_FP_FRAC_WORD_8(_z,4),        \
0276             _FP_FRAC_WORD_8(_z,3), 0,_c_f1,_c_f0,           \
0277             _FP_FRAC_WORD_8(_z,5),_FP_FRAC_WORD_8(_z,4),        \
0278             _FP_FRAC_WORD_8(_z,3));                 \
0279     __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,5),_FP_FRAC_WORD_8(_z,4),        \
0280             _FP_FRAC_WORD_8(_z,3), 0,_d_f1,_d_f0,           \
0281             _FP_FRAC_WORD_8(_z,5),_FP_FRAC_WORD_8(_z,4),        \
0282             _FP_FRAC_WORD_8(_z,3));                 \
0283     __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,5),_FP_FRAC_WORD_8(_z,4),        \
0284             _FP_FRAC_WORD_8(_z,3), 0,_e_f1,_e_f0,           \
0285             _FP_FRAC_WORD_8(_z,5),_FP_FRAC_WORD_8(_z,4),        \
0286             _FP_FRAC_WORD_8(_z,3));                 \
0287     doit(_b_f1, _b_f0, X##_f[2], Y##_f[2]);                 \
0288     doit(_c_f1, _c_f0, X##_f[1], Y##_f[3]);                 \
0289     doit(_d_f1, _d_f0, X##_f[3], Y##_f[1]);                 \
0290     doit(_e_f1, _e_f0, X##_f[2], Y##_f[3]);                 \
0291     doit(_f_f1, _f_f0, X##_f[3], Y##_f[2]);                 \
0292     __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,6),_FP_FRAC_WORD_8(_z,5),        \
0293             _FP_FRAC_WORD_8(_z,4), 0,_b_f1,_b_f0,           \
0294             0,_FP_FRAC_WORD_8(_z,5),_FP_FRAC_WORD_8(_z,4));     \
0295     __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,6),_FP_FRAC_WORD_8(_z,5),        \
0296             _FP_FRAC_WORD_8(_z,4), 0,_c_f1,_c_f0,           \
0297             _FP_FRAC_WORD_8(_z,6),_FP_FRAC_WORD_8(_z,5),        \
0298             _FP_FRAC_WORD_8(_z,4));                 \
0299     __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,6),_FP_FRAC_WORD_8(_z,5),        \
0300             _FP_FRAC_WORD_8(_z,4), 0,_d_f1,_d_f0,           \
0301             _FP_FRAC_WORD_8(_z,6),_FP_FRAC_WORD_8(_z,5),        \
0302             _FP_FRAC_WORD_8(_z,4));                 \
0303     __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,7),_FP_FRAC_WORD_8(_z,6),        \
0304             _FP_FRAC_WORD_8(_z,5), 0,_e_f1,_e_f0,           \
0305             0,_FP_FRAC_WORD_8(_z,6),_FP_FRAC_WORD_8(_z,5));     \
0306     __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,7),_FP_FRAC_WORD_8(_z,6),        \
0307             _FP_FRAC_WORD_8(_z,5), 0,_f_f1,_f_f0,           \
0308             _FP_FRAC_WORD_8(_z,7),_FP_FRAC_WORD_8(_z,6),        \
0309             _FP_FRAC_WORD_8(_z,5));                 \
0310     doit(_b_f1, _b_f0, X##_f[3], Y##_f[3]);                 \
0311     __FP_FRAC_ADD_2(_FP_FRAC_WORD_8(_z,7),_FP_FRAC_WORD_8(_z,6),        \
0312             _b_f1,_b_f0,                        \
0313             _FP_FRAC_WORD_8(_z,7),_FP_FRAC_WORD_8(_z,6));       \
0314                                         \
0315     /* Normalize since we know where the msb of the multiplicands       \
0316        were (bit B), we know that the msb of the of the product is      \
0317        at either 2B or 2B-1.  */                        \
0318     _FP_FRAC_SRS_8(_z, wfracbits-1, 2*wfracbits);               \
0319     __FP_FRAC_SET_4(R, _FP_FRAC_WORD_8(_z,3), _FP_FRAC_WORD_8(_z,2),        \
0320             _FP_FRAC_WORD_8(_z,1), _FP_FRAC_WORD_8(_z,0));      \
0321   } while (0)
0322 
0323 #define _FP_MUL_MEAT_4_gmp(wfracbits, R, X, Y)                  \
0324   do {                                      \
0325     _FP_FRAC_DECL_8(_z);                            \
0326                                         \
0327     mpn_mul_n(_z_f, _x_f, _y_f, 4);                     \
0328                                         \
0329     /* Normalize since we know where the msb of the multiplicands       \
0330        were (bit B), we know that the msb of the of the product is      \
0331        at either 2B or 2B-1.  */                        \
0332     _FP_FRAC_SRS_8(_z, wfracbits-1, 2*wfracbits);               \
0333     __FP_FRAC_SET_4(R, _FP_FRAC_WORD_8(_z,3), _FP_FRAC_WORD_8(_z,2),        \
0334             _FP_FRAC_WORD_8(_z,1), _FP_FRAC_WORD_8(_z,0));      \
0335   } while (0)
0336 
0337 /*
0338  * Helper utility for _FP_DIV_MEAT_4_udiv:
0339  * pppp = m * nnn
0340  */
0341 #define umul_ppppmnnn(p3,p2,p1,p0,m,n2,n1,n0)                   \
0342   do {                                      \
0343     UWtype _t;                                  \
0344     umul_ppmm(p1,p0,m,n0);                          \
0345     umul_ppmm(p2,_t,m,n1);                          \
0346     __FP_FRAC_ADDI_2(p2,p1,_t);                         \
0347     umul_ppmm(p3,_t,m,n2);                          \
0348     __FP_FRAC_ADDI_2(p3,p2,_t);                         \
0349   } while (0)
0350 
0351 /*
0352  * Division algorithms:
0353  */
0354 
0355 #define _FP_DIV_MEAT_4_udiv(fs, R, X, Y)                    \
0356   do {                                      \
0357     int _i;                                 \
0358     _FP_FRAC_DECL_4(_n); _FP_FRAC_DECL_4(_m);                   \
0359     _FP_FRAC_SET_4(_n, _FP_ZEROFRAC_4);                     \
0360     if (_FP_FRAC_GT_4(X, Y))                            \
0361       {                                     \
0362     _n_f[3] = X##_f[0] << (_FP_W_TYPE_SIZE - 1);                \
0363     _FP_FRAC_SRL_4(X, 1);                           \
0364       }                                     \
0365     else                                    \
0366       R##_e--;                                  \
0367                                         \
0368     /* Normalize, i.e. make the most significant bit of the             \
0369        denominator set. */                          \
0370     _FP_FRAC_SLL_4(Y, _FP_WFRACXBITS_##fs);                 \
0371                                         \
0372     for (_i = 3; ; _i--)                            \
0373       {                                     \
0374         if (X##_f[3] == Y##_f[3])                       \
0375           {                                 \
0376             /* This is a special case, not an optimization          \
0377                (X##_f[3]/Y##_f[3] would not fit into UWtype).           \
0378                As X## is guaranteed to be < Y,  R##_f[_i] can be either     \
0379                (UWtype)-1 or (UWtype)-2.  */                    \
0380             R##_f[_i] = -1;                         \
0381             if (!_i)                                \
0382           break;                                \
0383             __FP_FRAC_SUB_4(X##_f[3], X##_f[2], X##_f[1], X##_f[0],     \
0384                 Y##_f[2], Y##_f[1], Y##_f[0], 0,            \
0385                 X##_f[2], X##_f[1], X##_f[0], _n_f[_i]);        \
0386             _FP_FRAC_SUB_4(X, Y, X);                        \
0387             if (X##_f[3] > Y##_f[3])                        \
0388               {                                 \
0389                 R##_f[_i] = -2;                         \
0390                 _FP_FRAC_ADD_4(X, Y, X);                    \
0391               }                                 \
0392           }                                 \
0393         else                                    \
0394           {                                 \
0395             udiv_qrnnd(R##_f[_i], X##_f[3], X##_f[3], X##_f[2], Y##_f[3]);  \
0396             umul_ppppmnnn(_m_f[3], _m_f[2], _m_f[1], _m_f[0],           \
0397               R##_f[_i], Y##_f[2], Y##_f[1], Y##_f[0]);     \
0398             X##_f[2] = X##_f[1];                        \
0399             X##_f[1] = X##_f[0];                        \
0400             X##_f[0] = _n_f[_i];                        \
0401             if (_FP_FRAC_GT_4(_m, X))                       \
0402               {                                 \
0403                 R##_f[_i]--;                            \
0404                 _FP_FRAC_ADD_4(X, Y, X);                    \
0405                 if (_FP_FRAC_GE_4(X, Y) && _FP_FRAC_GT_4(_m, X))        \
0406                   {                             \
0407             R##_f[_i]--;                        \
0408             _FP_FRAC_ADD_4(X, Y, X);                    \
0409                   }                             \
0410               }                                 \
0411             _FP_FRAC_DEC_4(X, _m);                      \
0412             if (!_i)                                \
0413           {                                 \
0414         if (!_FP_FRAC_EQ_4(X, _m))                  \
0415           R##_f[0] |= _FP_WORK_STICKY;                  \
0416         break;                              \
0417           }                                 \
0418           }                                 \
0419       }                                     \
0420   } while (0)
0421 
0422 
0423 /*
0424  * Square root algorithms:
0425  * We have just one right now, maybe Newton approximation
0426  * should be added for those machines where division is fast.
0427  */
0428  
0429 #define _FP_SQRT_MEAT_4(R, S, T, X, q)              \
0430   do {                              \
0431     while (q)                           \
0432       {                             \
0433     T##_f[3] = S##_f[3] + q;                \
0434     if (T##_f[3] <= X##_f[3])               \
0435       {                         \
0436         S##_f[3] = T##_f[3] + q;                \
0437         X##_f[3] -= T##_f[3];               \
0438         R##_f[3] += q;                  \
0439       }                         \
0440     _FP_FRAC_SLL_4(X, 1);                   \
0441     q >>= 1;                        \
0442       }                             \
0443     q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1);         \
0444     while (q)                           \
0445       {                             \
0446     T##_f[2] = S##_f[2] + q;                \
0447     T##_f[3] = S##_f[3];                    \
0448     if (T##_f[3] < X##_f[3] ||              \
0449         (T##_f[3] == X##_f[3] && T##_f[2] <= X##_f[2])) \
0450       {                         \
0451         S##_f[2] = T##_f[2] + q;                \
0452         S##_f[3] += (T##_f[2] > S##_f[2]);          \
0453         __FP_FRAC_DEC_2(X##_f[3], X##_f[2],         \
0454                 T##_f[3], T##_f[2]);        \
0455         R##_f[2] += q;                  \
0456       }                         \
0457     _FP_FRAC_SLL_4(X, 1);                   \
0458     q >>= 1;                        \
0459       }                             \
0460     q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1);         \
0461     while (q)                           \
0462       {                             \
0463     T##_f[1] = S##_f[1] + q;                \
0464     T##_f[2] = S##_f[2];                    \
0465     T##_f[3] = S##_f[3];                    \
0466     if (T##_f[3] < X##_f[3] ||              \
0467         (T##_f[3] == X##_f[3] && (T##_f[2] < X##_f[2] ||    \
0468          (T##_f[2] == X##_f[2] && T##_f[1] <= X##_f[1]))))  \
0469       {                         \
0470         S##_f[1] = T##_f[1] + q;                \
0471         S##_f[2] += (T##_f[1] > S##_f[1]);          \
0472         S##_f[3] += (T##_f[2] > S##_f[2]);          \
0473         __FP_FRAC_DEC_3(X##_f[3], X##_f[2], X##_f[1],   \
0474                     T##_f[3], T##_f[2], T##_f[1]);  \
0475         R##_f[1] += q;                  \
0476       }                         \
0477     _FP_FRAC_SLL_4(X, 1);                   \
0478     q >>= 1;                        \
0479       }                             \
0480     q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1);         \
0481     while (q != _FP_WORK_ROUND)                 \
0482       {                             \
0483     T##_f[0] = S##_f[0] + q;                \
0484     T##_f[1] = S##_f[1];                    \
0485     T##_f[2] = S##_f[2];                    \
0486     T##_f[3] = S##_f[3];                    \
0487     if (_FP_FRAC_GE_4(X,T))                 \
0488       {                         \
0489         S##_f[0] = T##_f[0] + q;                \
0490         S##_f[1] += (T##_f[0] > S##_f[0]);          \
0491         S##_f[2] += (T##_f[1] > S##_f[1]);          \
0492         S##_f[3] += (T##_f[2] > S##_f[2]);          \
0493         _FP_FRAC_DEC_4(X, T);               \
0494         R##_f[0] += q;                  \
0495       }                         \
0496     _FP_FRAC_SLL_4(X, 1);                   \
0497     q >>= 1;                        \
0498       }                             \
0499     if (!_FP_FRAC_ZEROP_4(X))                   \
0500       {                             \
0501     if (_FP_FRAC_GT_4(X,S))                 \
0502       R##_f[0] |= _FP_WORK_ROUND;               \
0503     R##_f[0] |= _FP_WORK_STICKY;                \
0504       }                             \
0505   } while (0)
0506 
0507 
0508 /*
0509  * Internals 
0510  */
0511 
0512 #define __FP_FRAC_SET_4(X,I3,I2,I1,I0)                  \
0513   (X##_f[3] = I3, X##_f[2] = I2, X##_f[1] = I1, X##_f[0] = I0)
0514 
0515 #ifndef __FP_FRAC_ADD_3
0516 #define __FP_FRAC_ADD_3(r2,r1,r0,x2,x1,x0,y2,y1,y0)     \
0517   do {                              \
0518     int _c1, _c2;                           \
0519     r0 = x0 + y0;                       \
0520     _c1 = r0 < x0;                      \
0521     r1 = x1 + y1;                       \
0522     _c2 = r1 < x1;                      \
0523     r1 += _c1;                          \
0524     _c2 |= r1 < _c1;                        \
0525     r2 = x2 + y2 + _c2;                     \
0526   } while (0)
0527 #endif
0528 
0529 #ifndef __FP_FRAC_ADD_4
0530 #define __FP_FRAC_ADD_4(r3,r2,r1,r0,x3,x2,x1,x0,y3,y2,y1,y0)    \
0531   do {                              \
0532     int _c1, _c2, _c3;                      \
0533     r0 = x0 + y0;                       \
0534     _c1 = r0 < x0;                      \
0535     r1 = x1 + y1;                       \
0536     _c2 = r1 < x1;                      \
0537     r1 += _c1;                          \
0538     _c2 |= r1 < _c1;                        \
0539     r2 = x2 + y2;                       \
0540     _c3 = r2 < x2;                      \
0541     r2 += _c2;                          \
0542     _c3 |= r2 < _c2;                        \
0543     r3 = x3 + y3 + _c3;                     \
0544   } while (0)
0545 #endif
0546 
0547 #ifndef __FP_FRAC_SUB_3
0548 #define __FP_FRAC_SUB_3(r2,r1,r0,x2,x1,x0,y2,y1,y0)     \
0549   do {                              \
0550     int _c1, _c2;                           \
0551     r0 = x0 - y0;                       \
0552     _c1 = r0 > x0;                      \
0553     r1 = x1 - y1;                       \
0554     _c2 = r1 > x1;                      \
0555     r1 -= _c1;                          \
0556     _c2 |= r1 > _c1;                        \
0557     r2 = x2 - y2 - _c2;                     \
0558   } while (0)
0559 #endif
0560 
0561 #ifndef __FP_FRAC_SUB_4
0562 #define __FP_FRAC_SUB_4(r3,r2,r1,r0,x3,x2,x1,x0,y3,y2,y1,y0)    \
0563   do {                              \
0564     int _c1, _c2, _c3;                      \
0565     r0 = x0 - y0;                       \
0566     _c1 = r0 > x0;                      \
0567     r1 = x1 - y1;                       \
0568     _c2 = r1 > x1;                      \
0569     r1 -= _c1;                          \
0570     _c2 |= r1 > _c1;                        \
0571     r2 = x2 - y2;                       \
0572     _c3 = r2 > x2;                      \
0573     r2 -= _c2;                          \
0574     _c3 |= r2 > _c2;                        \
0575     r3 = x3 - y3 - _c3;                     \
0576   } while (0)
0577 #endif
0578 
0579 #ifndef __FP_FRAC_DEC_3
0580 #define __FP_FRAC_DEC_3(x2,x1,x0,y2,y1,y0)              \
0581   do {                                  \
0582     UWtype _t0, _t1, _t2;                       \
0583     _t0 = x0, _t1 = x1, _t2 = x2;                   \
0584     __FP_FRAC_SUB_3 (x2, x1, x0, _t2, _t1, _t0, y2, y1, y0);        \
0585   } while (0)
0586 #endif
0587 
0588 #ifndef __FP_FRAC_DEC_4
0589 #define __FP_FRAC_DEC_4(x3,x2,x1,x0,y3,y2,y1,y0)            \
0590   do {                                  \
0591     UWtype _t0, _t1, _t2, _t3;                      \
0592     _t0 = x0, _t1 = x1, _t2 = x2, _t3 = x3;             \
0593     __FP_FRAC_SUB_4 (x3,x2,x1,x0,_t3,_t2,_t1,_t0, y3,y2,y1,y0);     \
0594   } while (0)
0595 #endif
0596 
0597 #ifndef __FP_FRAC_ADDI_4
0598 #define __FP_FRAC_ADDI_4(x3,x2,x1,x0,i)                 \
0599   do {                                  \
0600     UWtype _t;                              \
0601     _t = ((x0 += i) < i);                       \
0602     x1 += _t; _t = (x1 < _t);                       \
0603     x2 += _t; _t = (x2 < _t);                       \
0604     x3 += _t;                               \
0605   } while (0)
0606 #endif
0607 
0608 /* Convert FP values between word sizes. This appears to be more
0609  * complicated than I'd have expected it to be, so these might be
0610  * wrong... These macros are in any case somewhat bogus because they
0611  * use information about what various FRAC_n variables look like 
0612  * internally [eg, that 2 word vars are X_f0 and x_f1]. But so do
0613  * the ones in op-2.h and op-1.h. 
0614  */
0615 #define _FP_FRAC_CONV_1_4(dfs, sfs, D, S)               \
0616    do {                                 \
0617      if (S##_c != FP_CLS_NAN)                       \
0618        _FP_FRAC_SRS_4(S, (_FP_WFRACBITS_##sfs - _FP_WFRACBITS_##dfs),   \
0619               _FP_WFRACBITS_##sfs);             \
0620      else                               \
0621        _FP_FRAC_SRL_4(S, (_FP_WFRACBITS_##sfs - _FP_WFRACBITS_##dfs));  \
0622      D##_f = S##_f[0];                          \
0623   } while (0)
0624 
0625 #define _FP_FRAC_CONV_2_4(dfs, sfs, D, S)               \
0626    do {                                 \
0627      if (S##_c != FP_CLS_NAN)                       \
0628        _FP_FRAC_SRS_4(S, (_FP_WFRACBITS_##sfs - _FP_WFRACBITS_##dfs),   \
0629               _FP_WFRACBITS_##sfs);             \
0630      else                               \
0631        _FP_FRAC_SRL_4(S, (_FP_WFRACBITS_##sfs - _FP_WFRACBITS_##dfs));  \
0632      D##_f0 = S##_f[0];                         \
0633      D##_f1 = S##_f[1];                         \
0634   } while (0)
0635 
0636 /* Assembly/disassembly for converting to/from integral types.  
0637  * No shifting or overflow handled here.
0638  */
0639 /* Put the FP value X into r, which is an integer of size rsize. */
0640 #define _FP_FRAC_ASSEMBLE_4(r, X, rsize)                \
0641   do {                                  \
0642     if (rsize <= _FP_W_TYPE_SIZE)                   \
0643       r = X##_f[0];                         \
0644     else if (rsize <= 2*_FP_W_TYPE_SIZE)                \
0645     {                                   \
0646       r = X##_f[1];                         \
0647       r <<= _FP_W_TYPE_SIZE;                        \
0648       r += X##_f[0];                            \
0649     }                                   \
0650     else                                \
0651     {                                   \
0652       /* I'm feeling lazy so we deal with int == 3words (implausible)*/ \
0653       /* and int == 4words as a single case.             */ \
0654       r = X##_f[3];                         \
0655       r <<= _FP_W_TYPE_SIZE;                        \
0656       r += X##_f[2];                            \
0657       r <<= _FP_W_TYPE_SIZE;                        \
0658       r += X##_f[1];                            \
0659       r <<= _FP_W_TYPE_SIZE;                        \
0660       r += X##_f[0];                            \
0661     }                                   \
0662   } while (0)
0663 
0664 /* "No disassemble Number Five!" */
0665 /* move an integer of size rsize into X's fractional part. We rely on
0666  * the _f[] array consisting of words of size _FP_W_TYPE_SIZE to avoid
0667  * having to mask the values we store into it.
0668  */
0669 #define _FP_FRAC_DISASSEMBLE_4(X, r, rsize)             \
0670   do {                                  \
0671     X##_f[0] = r;                           \
0672     X##_f[1] = (rsize <= _FP_W_TYPE_SIZE ? 0 : r >> _FP_W_TYPE_SIZE);   \
0673     X##_f[2] = (rsize <= 2*_FP_W_TYPE_SIZE ? 0 : r >> 2*_FP_W_TYPE_SIZE); \
0674     X##_f[3] = (rsize <= 3*_FP_W_TYPE_SIZE ? 0 : r >> 3*_FP_W_TYPE_SIZE); \
0675   } while (0)
0676 
0677 #define _FP_FRAC_CONV_4_1(dfs, sfs, D, S)               \
0678    do {                                 \
0679      D##_f[0] = S##_f;                          \
0680      D##_f[1] = D##_f[2] = D##_f[3] = 0;                \
0681      _FP_FRAC_SLL_4(D, (_FP_WFRACBITS_##dfs - _FP_WFRACBITS_##sfs));    \
0682    } while (0)
0683 
0684 #define _FP_FRAC_CONV_4_2(dfs, sfs, D, S)               \
0685    do {                                 \
0686      D##_f[0] = S##_f0;                         \
0687      D##_f[1] = S##_f1;                         \
0688      D##_f[2] = D##_f[3] = 0;                       \
0689      _FP_FRAC_SLL_4(D, (_FP_WFRACBITS_##dfs - _FP_WFRACBITS_##sfs));    \
0690    } while (0)
0691 
0692 #endif