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

   1 /* Software floating-point emulation.
   2    Basic four-word fraction declaration and manipulation.
   3    Copyright (C) 1997,1998,1999,2006,2007 Free Software Foundation, Inc.
   4    This file is part of the GNU C Library.
   5    Contributed by Richard Henderson (rth@cygnus.com),
   6                   Jakub Jelinek (jj@ultra.linux.cz),
   7                   David S. Miller (davem@redhat.com) and
   8                   Peter Maydell (pmaydell@chiark.greenend.org.uk).
   9
  10    The GNU C Library is free software; you can redistribute it and/or
  11    modify it under the terms of the GNU Lesser General Public
  12    License as published by the Free Software Foundation; either
  13    version 2.1 of the License, or (at your option) any later version.
  14
  15    In addition to the permissions in the GNU Lesser General Public
  16    License, the Free Software Foundation gives you unlimited
  17    permission to link the compiled version of this file into
  18    combinations with other programs, and to distribute those
  19    combinations without any restriction coming from the use of this
  20    file.  (The Lesser General Public License restrictions do apply in
  21    other respects; for example, they cover modification of the file,
  22    and distribution when not linked into a combine executable.)
  23
  24    The GNU C Library is distributed in the hope that it will be useful,
  25    but WITHOUT ANY WARRANTY; without even the implied warranty of
  26    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  27    Lesser General Public License for more details.
  28
  29    You should have received a copy of the GNU Lesser General Public
  30    License along with the GNU C Library; if not, write to the Free
  31    Software Foundation, 51 Franklin Street, Fifth Floor, Boston,
  32    MA 02110-1301, USA.  */
  33
  34 #define _FP_FRAC_DECL_4(X)      _FP_W_TYPE X##_f[4]
  35 #define _FP_FRAC_COPY_4(D,S)                    \
  36   (D##_f[0] = S##_f[0], D##_f[1] = S##_f[1],    \
  37    D##_f[2] = S##_f[2], D##_f[3] = S##_f[3])
  38 #define _FP_FRAC_SET_4(X,I)     __FP_FRAC_SET_4(X, I)
  39 #define _FP_FRAC_HIGH_4(X)      (X##_f[3])
  40 #define _FP_FRAC_LOW_4(X)       (X##_f[0])
  41 #define _FP_FRAC_WORD_4(X,w)    (X##_f[w])
  42
  43 #define _FP_FRAC_SLL_4(X,N)                                             \
  44   do {                                                                  \
  45     _FP_I_TYPE _up, _down, _skip, _i;                                   \
  46     _skip = (N) / _FP_W_TYPE_SIZE;                                      \
  47     _up = (N) % _FP_W_TYPE_SIZE;                                        \
  48     _down = _FP_W_TYPE_SIZE - _up;                                      \
  49     if (!_up)                                                           \
  50       for (_i = 3; _i >= _skip; --_i)                                   \
  51         X##_f[_i] = X##_f[_i-_skip];                                    \
  52     else                                                                \
  53       {                                                                 \
  54         for (_i = 3; _i > _skip; --_i)                                  \
  55           X##_f[_i] = X##_f[_i-_skip] << _up                            \
  56                       | X##_f[_i-_skip-1] >> _down;                     \
  57         X##_f[_i--] = X##_f[0] << _up;                                  \
  58       }                                                                 \
  59     for (; _i >= 0; --_i)                                               \
  60       X##_f[_i] = 0;                                                    \
  61   } while (0)
  62
  63 /* This one was broken too */
  64 #define _FP_FRAC_SRL_4(X,N)                                             \
  65   do {                                                                  \
  66     _FP_I_TYPE _up, _down, _skip, _i;                                   \
  67     _skip = (N) / _FP_W_TYPE_SIZE;                                      \
  68     _down = (N) % _FP_W_TYPE_SIZE;                                      \
  69     _up = _FP_W_TYPE_SIZE - _down;                                      \
  70     if (!_down)                                                         \
  71       for (_i = 0; _i <= 3-_skip; ++_i)                                 \
  72         X##_f[_i] = X##_f[_i+_skip];                                    \
  73     else                                                                \
  74       {                                                                 \
  75         for (_i = 0; _i < 3-_skip; ++_i)                                \
  76           X##_f[_i] = X##_f[_i+_skip] >> _down                          \
  77                       | X##_f[_i+_skip+1] << _up;                       \
  78         X##_f[_i++] = X##_f[3] >> _down;                                \
  79       }                                                                 \
  80     for (; _i < 4; ++_i)                                                \
  81       X##_f[_i] = 0;                                                    \
  82   } while (0)
  83
  84
  85 /* Right shift with sticky-lsb.
  86  * What this actually means is that we do a standard right-shift,
  87  * but that if any of the bits that fall off the right hand side
  88  * were one then we always set the LSbit.
  89  */
  90 #define _FP_FRAC_SRST_4(X,S,N,size)                     \
  91   do {                                                  \
  92     _FP_I_TYPE _up, _down, _skip, _i;                   \
  93     _FP_W_TYPE _s;                                      \
  94     _skip = (N) / _FP_W_TYPE_SIZE;                      \
  95     _down = (N) % _FP_W_TYPE_SIZE;                      \
  96     _up = _FP_W_TYPE_SIZE - _down;                      \
  97     for (_s = _i = 0; _i < _skip; ++_i)                 \
  98       _s |= X##_f[_i];                                  \
  99     if (!_down)                                         \
 100       for (_i = 0; _i <= 3-_skip; ++_i)                 \
 101         X##_f[_i] = X##_f[_i+_skip];                    \
 102     else                                                \
 103       {                                                 \
 104         _s |= X##_f[_i] << _up;                         \
 105         for (_i = 0; _i < 3-_skip; ++_i)                \
 106           X##_f[_i] = X##_f[_i+_skip] >> _down          \
 107                       | X##_f[_i+_skip+1] << _up;       \
 108         X##_f[_i++] = X##_f[3] >> _down;                \
 109       }                                                 \
 110     for (; _i < 4; ++_i)                                \
 111       X##_f[_i] = 0;                                    \
 112     S = (_s != 0);                                      \
 113   } while (0)
 114
 115 #define _FP_FRAC_SRS_4(X,N,size)                \
 116   do {                                          \
 117     int _sticky;                                \
 118     _FP_FRAC_SRST_4(X, _sticky, N, size);       \
 119     X##_f[0] |= _sticky;                        \
 120   } while (0)
 121
 122 #define _FP_FRAC_ADD_4(R,X,Y)                                           \
 123   __FP_FRAC_ADD_4(R##_f[3], R##_f[2], R##_f[1], R##_f[0],               \
 124                   X##_f[3], X##_f[2], X##_f[1], X##_f[0],               \
 125                   Y##_f[3], Y##_f[2], Y##_f[1], Y##_f[0])
 126
 127 #define _FP_FRAC_SUB_4(R,X,Y)                                           \
 128   __FP_FRAC_SUB_4(R##_f[3], R##_f[2], R##_f[1], R##_f[0],               \
 129                   X##_f[3], X##_f[2], X##_f[1], X##_f[0],               \
 130                   Y##_f[3], Y##_f[2], Y##_f[1], Y##_f[0])
 131
 132 #define _FP_FRAC_DEC_4(X,Y)                                             \
 133   __FP_FRAC_DEC_4(X##_f[3], X##_f[2], X##_f[1], X##_f[0],               \
 134                   Y##_f[3], Y##_f[2], Y##_f[1], Y##_f[0])
 135
 136 #define _FP_FRAC_ADDI_4(X,I)                                            \
 137   __FP_FRAC_ADDI_4(X##_f[3], X##_f[2], X##_f[1], X##_f[0], I)
 138
 139 #define _FP_ZEROFRAC_4  0,0,0,0
 140 #define _FP_MINFRAC_4   0,0,0,1
 141 #define _FP_MAXFRAC_4   (~(_FP_WS_TYPE)0), (~(_FP_WS_TYPE)0), (~(_FP_WS_TYPE)0), (~(_FP_WS_TYPE)0)
 142
 143 #define _FP_FRAC_ZEROP_4(X)     ((X##_f[0] | X##_f[1] | X##_f[2] | X##_f[3]) == 0)
 144 #define _FP_FRAC_NEGP_4(X)      ((_FP_WS_TYPE)X##_f[3] < 0)
 145 #define _FP_FRAC_OVERP_4(fs,X)  (_FP_FRAC_HIGH_##fs(X) & _FP_OVERFLOW_##fs)
 146 #define _FP_FRAC_CLEAR_OVERP_4(fs,X)  (_FP_FRAC_HIGH_##fs(X) &= ~_FP_OVERFLOW_##fs)
 147
 148 #define _FP_FRAC_EQ_4(X,Y)                              \
 149  (X##_f[0] == Y##_f[0] && X##_f[1] == Y##_f[1]          \
 150   && X##_f[2] == Y##_f[2] && X##_f[3] == Y##_f[3])
 151
 152 #define _FP_FRAC_GT_4(X,Y)                              \
 153  (X##_f[3] > Y##_f[3] ||                                \
 154   (X##_f[3] == Y##_f[3] && (X##_f[2] > Y##_f[2] ||      \
 155    (X##_f[2] == Y##_f[2] && (X##_f[1] > Y##_f[1] ||     \
 156     (X##_f[1] == Y##_f[1] && X##_f[0] > Y##_f[0])       \
 157    ))                                                   \
 158   ))                                                    \
 159  )
 160
 161 #define _FP_FRAC_GE_4(X,Y)                              \
 162  (X##_f[3] > Y##_f[3] ||                                \
 163   (X##_f[3] == Y##_f[3] && (X##_f[2] > Y##_f[2] ||      \
 164    (X##_f[2] == Y##_f[2] && (X##_f[1] > Y##_f[1] ||     \
 165     (X##_f[1] == Y##_f[1] && X##_f[0] >= Y##_f[0])      \
 166    ))                                                   \
 167   ))                                                    \
 168  )
 169
 170
 171 #define _FP_FRAC_CLZ_4(R,X)             \
 172   do {                                  \
 173     if (X##_f[3])                       \
 174     {                                   \
 175         __FP_CLZ(R,X##_f[3]);           \
 176     }                                   \
 177     else if (X##_f[2])                  \
 178     {                                   \
 179         __FP_CLZ(R,X##_f[2]);           \
 180         R += _FP_W_TYPE_SIZE;           \
 181     }                                   \
 182     else if (X##_f[1])                  \
 183     {                                   \
 184         __FP_CLZ(R,X##_f[1]);           \
 185         R += _FP_W_TYPE_SIZE*2;         \
 186     }                                   \
 187     else                                \
 188     {                                   \
 189         __FP_CLZ(R,X##_f[0]);           \
 190         R += _FP_W_TYPE_SIZE*3;         \
 191     }                                   \
 192   } while(0)
 193
 194
 195 #define _FP_UNPACK_RAW_4(fs, X, val)                            \
 196   do {                                                          \
 197     union _FP_UNION_##fs _flo; _flo.flt = (val);                \
 198     X##_f[0] = _flo.bits.frac0;                                 \
 199     X##_f[1] = _flo.bits.frac1;                                 \
 200     X##_f[2] = _flo.bits.frac2;                                 \
 201     X##_f[3] = _flo.bits.frac3;                                 \
 202     X##_e  = _flo.bits.exp;                                     \
 203     X##_s  = _flo.bits.sign;                                    \
 204   } while (0)
 205
 206 #define _FP_UNPACK_RAW_4_P(fs, X, val)                          \
 207   do {                                                          \
 208     union _FP_UNION_##fs *_flo =                                \
 209       (union _FP_UNION_##fs *)(val);                            \
 210                                                                 \
 211     X##_f[0] = _flo->bits.frac0;                                \
 212     X##_f[1] = _flo->bits.frac1;                                \
 213     X##_f[2] = _flo->bits.frac2;                                \
 214     X##_f[3] = _flo->bits.frac3;                                \
 215     X##_e  = _flo->bits.exp;                                    \
 216     X##_s  = _flo->bits.sign;                                   \
 217   } while (0)
 218
 219 #define _FP_PACK_RAW_4(fs, val, X)                              \
 220   do {                                                          \
 221     union _FP_UNION_##fs _flo;                                  \
 222     _flo.bits.frac0 = X##_f[0];                                 \
 223     _flo.bits.frac1 = X##_f[1];                                 \
 224     _flo.bits.frac2 = X##_f[2];                                 \
 225     _flo.bits.frac3 = X##_f[3];                                 \
 226     _flo.bits.exp   = X##_e;                                    \
 227     _flo.bits.sign  = X##_s;                                    \
 228     (val) = _flo.flt;                                           \
 229   } while (0)
 230
 231 #define _FP_PACK_RAW_4_P(fs, val, X)                            \
 232   do {                                                          \
 233     union _FP_UNION_##fs *_flo =                                \
 234       (union _FP_UNION_##fs *)(val);                            \
 235                                                                 \
 236     _flo->bits.frac0 = X##_f[0];                                \
 237     _flo->bits.frac1 = X##_f[1];                                \
 238     _flo->bits.frac2 = X##_f[2];                                \
 239     _flo->bits.frac3 = X##_f[3];                                \
 240     _flo->bits.exp   = X##_e;                                   \
 241     _flo->bits.sign  = X##_s;                                   \
 242   } while (0)
 243
 244 /*
 245  * Multiplication algorithms:
 246  */
 247
 248 /* Given a 1W * 1W => 2W primitive, do the extended multiplication.  */
 249
 250 #define _FP_MUL_MEAT_4_wide(wfracbits, R, X, Y, doit)                       \
 251   do {                                                                      \
 252     _FP_FRAC_DECL_8(_z); _FP_FRAC_DECL_2(_b); _FP_FRAC_DECL_2(_c);          \
 253     _FP_FRAC_DECL_2(_d); _FP_FRAC_DECL_2(_e); _FP_FRAC_DECL_2(_f);          \
 254                                                                             \
 255     doit(_FP_FRAC_WORD_8(_z,1), _FP_FRAC_WORD_8(_z,0), X##_f[0], Y##_f[0]); \
 256     doit(_b_f1, _b_f0, X##_f[0], Y##_f[1]);                                 \
 257     doit(_c_f1, _c_f0, X##_f[1], Y##_f[0]);                                 \
 258     doit(_d_f1, _d_f0, X##_f[1], Y##_f[1]);                                 \
 259     doit(_e_f1, _e_f0, X##_f[0], Y##_f[2]);                                 \
 260     doit(_f_f1, _f_f0, X##_f[2], Y##_f[0]);                                 \
 261     __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,3),_FP_FRAC_WORD_8(_z,2),            \
 262                     _FP_FRAC_WORD_8(_z,1), 0,_b_f1,_b_f0,                   \
 263                     0,0,_FP_FRAC_WORD_8(_z,1));                             \
 264     __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,3),_FP_FRAC_WORD_8(_z,2),            \
 265                     _FP_FRAC_WORD_8(_z,1), 0,_c_f1,_c_f0,                   \
 266                     _FP_FRAC_WORD_8(_z,3),_FP_FRAC_WORD_8(_z,2),            \
 267                     _FP_FRAC_WORD_8(_z,1));                                 \
 268     __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,4),_FP_FRAC_WORD_8(_z,3),            \
 269                     _FP_FRAC_WORD_8(_z,2), 0,_d_f1,_d_f0,                   \
 270                     0,_FP_FRAC_WORD_8(_z,3),_FP_FRAC_WORD_8(_z,2));         \
 271     __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,4),_FP_FRAC_WORD_8(_z,3),            \
 272                     _FP_FRAC_WORD_8(_z,2), 0,_e_f1,_e_f0,                   \
 273                     _FP_FRAC_WORD_8(_z,4),_FP_FRAC_WORD_8(_z,3),            \
 274                     _FP_FRAC_WORD_8(_z,2));                                 \
 275     __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,4),_FP_FRAC_WORD_8(_z,3),            \
 276                     _FP_FRAC_WORD_8(_z,2), 0,_f_f1,_f_f0,                   \
 277                     _FP_FRAC_WORD_8(_z,4),_FP_FRAC_WORD_8(_z,3),            \
 278                     _FP_FRAC_WORD_8(_z,2));                                 \
 279     doit(_b_f1, _b_f0, X##_f[0], Y##_f[3]);                                 \
 280     doit(_c_f1, _c_f0, X##_f[3], Y##_f[0]);                                 \
 281     doit(_d_f1, _d_f0, X##_f[1], Y##_f[2]);                                 \
 282     doit(_e_f1, _e_f0, X##_f[2], Y##_f[1]);                                 \
 283     __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,5),_FP_FRAC_WORD_8(_z,4),            \
 284                     _FP_FRAC_WORD_8(_z,3), 0,_b_f1,_b_f0,                   \
 285                     0,_FP_FRAC_WORD_8(_z,4),_FP_FRAC_WORD_8(_z,3));         \
 286     __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,5),_FP_FRAC_WORD_8(_z,4),            \
 287                     _FP_FRAC_WORD_8(_z,3), 0,_c_f1,_c_f0,                   \
 288                     _FP_FRAC_WORD_8(_z,5),_FP_FRAC_WORD_8(_z,4),            \
 289                     _FP_FRAC_WORD_8(_z,3));                                 \
 290     __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,5),_FP_FRAC_WORD_8(_z,4),            \
 291                     _FP_FRAC_WORD_8(_z,3), 0,_d_f1,_d_f0,                   \
 292                     _FP_FRAC_WORD_8(_z,5),_FP_FRAC_WORD_8(_z,4),            \
 293                     _FP_FRAC_WORD_8(_z,3));                                 \
 294     __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,5),_FP_FRAC_WORD_8(_z,4),            \
 295                     _FP_FRAC_WORD_8(_z,3), 0,_e_f1,_e_f0,                   \
 296                     _FP_FRAC_WORD_8(_z,5),_FP_FRAC_WORD_8(_z,4),            \
 297                     _FP_FRAC_WORD_8(_z,3));                                 \
 298     doit(_b_f1, _b_f0, X##_f[2], Y##_f[2]);                                 \
 299     doit(_c_f1, _c_f0, X##_f[1], Y##_f[3]);                                 \
 300     doit(_d_f1, _d_f0, X##_f[3], Y##_f[1]);                                 \
 301     doit(_e_f1, _e_f0, X##_f[2], Y##_f[3]);                                 \
 302     doit(_f_f1, _f_f0, X##_f[3], Y##_f[2]);                                 \
 303     __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,6),_FP_FRAC_WORD_8(_z,5),            \
 304                     _FP_FRAC_WORD_8(_z,4), 0,_b_f1,_b_f0,                   \
 305                     0,_FP_FRAC_WORD_8(_z,5),_FP_FRAC_WORD_8(_z,4));         \
 306     __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,6),_FP_FRAC_WORD_8(_z,5),            \
 307                     _FP_FRAC_WORD_8(_z,4), 0,_c_f1,_c_f0,                   \
 308                     _FP_FRAC_WORD_8(_z,6),_FP_FRAC_WORD_8(_z,5),            \
 309                     _FP_FRAC_WORD_8(_z,4));                                 \
 310     __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,6),_FP_FRAC_WORD_8(_z,5),            \
 311                     _FP_FRAC_WORD_8(_z,4), 0,_d_f1,_d_f0,                   \
 312                     _FP_FRAC_WORD_8(_z,6),_FP_FRAC_WORD_8(_z,5),            \
 313                     _FP_FRAC_WORD_8(_z,4));                                 \
 314     __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,7),_FP_FRAC_WORD_8(_z,6),            \
 315                     _FP_FRAC_WORD_8(_z,5), 0,_e_f1,_e_f0,                   \
 316                     0,_FP_FRAC_WORD_8(_z,6),_FP_FRAC_WORD_8(_z,5));         \
 317     __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,7),_FP_FRAC_WORD_8(_z,6),            \
 318                     _FP_FRAC_WORD_8(_z,5), 0,_f_f1,_f_f0,                   \
 319                     _FP_FRAC_WORD_8(_z,7),_FP_FRAC_WORD_8(_z,6),            \
 320                     _FP_FRAC_WORD_8(_z,5));                                 \
 321     doit(_b_f1, _b_f0, X##_f[3], Y##_f[3]);                                 \
 322     __FP_FRAC_ADD_2(_FP_FRAC_WORD_8(_z,7),_FP_FRAC_WORD_8(_z,6),            \
 323                     _b_f1,_b_f0,                                            \
 324                     _FP_FRAC_WORD_8(_z,7),_FP_FRAC_WORD_8(_z,6));           \
 325                                                                             \
 326     /* Normalize since we know where the msb of the multiplicands           \
 327        were (bit B), we know that the msb of the of the product is          \
 328        at either 2B or 2B-1.  */                                            \
 329     _FP_FRAC_SRS_8(_z, wfracbits-1, 2*wfracbits);                           \
 330     __FP_FRAC_SET_4(R, _FP_FRAC_WORD_8(_z,3), _FP_FRAC_WORD_8(_z,2),        \
 331                     _FP_FRAC_WORD_8(_z,1), _FP_FRAC_WORD_8(_z,0));          \
 332   } while (0)
 333
 334 #define _FP_MUL_MEAT_4_gmp(wfracbits, R, X, Y)                              \
 335   do {                                                                      \
 336     _FP_FRAC_DECL_8(_z);                                                    \
 337                                                                             \
 338     mpn_mul_n(_z_f, _x_f, _y_f, 4);                                         \
 339                                                                             \
 340     /* Normalize since we know where the msb of the multiplicands           \
 341        were (bit B), we know that the msb of the of the product is          \
 342        at either 2B or 2B-1.  */                                            \
 343     _FP_FRAC_SRS_8(_z, wfracbits-1, 2*wfracbits);                           \
 344     __FP_FRAC_SET_4(R, _FP_FRAC_WORD_8(_z,3), _FP_FRAC_WORD_8(_z,2),        \
 345                     _FP_FRAC_WORD_8(_z,1), _FP_FRAC_WORD_8(_z,0));          \
 346   } while (0)
 347
 348 /*
 349  * Helper utility for _FP_DIV_MEAT_4_udiv:
 350  * pppp = m * nnn
 351  */
 352 #define umul_ppppmnnn(p3,p2,p1,p0,m,n2,n1,n0)                               \
 353   do {                                                                      \
 354     UWtype _t;                                                              \
 355     umul_ppmm(p1,p0,m,n0);                                                  \
 356     umul_ppmm(p2,_t,m,n1);                                                  \
 357     __FP_FRAC_ADDI_2(p2,p1,_t);                                             \
 358     umul_ppmm(p3,_t,m,n2);                                                  \
 359     __FP_FRAC_ADDI_2(p3,p2,_t);                                             \
 360   } while (0)
 361
 362 /*
 363  * Division algorithms:
 364  */
 365
 366 #define _FP_DIV_MEAT_4_udiv(fs, R, X, Y)                                    \
 367   do {                                                                      \
 368     int _i;                                                                 \
 369     _FP_FRAC_DECL_4(_n); _FP_FRAC_DECL_4(_m);                               \
 370     _FP_FRAC_SET_4(_n, _FP_ZEROFRAC_4);                                     \
 371     if (_FP_FRAC_GT_4(X, Y))                                                \
 372       {                                                                     \
 373         _n_f[3] = X##_f[0] << (_FP_W_TYPE_SIZE - 1);                        \
 374         _FP_FRAC_SRL_4(X, 1);                                               \
 375       }                                                                     \
 376     else                                                                    \
 377       R##_e--;                                                              \
 378                                                                             \
 379     /* Normalize, i.e. make the most significant bit of the                 \
 380        denominator set. */                                                  \
 381     _FP_FRAC_SLL_4(Y, _FP_WFRACXBITS_##fs);                                 \
 382                                                                             \
 383     for (_i = 3; ; _i--)                                                    \
 384       {                                                                     \
 385         if (X##_f[3] == Y##_f[3])                                           \
 386           {                                                                 \
 387             /* This is a special case, not an optimization                  \
 388                (X##_f[3]/Y##_f[3] would not fit into UWtype).               \
 389                As X## is guaranteed to be < Y,  R##_f[_i] can be either     \
 390                (UWtype)-1 or (UWtype)-2.  */                                \
 391             R##_f[_i] = -1;                                                 \
 392             if (!_i)                                                        \
 393               break;                                                        \
 394             __FP_FRAC_SUB_4(X##_f[3], X##_f[2], X##_f[1], X##_f[0],         \
 395                             Y##_f[2], Y##_f[1], Y##_f[0], 0,                \
 396                             X##_f[2], X##_f[1], X##_f[0], _n_f[_i]);        \
 397             _FP_FRAC_SUB_4(X, Y, X);                                        \
 398             if (X##_f[3] > Y##_f[3])                                        \
 399               {                                                             \
 400                 R##_f[_i] = -2;                                             \
 401                 _FP_FRAC_ADD_4(X, Y, X);                                    \
 402               }                                                             \
 403           }                                                                 \
 404         else                                                                \
 405           {                                                                 \
 406             udiv_qrnnd(R##_f[_i], X##_f[3], X##_f[3], X##_f[2], Y##_f[3]);  \
 407             umul_ppppmnnn(_m_f[3], _m_f[2], _m_f[1], _m_f[0],               \
 408                           R##_f[_i], Y##_f[2], Y##_f[1], Y##_f[0]);         \
 409             X##_f[2] = X##_f[1];                                            \
 410             X##_f[1] = X##_f[0];                                            \
 411             X##_f[0] = _n_f[_i];                                            \
 412             if (_FP_FRAC_GT_4(_m, X))                                       \
 413               {                                                             \
 414                 R##_f[_i]--;                                                \
 415                 _FP_FRAC_ADD_4(X, Y, X);                                    \
 416                 if (_FP_FRAC_GE_4(X, Y) && _FP_FRAC_GT_4(_m, X))            \
 417                   {                                                         \
 418                     R##_f[_i]--;                                            \
 419                     _FP_FRAC_ADD_4(X, Y, X);                                \
 420                   }                                                         \
 421               }                                                             \
 422             _FP_FRAC_DEC_4(X, _m);                                          \
 423             if (!_i)                                                        \
 424               {                                                             \
 425                 if (!_FP_FRAC_EQ_4(X, _m))                                  \
 426                   R##_f[0] |= _FP_WORK_STICKY;                              \
 427                 break;                                                      \
 428               }                                                             \
 429           }                                                                 \
 430       }                                                                     \
 431   } while (0)
 432
 433
 434 /*
 435  * Square root algorithms:
 436  * We have just one right now, maybe Newton approximation
 437  * should be added for those machines where division is fast.
 438  */
 439
 440 #define _FP_SQRT_MEAT_4(R, S, T, X, q)                          \
 441   do {                                                          \
 442     while (q)                                                   \
 443       {                                                         \
 444         T##_f[3] = S##_f[3] + q;                                \
 445         if (T##_f[3] <= X##_f[3])                               \
 446           {                                                     \
 447             S##_f[3] = T##_f[3] + q;                            \
 448             X##_f[3] -= T##_f[3];                               \
 449             R##_f[3] += q;                                      \
 450           }                                                     \
 451         _FP_FRAC_SLL_4(X, 1);                                   \
 452         q >>= 1;                                                \
 453       }                                                         \
 454     q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1);                 \
 455     while (q)                                                   \
 456       {                                                         \
 457         T##_f[2] = S##_f[2] + q;                                \
 458         T##_f[3] = S##_f[3];                                    \
 459         if (T##_f[3] < X##_f[3] ||                              \
 460             (T##_f[3] == X##_f[3] && T##_f[2] <= X##_f[2]))     \
 461           {                                                     \
 462             S##_f[2] = T##_f[2] + q;                            \
 463             S##_f[3] += (T##_f[2] > S##_f[2]);                  \
 464             __FP_FRAC_DEC_2(X##_f[3], X##_f[2],                 \
 465                             T##_f[3], T##_f[2]);                \
 466             R##_f[2] += q;                                      \
 467           }                                                     \
 468         _FP_FRAC_SLL_4(X, 1);                                   \
 469         q >>= 1;                                                \
 470       }                                                         \
 471     q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1);                 \
 472     while (q)                                                   \
 473       {                                                         \
 474         T##_f[1] = S##_f[1] + q;                                \
 475         T##_f[2] = S##_f[2];                                    \
 476         T##_f[3] = S##_f[3];                                    \
 477         if (T##_f[3] < X##_f[3] ||                              \
 478             (T##_f[3] == X##_f[3] && (T##_f[2] < X##_f[2] ||    \
 479              (T##_f[2] == X##_f[2] && T##_f[1] <= X##_f[1]))))  \
 480           {                                                     \
 481             S##_f[1] = T##_f[1] + q;                            \
 482             S##_f[2] += (T##_f[1] > S##_f[1]);                  \
 483             S##_f[3] += (T##_f[2] > S##_f[2]);                  \
 484             __FP_FRAC_DEC_3(X##_f[3], X##_f[2], X##_f[1],       \
 485                             T##_f[3], T##_f[2], T##_f[1]);      \
 486             R##_f[1] += q;                                      \
 487           }                                                     \
 488         _FP_FRAC_SLL_4(X, 1);                                   \
 489         q >>= 1;                                                \
 490       }                                                         \
 491     q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1);                 \
 492     while (q != _FP_WORK_ROUND)                                 \
 493       {                                                         \
 494         T##_f[0] = S##_f[0] + q;                                \
 495         T##_f[1] = S##_f[1];                                    \
 496         T##_f[2] = S##_f[2];                                    \
 497         T##_f[3] = S##_f[3];                                    \
 498         if (_FP_FRAC_GE_4(X,T))                                 \
 499           {                                                     \
 500             S##_f[0] = T##_f[0] + q;                            \
 501             S##_f[1] += (T##_f[0] > S##_f[0]);                  \
 502             S##_f[2] += (T##_f[1] > S##_f[1]);                  \
 503             S##_f[3] += (T##_f[2] > S##_f[2]);                  \
 504             _FP_FRAC_DEC_4(X, T);                               \
 505             R##_f[0] += q;                                      \
 506           }                                                     \
 507         _FP_FRAC_SLL_4(X, 1);                                   \
 508         q >>= 1;                                                \
 509       }                                                         \
 510     if (!_FP_FRAC_ZEROP_4(X))                                   \
 511       {                                                         \
 512         if (_FP_FRAC_GT_4(X,S))                                 \
 513           R##_f[0] |= _FP_WORK_ROUND;                           \
 514         R##_f[0] |= _FP_WORK_STICKY;                            \
 515       }                                                         \
 516   } while (0)
 517
 518
 519 /*
 520  * Internals
 521  */
 522
 523 #define __FP_FRAC_SET_4(X,I3,I2,I1,I0)                                  \
 524   (X##_f[3] = I3, X##_f[2] = I2, X##_f[1] = I1, X##_f[0] = I0)
 525
 526 #ifndef __FP_FRAC_ADD_3
 527 #define __FP_FRAC_ADD_3(r2,r1,r0,x2,x1,x0,y2,y1,y0)             \
 528   do {                                                          \
 529     _FP_W_TYPE _c1, _c2;                                        \
 530     r0 = x0 + y0;                                               \
 531     _c1 = r0 < x0;                                              \
 532     r1 = x1 + y1;                                               \
 533     _c2 = r1 < x1;                                              \
 534     r1 += _c1;                                                  \
 535     _c2 |= r1 < _c1;                                            \
 536     r2 = x2 + y2 + _c2;                                         \
 537   } while (0)
 538 #endif
 539
 540 #ifndef __FP_FRAC_ADD_4
 541 #define __FP_FRAC_ADD_4(r3,r2,r1,r0,x3,x2,x1,x0,y3,y2,y1,y0)    \
 542   do {                                                          \
 543     _FP_W_TYPE _c1, _c2, _c3;                                   \
 544     r0 = x0 + y0;                                               \
 545     _c1 = r0 < x0;                                              \
 546     r1 = x1 + y1;                                               \
 547     _c2 = r1 < x1;                                              \
 548     r1 += _c1;                                                  \
 549     _c2 |= r1 < _c1;                                            \
 550     r2 = x2 + y2;                                               \
 551     _c3 = r2 < x2;                                              \
 552     r2 += _c2;                                                  \
 553     _c3 |= r2 < _c2;                                            \
 554     r3 = x3 + y3 + _c3;                                         \
 555   } while (0)
 556 #endif
 557
 558 #ifndef __FP_FRAC_SUB_3
 559 #define __FP_FRAC_SUB_3(r2,r1,r0,x2,x1,x0,y2,y1,y0)             \
 560   do {                                                          \
 561     _FP_W_TYPE _c1, _c2;                                        \
 562     r0 = x0 - y0;                                               \
 563     _c1 = r0 > x0;                                              \
 564     r1 = x1 - y1;                                               \
 565     _c2 = r1 > x1;                                              \
 566     r1 -= _c1;                                                  \
 567     _c2 |= _c1 && (y1 == x1);                                   \
 568     r2 = x2 - y2 - _c2;                                         \
 569   } while (0)
 570 #endif
 571
 572 #ifndef __FP_FRAC_SUB_4
 573 #define __FP_FRAC_SUB_4(r3,r2,r1,r0,x3,x2,x1,x0,y3,y2,y1,y0)    \
 574   do {                                                          \
 575     _FP_W_TYPE _c1, _c2, _c3;                                   \
 576     r0 = x0 - y0;                                               \
 577     _c1 = r0 > x0;                                              \
 578     r1 = x1 - y1;                                               \
 579     _c2 = r1 > x1;                                              \
 580     r1 -= _c1;                                                  \
 581     _c2 |= _c1 && (y1 == x1);                                   \
 582     r2 = x2 - y2;                                               \
 583     _c3 = r2 > x2;                                              \
 584     r2 -= _c2;                                                  \
 585     _c3 |= _c2 && (y2 == x2);                                   \
 586     r3 = x3 - y3 - _c3;                                         \
 587   } while (0)
 588 #endif
 589
 590 #ifndef __FP_FRAC_DEC_3
 591 #define __FP_FRAC_DEC_3(x2,x1,x0,y2,y1,y0)                              \
 592   do {                                                                  \
 593     UWtype _t0, _t1, _t2;                                               \
 594     _t0 = x0, _t1 = x1, _t2 = x2;                                       \
 595     __FP_FRAC_SUB_3 (x2, x1, x0, _t2, _t1, _t0, y2, y1, y0);            \
 596   } while (0)
 597 #endif
 598
 599 #ifndef __FP_FRAC_DEC_4
 600 #define __FP_FRAC_DEC_4(x3,x2,x1,x0,y3,y2,y1,y0)                        \
 601   do {                                                                  \
 602     UWtype _t0, _t1, _t2, _t3;                                          \
 603     _t0 = x0, _t1 = x1, _t2 = x2, _t3 = x3;                             \
 604     __FP_FRAC_SUB_4 (x3,x2,x1,x0,_t3,_t2,_t1,_t0, y3,y2,y1,y0);         \
 605   } while (0)
 606 #endif
 607
 608 #ifndef __FP_FRAC_ADDI_4
 609 #define __FP_FRAC_ADDI_4(x3,x2,x1,x0,i)                                 \
 610   do {                                                                  \
 611     UWtype _t;                                                          \
 612     _t = ((x0 += i) < i);                                               \
 613     x1 += _t; _t = (x1 < _t);                                           \
 614     x2 += _t; _t = (x2 < _t);                                           \
 615     x3 += _t;                                                           \
 616   } while (0)
 617 #endif
 618
 619 /* Convert FP values between word sizes. This appears to be more
 620  * complicated than I'd have expected it to be, so these might be
 621  * wrong... These macros are in any case somewhat bogus because they
 622  * use information about what various FRAC_n variables look like
 623  * internally [eg, that 2 word vars are X_f0 and x_f1]. But so do
 624  * the ones in op-2.h and op-1.h.
 625  */
 626 #define _FP_FRAC_COPY_1_4(D, S)         (D##_f = S##_f[0])
 627
 628 #define _FP_FRAC_COPY_2_4(D, S)                 \
 629 do {                                            \
 630   D##_f0 = S##_f[0];                            \
 631   D##_f1 = S##_f[1];                            \
 632 } while (0)
 633
 634 /* Assembly/disassembly for converting to/from integral types.
 635  * No shifting or overflow handled here.
 636  */
 637 /* Put the FP value X into r, which is an integer of size rsize. */
 638 #define _FP_FRAC_ASSEMBLE_4(r, X, rsize)                                \
 639   do {                                                                  \
 640     if (rsize <= _FP_W_TYPE_SIZE)                                       \
 641       r = X##_f[0];                                                     \
 642     else if (rsize <= 2*_FP_W_TYPE_SIZE)                                \
 643     {                                                                   \
 644       r = X##_f[1];                                                     \
 645       r <<= _FP_W_TYPE_SIZE;                                            \
 646       r += X##_f[0];                                                    \
 647     }                                                                   \
 648     else                                                                \
 649     {                                                                   \
 650       /* I'm feeling lazy so we deal with int == 3words (implausible)*/ \
 651       /* and int == 4words as a single case.                     */     \
 652       r = X##_f[3];                                                     \
 653       r <<= _FP_W_TYPE_SIZE;                                            \
 654       r += X##_f[2];                                                    \
 655       r <<= _FP_W_TYPE_SIZE;                                            \
 656       r += X##_f[1];                                                    \
 657       r <<= _FP_W_TYPE_SIZE;                                            \
 658       r += X##_f[0];                                                    \
 659     }                                                                   \
 660   } while (0)
 661
 662 /* "No disassemble Number Five!" */
 663 /* move an integer of size rsize into X's fractional part. We rely on
 664  * the _f[] array consisting of words of size _FP_W_TYPE_SIZE to avoid
 665  * having to mask the values we store into it.
 666  */
 667 #define _FP_FRAC_DISASSEMBLE_4(X, r, rsize)                             \
 668   do {                                                                  \
 669     X##_f[0] = r;                                                       \
 670     X##_f[1] = (rsize <= _FP_W_TYPE_SIZE ? 0 : r >> _FP_W_TYPE_SIZE);   \
 671     X##_f[2] = (rsize <= 2*_FP_W_TYPE_SIZE ? 0 : r >> 2*_FP_W_TYPE_SIZE); \
 672     X##_f[3] = (rsize <= 3*_FP_W_TYPE_SIZE ? 0 : r >> 3*_FP_W_TYPE_SIZE); \
 673   } while (0);
 674
 675 #define _FP_FRAC_COPY_4_1(D, S)                 \
 676 do {                                            \
 677   D##_f[0] = S##_f;                             \
 678   D##_f[1] = D##_f[2] = D##_f[3] = 0;           \
 679 } while (0)
 680
 681 #define _FP_FRAC_COPY_4_2(D, S)                 \
 682 do {                                            \
 683   D##_f[0] = S##_f0;                            \
 684   D##_f[1] = S##_f1;                            \
 685   D##_f[2] = D##_f[3] = 0;                      \
 686 } while (0)
 687
 688 #define _FP_FRAC_COPY_4_4(D,S)  _FP_FRAC_COPY_4(D,S)