soft-fp/op-2.h

   1 /* Software floating-point emulation.
   2    Basic two-word fraction declaration and manipulation.
   3    Copyright (C) 1997-2015 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, see
  31    <http://www.gnu.org/licenses/>.  */
  32
  33 #define _FP_FRAC_DECL_2(X)      _FP_W_TYPE X##_f0, X##_f1
  34 #define _FP_FRAC_COPY_2(D, S)   (D##_f0 = S##_f0, D##_f1 = S##_f1)
  35 #define _FP_FRAC_SET_2(X, I)    __FP_FRAC_SET_2 (X, I)
  36 #define _FP_FRAC_HIGH_2(X)      (X##_f1)
  37 #define _FP_FRAC_LOW_2(X)       (X##_f0)
  38 #define _FP_FRAC_WORD_2(X, w)   (X##_f##w)
  39
  40 #define _FP_FRAC_SLL_2(X, N)                                            \
  41   (void) (((N) < _FP_W_TYPE_SIZE)                                       \
  42           ? ({                                                          \
  43               if (__builtin_constant_p (N) && (N) == 1)                 \
  44                 {                                                       \
  45                   X##_f1 = X##_f1 + X##_f1 + (((_FP_WS_TYPE) (X##_f0)) < 0); \
  46                   X##_f0 += X##_f0;                                     \
  47                 }                                                       \
  48               else                                                      \
  49                 {                                                       \
  50                   X##_f1 = X##_f1 << (N) | X##_f0 >> (_FP_W_TYPE_SIZE - (N)); \
  51                   X##_f0 <<= (N);                                       \
  52                 }                                                       \
  53               0;                                                        \
  54             })                                                          \
  55           : ({                                                          \
  56               X##_f1 = X##_f0 << ((N) - _FP_W_TYPE_SIZE);               \
  57               X##_f0 = 0;                                               \
  58             }))
  59
  60
  61 #define _FP_FRAC_SRL_2(X, N)                                            \
  62   (void) (((N) < _FP_W_TYPE_SIZE)                                       \
  63           ? ({                                                          \
  64               X##_f0 = X##_f0 >> (N) | X##_f1 << (_FP_W_TYPE_SIZE - (N)); \
  65               X##_f1 >>= (N);                                           \
  66             })                                                          \
  67           : ({                                                          \
  68               X##_f0 = X##_f1 >> ((N) - _FP_W_TYPE_SIZE);               \
  69               X##_f1 = 0;                                               \
  70             }))
  71
  72 /* Right shift with sticky-lsb.  */
  73 #define _FP_FRAC_SRST_2(X, S, N, sz)                                    \
  74   (void) (((N) < _FP_W_TYPE_SIZE)                                       \
  75           ? ({                                                          \
  76               S = (__builtin_constant_p (N) && (N) == 1                 \
  77                    ? X##_f0 & 1                                         \
  78                    : (X##_f0 << (_FP_W_TYPE_SIZE - (N))) != 0);         \
  79               X##_f0 = (X##_f1 << (_FP_W_TYPE_SIZE - (N)) | X##_f0 >> (N)); \
  80               X##_f1 >>= (N);                                           \
  81             })                                                          \
  82           : ({                                                          \
  83               S = ((((N) == _FP_W_TYPE_SIZE                             \
  84                      ? 0                                                \
  85                      : (X##_f1 << (2*_FP_W_TYPE_SIZE - (N))))           \
  86                     | X##_f0) != 0);                                    \
  87               X##_f0 = (X##_f1 >> ((N) - _FP_W_TYPE_SIZE));             \
  88               X##_f1 = 0;                                               \
  89             }))
  90
  91 #define _FP_FRAC_SRS_2(X, N, sz)                                        \
  92   (void) (((N) < _FP_W_TYPE_SIZE)                                       \
  93           ? ({                                                          \
  94               X##_f0 = (X##_f1 << (_FP_W_TYPE_SIZE - (N)) | X##_f0 >> (N) \
  95                         | (__builtin_constant_p (N) && (N) == 1         \
  96                            ? X##_f0 & 1                                 \
  97                            : (X##_f0 << (_FP_W_TYPE_SIZE - (N))) != 0)); \
  98               X##_f1 >>= (N);                                           \
  99             })                                                          \
 100           : ({                                                          \
 101               X##_f0 = (X##_f1 >> ((N) - _FP_W_TYPE_SIZE)               \
 102                         | ((((N) == _FP_W_TYPE_SIZE                     \
 103                              ? 0                                        \
 104                              : (X##_f1 << (2*_FP_W_TYPE_SIZE - (N))))   \
 105                             | X##_f0) != 0));                           \
 106               X##_f1 = 0;                                               \
 107             }))
 108
 109 #define _FP_FRAC_ADDI_2(X, I)   \
 110   __FP_FRAC_ADDI_2 (X##_f1, X##_f0, I)
 111
 112 #define _FP_FRAC_ADD_2(R, X, Y) \
 113   __FP_FRAC_ADD_2 (R##_f1, R##_f0, X##_f1, X##_f0, Y##_f1, Y##_f0)
 114
 115 #define _FP_FRAC_SUB_2(R, X, Y) \
 116   __FP_FRAC_SUB_2 (R##_f1, R##_f0, X##_f1, X##_f0, Y##_f1, Y##_f0)
 117
 118 #define _FP_FRAC_DEC_2(X, Y)    \
 119   __FP_FRAC_DEC_2 (X##_f1, X##_f0, Y##_f1, Y##_f0)
 120
 121 #define _FP_FRAC_CLZ_2(R, X)                    \
 122   do                                            \
 123     {                                           \
 124       if (X##_f1)                               \
 125         __FP_CLZ ((R), X##_f1);                 \
 126       else                                      \
 127         {                                       \
 128           __FP_CLZ ((R), X##_f0);               \
 129           (R) += _FP_W_TYPE_SIZE;               \
 130         }                                       \
 131     }                                           \
 132   while (0)
 133
 134 /* Predicates.  */
 135 #define _FP_FRAC_NEGP_2(X)      ((_FP_WS_TYPE) X##_f1 < 0)
 136 #define _FP_FRAC_ZEROP_2(X)     ((X##_f1 | X##_f0) == 0)
 137 #define _FP_FRAC_OVERP_2(fs, X) (_FP_FRAC_HIGH_##fs (X) & _FP_OVERFLOW_##fs)
 138 #define _FP_FRAC_CLEAR_OVERP_2(fs, X)   (_FP_FRAC_HIGH_##fs (X) &= ~_FP_OVERFLOW_##fs)
 139 #define _FP_FRAC_HIGHBIT_DW_2(fs, X)    \
 140   (_FP_FRAC_HIGH_DW_##fs (X) & _FP_HIGHBIT_DW_##fs)
 141 #define _FP_FRAC_EQ_2(X, Y)     (X##_f1 == Y##_f1 && X##_f0 == Y##_f0)
 142 #define _FP_FRAC_GT_2(X, Y)     \
 143   (X##_f1 > Y##_f1 || (X##_f1 == Y##_f1 && X##_f0 > Y##_f0))
 144 #define _FP_FRAC_GE_2(X, Y)     \
 145   (X##_f1 > Y##_f1 || (X##_f1 == Y##_f1 && X##_f0 >= Y##_f0))
 146
 147 #define _FP_ZEROFRAC_2          0, 0
 148 #define _FP_MINFRAC_2           0, 1
 149 #define _FP_MAXFRAC_2           (~(_FP_WS_TYPE) 0), (~(_FP_WS_TYPE) 0)
 150
 151 /* Internals.  */
 152
 153 #define __FP_FRAC_SET_2(X, I1, I0)      (X##_f0 = I0, X##_f1 = I1)
 154
 155 #define __FP_CLZ_2(R, xh, xl)                   \
 156   do                                            \
 157     {                                           \
 158       if (xh)                                   \
 159         __FP_CLZ ((R), xh);                     \
 160       else                                      \
 161         {                                       \
 162           __FP_CLZ ((R), xl);                   \
 163           (R) += _FP_W_TYPE_SIZE;               \
 164         }                                       \
 165     }                                           \
 166   while (0)
 167
 168 #if 0
 169
 170 # ifndef __FP_FRAC_ADDI_2
 171 #  define __FP_FRAC_ADDI_2(xh, xl, i)   \
 172   (xh += ((xl += i) < i))
 173 # endif
 174 # ifndef __FP_FRAC_ADD_2
 175 #  define __FP_FRAC_ADD_2(rh, rl, xh, xl, yh, yl)       \
 176   (rh = xh + yh + ((rl = xl + yl) < xl))
 177 # endif
 178 # ifndef __FP_FRAC_SUB_2
 179 #  define __FP_FRAC_SUB_2(rh, rl, xh, xl, yh, yl)       \
 180   (rh = xh - yh - ((rl = xl - yl) > xl))
 181 # endif
 182 # ifndef __FP_FRAC_DEC_2
 183 #  define __FP_FRAC_DEC_2(xh, xl, yh, yl)               \
 184   do                                                    \
 185     {                                                   \
 186       UWtype __FP_FRAC_DEC_2_t = xl;                    \
 187       xh -= yh + ((xl -= yl) > __FP_FRAC_DEC_2_t);      \
 188     }                                                   \
 189   while (0)
 190 # endif
 191
 192 #else
 193
 194 # undef __FP_FRAC_ADDI_2
 195 # define __FP_FRAC_ADDI_2(xh, xl, i)    add_ssaaaa (xh, xl, xh, xl, 0, i)
 196 # undef __FP_FRAC_ADD_2
 197 # define __FP_FRAC_ADD_2                add_ssaaaa
 198 # undef __FP_FRAC_SUB_2
 199 # define __FP_FRAC_SUB_2                sub_ddmmss
 200 # undef __FP_FRAC_DEC_2
 201 # define __FP_FRAC_DEC_2(xh, xl, yh, yl)        \
 202   sub_ddmmss (xh, xl, xh, xl, yh, yl)
 203
 204 #endif
 205
 206 /* Unpack the raw bits of a native fp value.  Do not classify or
 207    normalize the data.  */
 208
 209 #define _FP_UNPACK_RAW_2(fs, X, val)                    \
 210   do                                                    \
 211     {                                                   \
 212       union _FP_UNION_##fs _FP_UNPACK_RAW_2_flo;        \
 213       _FP_UNPACK_RAW_2_flo.flt = (val);                 \
 214                                                         \
 215       X##_f0 = _FP_UNPACK_RAW_2_flo.bits.frac0;         \
 216       X##_f1 = _FP_UNPACK_RAW_2_flo.bits.frac1;         \
 217       X##_e  = _FP_UNPACK_RAW_2_flo.bits.exp;           \
 218       X##_s  = _FP_UNPACK_RAW_2_flo.bits.sign;          \
 219     }                                                   \
 220   while (0)
 221
 222 #define _FP_UNPACK_RAW_2_P(fs, X, val)                  \
 223   do                                                    \
 224     {                                                   \
 225       union _FP_UNION_##fs *_FP_UNPACK_RAW_2_P_flo      \
 226         = (union _FP_UNION_##fs *) (val);               \
 227                                                         \
 228       X##_f0 = _FP_UNPACK_RAW_2_P_flo->bits.frac0;      \
 229       X##_f1 = _FP_UNPACK_RAW_2_P_flo->bits.frac1;      \
 230       X##_e  = _FP_UNPACK_RAW_2_P_flo->bits.exp;        \
 231       X##_s  = _FP_UNPACK_RAW_2_P_flo->bits.sign;       \
 232     }                                                   \
 233   while (0)
 234
 235
 236 /* Repack the raw bits of a native fp value.  */
 237
 238 #define _FP_PACK_RAW_2(fs, val, X)              \
 239   do                                            \
 240     {                                           \
 241       union _FP_UNION_##fs _FP_PACK_RAW_2_flo;  \
 242                                                 \
 243       _FP_PACK_RAW_2_flo.bits.frac0 = X##_f0;   \
 244       _FP_PACK_RAW_2_flo.bits.frac1 = X##_f1;   \
 245       _FP_PACK_RAW_2_flo.bits.exp   = X##_e;    \
 246       _FP_PACK_RAW_2_flo.bits.sign  = X##_s;    \
 247                                                 \
 248       (val) = _FP_PACK_RAW_2_flo.flt;           \
 249     }                                           \
 250   while (0)
 251
 252 #define _FP_PACK_RAW_2_P(fs, val, X)                    \
 253   do                                                    \
 254     {                                                   \
 255       union _FP_UNION_##fs *_FP_PACK_RAW_2_P_flo        \
 256         = (union _FP_UNION_##fs *) (val);               \
 257                                                         \
 258       _FP_PACK_RAW_2_P_flo->bits.frac0 = X##_f0;        \
 259       _FP_PACK_RAW_2_P_flo->bits.frac1 = X##_f1;        \
 260       _FP_PACK_RAW_2_P_flo->bits.exp   = X##_e;         \
 261       _FP_PACK_RAW_2_P_flo->bits.sign  = X##_s;         \
 262     }                                                   \
 263   while (0)
 264
 265
 266 /* Multiplication algorithms: */
 267
 268 /* Given a 1W * 1W => 2W primitive, do the extended multiplication.  */
 269
 270 #define _FP_MUL_MEAT_DW_2_wide(wfracbits, R, X, Y, doit)                \
 271   do                                                                    \
 272     {                                                                   \
 273       _FP_FRAC_DECL_2 (_FP_MUL_MEAT_DW_2_wide_b);                       \
 274       _FP_FRAC_DECL_2 (_FP_MUL_MEAT_DW_2_wide_c);                       \
 275                                                                         \
 276       doit (_FP_FRAC_WORD_4 (R, 1), _FP_FRAC_WORD_4 (R, 0),             \
 277             X##_f0, Y##_f0);                                            \
 278       doit (_FP_MUL_MEAT_DW_2_wide_b_f1, _FP_MUL_MEAT_DW_2_wide_b_f0,   \
 279             X##_f0, Y##_f1);                                            \
 280       doit (_FP_MUL_MEAT_DW_2_wide_c_f1, _FP_MUL_MEAT_DW_2_wide_c_f0,   \
 281             X##_f1, Y##_f0);                                            \
 282       doit (_FP_FRAC_WORD_4 (R, 3), _FP_FRAC_WORD_4 (R, 2),             \
 283             X##_f1, Y##_f1);                                            \
 284                                                                         \
 285       __FP_FRAC_ADD_3 (_FP_FRAC_WORD_4 (R, 3), _FP_FRAC_WORD_4 (R, 2),  \
 286                        _FP_FRAC_WORD_4 (R, 1), 0,                       \
 287                        _FP_MUL_MEAT_DW_2_wide_b_f1,                     \
 288                        _FP_MUL_MEAT_DW_2_wide_b_f0,                     \
 289                        _FP_FRAC_WORD_4 (R, 3), _FP_FRAC_WORD_4 (R, 2),  \
 290                        _FP_FRAC_WORD_4 (R, 1));                         \
 291       __FP_FRAC_ADD_3 (_FP_FRAC_WORD_4 (R, 3), _FP_FRAC_WORD_4 (R, 2),  \
 292                        _FP_FRAC_WORD_4 (R, 1), 0,                       \
 293                        _FP_MUL_MEAT_DW_2_wide_c_f1,                     \
 294                        _FP_MUL_MEAT_DW_2_wide_c_f0,                     \
 295                        _FP_FRAC_WORD_4 (R, 3), _FP_FRAC_WORD_4 (R, 2),  \
 296                        _FP_FRAC_WORD_4 (R, 1));                         \
 297     }                                                                   \
 298   while (0)
 299
 300 #define _FP_MUL_MEAT_2_wide(wfracbits, R, X, Y, doit)                   \
 301   do                                                                    \
 302     {                                                                   \
 303       _FP_FRAC_DECL_4 (_FP_MUL_MEAT_2_wide_z);                          \
 304                                                                         \
 305       _FP_MUL_MEAT_DW_2_wide ((wfracbits), _FP_MUL_MEAT_2_wide_z,       \
 306                               X, Y, doit);                              \
 307                                                                         \
 308       /* Normalize since we know where the msb of the multiplicands     \
 309          were (bit B), we know that the msb of the of the product is    \
 310          at either 2B or 2B-1.  */                                      \
 311       _FP_FRAC_SRS_4 (_FP_MUL_MEAT_2_wide_z, (wfracbits)-1,             \
 312                       2*(wfracbits));                                   \
 313       R##_f0 = _FP_FRAC_WORD_4 (_FP_MUL_MEAT_2_wide_z, 0);              \
 314       R##_f1 = _FP_FRAC_WORD_4 (_FP_MUL_MEAT_2_wide_z, 1);              \
 315     }                                                                   \
 316   while (0)
 317
 318 /* Given a 1W * 1W => 2W primitive, do the extended multiplication.
 319    Do only 3 multiplications instead of four. This one is for machines
 320    where multiplication is much more expensive than subtraction.  */
 321
 322 #define _FP_MUL_MEAT_DW_2_wide_3mul(wfracbits, R, X, Y, doit)           \
 323   do                                                                    \
 324     {                                                                   \
 325       _FP_FRAC_DECL_2 (_FP_MUL_MEAT_DW_2_wide_3mul_b);                  \
 326       _FP_FRAC_DECL_2 (_FP_MUL_MEAT_DW_2_wide_3mul_c);                  \
 327       _FP_W_TYPE _FP_MUL_MEAT_DW_2_wide_3mul_d;                         \
 328       int _FP_MUL_MEAT_DW_2_wide_3mul_c1;                               \
 329       int _FP_MUL_MEAT_DW_2_wide_3mul_c2;                               \
 330                                                                         \
 331       _FP_MUL_MEAT_DW_2_wide_3mul_b_f0 = X##_f0 + X##_f1;               \
 332       _FP_MUL_MEAT_DW_2_wide_3mul_c1                                    \
 333         = _FP_MUL_MEAT_DW_2_wide_3mul_b_f0 < X##_f0;                    \
 334       _FP_MUL_MEAT_DW_2_wide_3mul_b_f1 = Y##_f0 + Y##_f1;               \
 335       _FP_MUL_MEAT_DW_2_wide_3mul_c2                                    \
 336         = _FP_MUL_MEAT_DW_2_wide_3mul_b_f1 < Y##_f0;                    \
 337       doit (_FP_MUL_MEAT_DW_2_wide_3mul_d, _FP_FRAC_WORD_4 (R, 0),      \
 338             X##_f0, Y##_f0);                                            \
 339       doit (_FP_FRAC_WORD_4 (R, 2), _FP_FRAC_WORD_4 (R, 1),             \
 340             _FP_MUL_MEAT_DW_2_wide_3mul_b_f0,                           \
 341             _FP_MUL_MEAT_DW_2_wide_3mul_b_f1);                          \
 342       doit (_FP_MUL_MEAT_DW_2_wide_3mul_c_f1,                           \
 343             _FP_MUL_MEAT_DW_2_wide_3mul_c_f0, X##_f1, Y##_f1);          \
 344                                                                         \
 345       _FP_MUL_MEAT_DW_2_wide_3mul_b_f0                                  \
 346         &= -_FP_MUL_MEAT_DW_2_wide_3mul_c2;                             \
 347       _FP_MUL_MEAT_DW_2_wide_3mul_b_f1                                  \
 348         &= -_FP_MUL_MEAT_DW_2_wide_3mul_c1;                             \
 349       __FP_FRAC_ADD_3 (_FP_FRAC_WORD_4 (R, 3), _FP_FRAC_WORD_4 (R, 2),  \
 350                        _FP_FRAC_WORD_4 (R, 1),                          \
 351                        (_FP_MUL_MEAT_DW_2_wide_3mul_c1                  \
 352                         & _FP_MUL_MEAT_DW_2_wide_3mul_c2), 0,           \
 353                        _FP_MUL_MEAT_DW_2_wide_3mul_d,                   \
 354                        0, _FP_FRAC_WORD_4 (R, 2), _FP_FRAC_WORD_4 (R, 1)); \
 355       __FP_FRAC_ADDI_2 (_FP_FRAC_WORD_4 (R, 3), _FP_FRAC_WORD_4 (R, 2), \
 356                         _FP_MUL_MEAT_DW_2_wide_3mul_b_f0);              \
 357       __FP_FRAC_ADDI_2 (_FP_FRAC_WORD_4 (R, 3), _FP_FRAC_WORD_4 (R, 2), \
 358                         _FP_MUL_MEAT_DW_2_wide_3mul_b_f1);              \
 359       __FP_FRAC_DEC_3 (_FP_FRAC_WORD_4 (R, 3), _FP_FRAC_WORD_4 (R, 2),  \
 360                        _FP_FRAC_WORD_4 (R, 1),                          \
 361                        0, _FP_MUL_MEAT_DW_2_wide_3mul_d,                \
 362                        _FP_FRAC_WORD_4 (R, 0));                         \
 363       __FP_FRAC_DEC_3 (_FP_FRAC_WORD_4 (R, 3), _FP_FRAC_WORD_4 (R, 2),  \
 364                        _FP_FRAC_WORD_4 (R, 1), 0,                       \
 365                        _FP_MUL_MEAT_DW_2_wide_3mul_c_f1,                \
 366                        _FP_MUL_MEAT_DW_2_wide_3mul_c_f0);               \
 367       __FP_FRAC_ADD_2 (_FP_FRAC_WORD_4 (R, 3), _FP_FRAC_WORD_4 (R, 2),  \
 368                        _FP_MUL_MEAT_DW_2_wide_3mul_c_f1,                \
 369                        _FP_MUL_MEAT_DW_2_wide_3mul_c_f0,                \
 370                        _FP_FRAC_WORD_4 (R, 3), _FP_FRAC_WORD_4 (R, 2)); \
 371     }                                                                   \
 372   while (0)
 373
 374 #define _FP_MUL_MEAT_2_wide_3mul(wfracbits, R, X, Y, doit)              \
 375   do                                                                    \
 376     {                                                                   \
 377       _FP_FRAC_DECL_4 (_FP_MUL_MEAT_2_wide_3mul_z);                     \
 378                                                                         \
 379       _FP_MUL_MEAT_DW_2_wide_3mul ((wfracbits),                         \
 380                                    _FP_MUL_MEAT_2_wide_3mul_z,          \
 381                                    X, Y, doit);                         \
 382                                                                         \
 383       /* Normalize since we know where the msb of the multiplicands     \
 384          were (bit B), we know that the msb of the of the product is    \
 385          at either 2B or 2B-1.  */                                      \
 386       _FP_FRAC_SRS_4 (_FP_MUL_MEAT_2_wide_3mul_z,                       \
 387                       (wfracbits)-1, 2*(wfracbits));                    \
 388       R##_f0 = _FP_FRAC_WORD_4 (_FP_MUL_MEAT_2_wide_3mul_z, 0);         \
 389       R##_f1 = _FP_FRAC_WORD_4 (_FP_MUL_MEAT_2_wide_3mul_z, 1);         \
 390     }                                                                   \
 391   while (0)
 392
 393 #define _FP_MUL_MEAT_DW_2_gmp(wfracbits, R, X, Y)       \
 394   do                                                    \
 395     {                                                   \
 396       _FP_W_TYPE _FP_MUL_MEAT_DW_2_gmp_x[2];            \
 397       _FP_W_TYPE _FP_MUL_MEAT_DW_2_gmp_y[2];            \
 398       _FP_MUL_MEAT_DW_2_gmp_x[0] = X##_f0;              \
 399       _FP_MUL_MEAT_DW_2_gmp_x[1] = X##_f1;              \
 400       _FP_MUL_MEAT_DW_2_gmp_y[0] = Y##_f0;              \
 401       _FP_MUL_MEAT_DW_2_gmp_y[1] = Y##_f1;              \
 402                                                         \
 403       mpn_mul_n (R##_f, _FP_MUL_MEAT_DW_2_gmp_x,        \
 404                  _FP_MUL_MEAT_DW_2_gmp_y, 2);           \
 405     }                                                   \
 406   while (0)
 407
 408 #define _FP_MUL_MEAT_2_gmp(wfracbits, R, X, Y)                          \
 409   do                                                                    \
 410     {                                                                   \
 411       _FP_FRAC_DECL_4 (_FP_MUL_MEAT_2_gmp_z);                           \
 412                                                                         \
 413       _FP_MUL_MEAT_DW_2_gmp ((wfracbits), _FP_MUL_MEAT_2_gmp_z, X, Y);  \
 414                                                                         \
 415       /* Normalize since we know where the msb of the multiplicands     \
 416          were (bit B), we know that the msb of the of the product is    \
 417          at either 2B or 2B-1.  */                                      \
 418       _FP_FRAC_SRS_4 (_FP_MUL_MEAT_2_gmp_z, (wfracbits)-1,              \
 419                       2*(wfracbits));                                   \
 420       R##_f0 = _FP_MUL_MEAT_2_gmp_z_f[0];                               \
 421       R##_f1 = _FP_MUL_MEAT_2_gmp_z_f[1];                               \
 422     }                                                                   \
 423   while (0)
 424
 425 /* Do at most 120x120=240 bits multiplication using double floating
 426    point multiplication.  This is useful if floating point
 427    multiplication has much bigger throughput than integer multiply.
 428    It is supposed to work for _FP_W_TYPE_SIZE 64 and wfracbits
 429    between 106 and 120 only.
 430    Caller guarantees that X and Y has (1LLL << (wfracbits - 1)) set.
 431    SETFETZ is a macro which will disable all FPU exceptions and set rounding
 432    towards zero,  RESETFE should optionally reset it back.  */
 433
 434 #define _FP_MUL_MEAT_2_120_240_double(wfracbits, R, X, Y, setfetz, resetfe) \
 435   do                                                                    \
 436     {                                                                   \
 437       static const double _const[] =                                    \
 438         {                                                               \
 439           /* 2^-24 */ 5.9604644775390625e-08,                           \
 440           /* 2^-48 */ 3.5527136788005009e-15,                           \
 441           /* 2^-72 */ 2.1175823681357508e-22,                           \
 442           /* 2^-96 */ 1.2621774483536189e-29,                           \
 443           /* 2^28 */ 2.68435456e+08,                                    \
 444           /* 2^4 */ 1.600000e+01,                                       \
 445           /* 2^-20 */ 9.5367431640625e-07,                              \
 446           /* 2^-44 */ 5.6843418860808015e-14,                           \
 447           /* 2^-68 */ 3.3881317890172014e-21,                           \
 448           /* 2^-92 */ 2.0194839173657902e-28,                           \
 449           /* 2^-116 */ 1.2037062152420224e-35                           \
 450         };                                                              \
 451       double _a240, _b240, _c240, _d240, _e240, _f240,                  \
 452         _g240, _h240, _i240, _j240, _k240;                              \
 453       union { double d; UDItype i; } _l240, _m240, _n240, _o240,        \
 454                                        _p240, _q240, _r240, _s240;      \
 455       UDItype _t240, _u240, _v240, _w240, _x240, _y240 = 0;             \
 456                                                                         \
 457       if ((wfracbits) < 106 || (wfracbits) > 120)                       \
 458         abort ();                                                       \
 459                                                                         \
 460       setfetz;                                                          \
 461                                                                         \
 462       _e240 = (double) (long) (X##_f0 & 0xffffff);                      \
 463       _j240 = (double) (long) (Y##_f0 & 0xffffff);                      \
 464       _d240 = (double) (long) ((X##_f0 >> 24) & 0xffffff);              \
 465       _i240 = (double) (long) ((Y##_f0 >> 24) & 0xffffff);              \
 466       _c240 = (double) (long) (((X##_f1 << 16) & 0xffffff) | (X##_f0 >> 48)); \
 467       _h240 = (double) (long) (((Y##_f1 << 16) & 0xffffff) | (Y##_f0 >> 48)); \
 468       _b240 = (double) (long) ((X##_f1 >> 8) & 0xffffff);               \
 469       _g240 = (double) (long) ((Y##_f1 >> 8) & 0xffffff);               \
 470       _a240 = (double) (long) (X##_f1 >> 32);                           \
 471       _f240 = (double) (long) (Y##_f1 >> 32);                           \
 472       _e240 *= _const[3];                                               \
 473       _j240 *= _const[3];                                               \
 474       _d240 *= _const[2];                                               \
 475       _i240 *= _const[2];                                               \
 476       _c240 *= _const[1];                                               \
 477       _h240 *= _const[1];                                               \
 478       _b240 *= _const[0];                                               \
 479       _g240 *= _const[0];                                               \
 480       _s240.d =                                                       _e240*_j240; \
 481       _r240.d =                                         _d240*_j240 + _e240*_i240; \
 482       _q240.d =                           _c240*_j240 + _d240*_i240 + _e240*_h240; \
 483       _p240.d =             _b240*_j240 + _c240*_i240 + _d240*_h240 + _e240*_g240; \
 484       _o240.d = _a240*_j240 + _b240*_i240 + _c240*_h240 + _d240*_g240 + _e240*_f240; \
 485       _n240.d = _a240*_i240 + _b240*_h240 + _c240*_g240 + _d240*_f240;  \
 486       _m240.d = _a240*_h240 + _b240*_g240 + _c240*_f240;                \
 487       _l240.d = _a240*_g240 + _b240*_f240;                              \
 488       _k240 =   _a240*_f240;                                            \
 489       _r240.d += _s240.d;                                               \
 490       _q240.d += _r240.d;                                               \
 491       _p240.d += _q240.d;                                               \
 492       _o240.d += _p240.d;                                               \
 493       _n240.d += _o240.d;                                               \
 494       _m240.d += _n240.d;                                               \
 495       _l240.d += _m240.d;                                               \
 496       _k240 += _l240.d;                                                 \
 497       _s240.d -= ((_const[10]+_s240.d)-_const[10]);                     \
 498       _r240.d -= ((_const[9]+_r240.d)-_const[9]);                       \
 499       _q240.d -= ((_const[8]+_q240.d)-_const[8]);                       \
 500       _p240.d -= ((_const[7]+_p240.d)-_const[7]);                       \
 501       _o240.d += _const[7];                                             \
 502       _n240.d += _const[6];                                             \
 503       _m240.d += _const[5];                                             \
 504       _l240.d += _const[4];                                             \
 505       if (_s240.d != 0.0)                                               \
 506         _y240 = 1;                                                      \
 507       if (_r240.d != 0.0)                                               \
 508         _y240 = 1;                                                      \
 509       if (_q240.d != 0.0)                                               \
 510         _y240 = 1;                                                      \
 511       if (_p240.d != 0.0)                                               \
 512         _y240 = 1;                                                      \
 513       _t240 = (DItype) _k240;                                           \
 514       _u240 = _l240.i;                                                  \
 515       _v240 = _m240.i;                                                  \
 516       _w240 = _n240.i;                                                  \
 517       _x240 = _o240.i;                                                  \
 518       R##_f1 = ((_t240 << (128 - (wfracbits - 1)))                      \
 519                 | ((_u240 & 0xffffff) >> ((wfracbits - 1) - 104)));     \
 520       R##_f0 = (((_u240 & 0xffffff) << (168 - (wfracbits - 1)))         \
 521                 | ((_v240 & 0xffffff) << (144 - (wfracbits - 1)))       \
 522                 | ((_w240 & 0xffffff) << (120 - (wfracbits - 1)))       \
 523                 | ((_x240 & 0xffffff) >> ((wfracbits - 1) - 96))        \
 524                 | _y240);                                               \
 525       resetfe;                                                          \
 526     }                                                                   \
 527   while (0)
 528
 529 /* Division algorithms: */
 530
 531 #define _FP_DIV_MEAT_2_udiv(fs, R, X, Y)                                \
 532   do                                                                    \
 533     {                                                                   \
 534       _FP_W_TYPE _FP_DIV_MEAT_2_udiv_n_f2;                              \
 535       _FP_W_TYPE _FP_DIV_MEAT_2_udiv_n_f1;                              \
 536       _FP_W_TYPE _FP_DIV_MEAT_2_udiv_n_f0;                              \
 537       _FP_W_TYPE _FP_DIV_MEAT_2_udiv_r_f1;                              \
 538       _FP_W_TYPE _FP_DIV_MEAT_2_udiv_r_f0;                              \
 539       _FP_W_TYPE _FP_DIV_MEAT_2_udiv_m_f1;                              \
 540       _FP_W_TYPE _FP_DIV_MEAT_2_udiv_m_f0;                              \
 541       if (_FP_FRAC_GE_2 (X, Y))                                         \
 542         {                                                               \
 543           _FP_DIV_MEAT_2_udiv_n_f2 = X##_f1 >> 1;                       \
 544           _FP_DIV_MEAT_2_udiv_n_f1                                      \
 545             = X##_f1 << (_FP_W_TYPE_SIZE - 1) | X##_f0 >> 1;            \
 546           _FP_DIV_MEAT_2_udiv_n_f0                                      \
 547             = X##_f0 << (_FP_W_TYPE_SIZE - 1);                          \
 548         }                                                               \
 549       else                                                              \
 550         {                                                               \
 551           R##_e--;                                                      \
 552           _FP_DIV_MEAT_2_udiv_n_f2 = X##_f1;                            \
 553           _FP_DIV_MEAT_2_udiv_n_f1 = X##_f0;                            \
 554           _FP_DIV_MEAT_2_udiv_n_f0 = 0;                                 \
 555         }                                                               \
 556                                                                         \
 557       /* Normalize, i.e. make the most significant bit of the           \
 558          denominator set.  */                                           \
 559       _FP_FRAC_SLL_2 (Y, _FP_WFRACXBITS_##fs);                          \
 560                                                                         \
 561       udiv_qrnnd (R##_f1, _FP_DIV_MEAT_2_udiv_r_f1,                     \
 562                   _FP_DIV_MEAT_2_udiv_n_f2, _FP_DIV_MEAT_2_udiv_n_f1,   \
 563                   Y##_f1);                                              \
 564       umul_ppmm (_FP_DIV_MEAT_2_udiv_m_f1, _FP_DIV_MEAT_2_udiv_m_f0,    \
 565                  R##_f1, Y##_f0);                                       \
 566       _FP_DIV_MEAT_2_udiv_r_f0 = _FP_DIV_MEAT_2_udiv_n_f0;              \
 567       if (_FP_FRAC_GT_2 (_FP_DIV_MEAT_2_udiv_m, _FP_DIV_MEAT_2_udiv_r)) \
 568         {                                                               \
 569           R##_f1--;                                                     \
 570           _FP_FRAC_ADD_2 (_FP_DIV_MEAT_2_udiv_r, Y,                     \
 571                           _FP_DIV_MEAT_2_udiv_r);                       \
 572           if (_FP_FRAC_GE_2 (_FP_DIV_MEAT_2_udiv_r, Y)                  \
 573               && _FP_FRAC_GT_2 (_FP_DIV_MEAT_2_udiv_m,                  \
 574                                 _FP_DIV_MEAT_2_udiv_r))                 \
 575             {                                                           \
 576               R##_f1--;                                                 \
 577               _FP_FRAC_ADD_2 (_FP_DIV_MEAT_2_udiv_r, Y,                 \
 578                               _FP_DIV_MEAT_2_udiv_r);                   \
 579             }                                                           \
 580         }                                                               \
 581       _FP_FRAC_DEC_2 (_FP_DIV_MEAT_2_udiv_r, _FP_DIV_MEAT_2_udiv_m);    \
 582                                                                         \
 583       if (_FP_DIV_MEAT_2_udiv_r_f1 == Y##_f1)                           \
 584         {                                                               \
 585           /* This is a special case, not an optimization                \
 586              (_FP_DIV_MEAT_2_udiv_r/Y##_f1 would not fit into UWtype).  \
 587              As _FP_DIV_MEAT_2_udiv_r is guaranteed to be < Y,          \
 588              R##_f0 can be either (UWtype)-1 or (UWtype)-2.  But as we  \
 589              know what kind of bits it is (sticky, guard, round),       \
 590              we don't care.  We also don't care what the reminder is,   \
 591              because the guard bit will be set anyway.  -jj */          \
 592           R##_f0 = -1;                                                  \
 593         }                                                               \
 594       else                                                              \
 595         {                                                               \
 596           udiv_qrnnd (R##_f0, _FP_DIV_MEAT_2_udiv_r_f1,                 \
 597                       _FP_DIV_MEAT_2_udiv_r_f1,                         \
 598                       _FP_DIV_MEAT_2_udiv_r_f0, Y##_f1);                \
 599           umul_ppmm (_FP_DIV_MEAT_2_udiv_m_f1,                          \
 600                      _FP_DIV_MEAT_2_udiv_m_f0, R##_f0, Y##_f0);         \
 601           _FP_DIV_MEAT_2_udiv_r_f0 = 0;                                 \
 602           if (_FP_FRAC_GT_2 (_FP_DIV_MEAT_2_udiv_m,                     \
 603                              _FP_DIV_MEAT_2_udiv_r))                    \
 604             {                                                           \
 605               R##_f0--;                                                 \
 606               _FP_FRAC_ADD_2 (_FP_DIV_MEAT_2_udiv_r, Y,                 \
 607                               _FP_DIV_MEAT_2_udiv_r);                   \
 608               if (_FP_FRAC_GE_2 (_FP_DIV_MEAT_2_udiv_r, Y)              \
 609                   && _FP_FRAC_GT_2 (_FP_DIV_MEAT_2_udiv_m,              \
 610                                     _FP_DIV_MEAT_2_udiv_r))             \
 611                 {                                                       \
 612                   R##_f0--;                                             \
 613                   _FP_FRAC_ADD_2 (_FP_DIV_MEAT_2_udiv_r, Y,             \
 614                                   _FP_DIV_MEAT_2_udiv_r);               \
 615                 }                                                       \
 616             }                                                           \
 617           if (!_FP_FRAC_EQ_2 (_FP_DIV_MEAT_2_udiv_r,                    \
 618                               _FP_DIV_MEAT_2_udiv_m))                   \
 619             R##_f0 |= _FP_WORK_STICKY;                                  \
 620         }                                                               \
 621     }                                                                   \
 622   while (0)
 623
 624
 625 /* Square root algorithms:
 626    We have just one right now, maybe Newton approximation
 627    should be added for those machines where division is fast.  */
 628
 629 #define _FP_SQRT_MEAT_2(R, S, T, X, q)                          \
 630   do                                                            \
 631     {                                                           \
 632       while (q)                                                 \
 633         {                                                       \
 634           T##_f1 = S##_f1 + (q);                                \
 635           if (T##_f1 <= X##_f1)                                 \
 636             {                                                   \
 637               S##_f1 = T##_f1 + (q);                            \
 638               X##_f1 -= T##_f1;                                 \
 639               R##_f1 += (q);                                    \
 640             }                                                   \
 641           _FP_FRAC_SLL_2 (X, 1);                                \
 642           (q) >>= 1;                                            \
 643         }                                                       \
 644       (q) = (_FP_W_TYPE) 1 << (_FP_W_TYPE_SIZE - 1);            \
 645       while ((q) != _FP_WORK_ROUND)                             \
 646         {                                                       \
 647           T##_f0 = S##_f0 + (q);                                \
 648           T##_f1 = S##_f1;                                      \
 649           if (T##_f1 < X##_f1                                   \
 650               || (T##_f1 == X##_f1 && T##_f0 <= X##_f0))        \
 651             {                                                   \
 652               S##_f0 = T##_f0 + (q);                            \
 653               S##_f1 += (T##_f0 > S##_f0);                      \
 654               _FP_FRAC_DEC_2 (X, T);                            \
 655               R##_f0 += (q);                                    \
 656             }                                                   \
 657           _FP_FRAC_SLL_2 (X, 1);                                \
 658           (q) >>= 1;                                            \
 659         }                                                       \
 660       if (X##_f0 | X##_f1)                                      \
 661         {                                                       \
 662           if (S##_f1 < X##_f1                                   \
 663               || (S##_f1 == X##_f1 && S##_f0 < X##_f0))         \
 664             R##_f0 |= _FP_WORK_ROUND;                           \
 665           R##_f0 |= _FP_WORK_STICKY;                            \
 666         }                                                       \
 667     }                                                           \
 668   while (0)
 669
 670
 671 /* Assembly/disassembly for converting to/from integral types.
 672    No shifting or overflow handled here.  */
 673
 674 #define _FP_FRAC_ASSEMBLE_2(r, X, rsize)        \
 675   (void) (((rsize) <= _FP_W_TYPE_SIZE)          \
 676           ? ({ (r) = X##_f0; })                 \
 677           : ({                                  \
 678               (r) = X##_f1;                     \
 679               (r) <<= _FP_W_TYPE_SIZE;          \
 680               (r) += X##_f0;                    \
 681             }))
 682
 683 #define _FP_FRAC_DISASSEMBLE_2(X, r, rsize)     \
 684   do                                            \
 685     {                                           \
 686       X##_f0 = (r);                             \
 687       X##_f1 = ((rsize) <= _FP_W_TYPE_SIZE      \
 688                 ? 0                             \
 689                 : (r) >> _FP_W_TYPE_SIZE);      \
 690     }                                           \
 691   while (0)
 692
 693 /* Convert FP values between word sizes.  */
 694
 695 #define _FP_FRAC_COPY_1_2(D, S)         (D##_f = S##_f0)
 696
 697 #define _FP_FRAC_COPY_2_1(D, S)         ((D##_f0 = S##_f), (D##_f1 = 0))
 698
 699 #define _FP_FRAC_COPY_2_2(D, S)         _FP_FRAC_COPY_2 (D, S)