include/math-emu/extended.h

   1 /* Software floating-point emulation.
   2    Definitions for IEEE Extended Precision.
   3    Copyright (C) 1999 Free Software Foundation, Inc.
   4    This file is part of the GNU C Library.
   5    Contributed by Jakub Jelinek (jj@ultra.linux.cz).
   6
   7    The GNU C Library is free software; you can redistribute it and/or
   8    modify it under the terms of the GNU Library General Public License as
   9    published by the Free Software Foundation; either version 2 of the
  10    License, or (at your option) any later version.
  11
  12    The GNU C Library is distributed in the hope that it will be useful,
  13    but WITHOUT ANY WARRANTY; without even the implied warranty of
  14    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15    Library General Public License for more details.
  16
  17    You should have received a copy of the GNU Library General Public
  18    License along with the GNU C Library; see the file COPYING.LIB.  If
  19    not, write to the Free Software Foundation, Inc.,
  20    59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  */
  21
  22
  23 #ifndef    __MATH_EMU_EXTENDED_H__
  24 #define    __MATH_EMU_EXTENDED_H__
  25
  26 #if _FP_W_TYPE_SIZE < 32
  27 #error "Here's a nickel, kid. Go buy yourself a real computer."
  28 #endif
  29
  30 #if _FP_W_TYPE_SIZE < 64
  31 #define _FP_FRACTBITS_E         (4*_FP_W_TYPE_SIZE)
  32 #else
  33 #define _FP_FRACTBITS_E         (2*_FP_W_TYPE_SIZE)
  34 #endif
  35
  36 #define _FP_FRACBITS_E          64
  37 #define _FP_FRACXBITS_E         (_FP_FRACTBITS_E - _FP_FRACBITS_E)
  38 #define _FP_WFRACBITS_E         (_FP_WORKBITS + _FP_FRACBITS_E)
  39 #define _FP_WFRACXBITS_E        (_FP_FRACTBITS_E - _FP_WFRACBITS_E)
  40 #define _FP_EXPBITS_E           15
  41 #define _FP_EXPBIAS_E           16383
  42 #define _FP_EXPMAX_E            32767
  43
  44 #define _FP_QNANBIT_E           \
  45         ((_FP_W_TYPE)1 << (_FP_FRACBITS_E-2) % _FP_W_TYPE_SIZE)
  46 #define _FP_IMPLBIT_E           \
  47         ((_FP_W_TYPE)1 << (_FP_FRACBITS_E-1) % _FP_W_TYPE_SIZE)
  48 #define _FP_OVERFLOW_E          \
  49         ((_FP_W_TYPE)1 << (_FP_WFRACBITS_E % _FP_W_TYPE_SIZE))
  50
  51 #if _FP_W_TYPE_SIZE < 64
  52
  53 union _FP_UNION_E
  54 {
  55    long double flt;
  56    struct
  57    {
  58 #if __BYTE_ORDER == __BIG_ENDIAN
  59       unsigned long pad1 : _FP_W_TYPE_SIZE;
  60       unsigned long pad2 : (_FP_W_TYPE_SIZE - 1 - _FP_EXPBITS_E);
  61       unsigned long sign : 1;
  62       unsigned long exp : _FP_EXPBITS_E;
  63       unsigned long frac1 : _FP_W_TYPE_SIZE;
  64       unsigned long frac0 : _FP_W_TYPE_SIZE;
  65 #else
  66       unsigned long frac0 : _FP_W_TYPE_SIZE;
  67       unsigned long frac1 : _FP_W_TYPE_SIZE;
  68       unsigned exp : _FP_EXPBITS_E;
  69       unsigned sign : 1;
  70 #endif /* not bigendian */
  71    } bits __attribute__((packed));
  72 };
  73
  74
  75 #define FP_DECL_E(X)            _FP_DECL(4,X)
  76
  77 #define FP_UNPACK_RAW_E(X, val)                         \
  78   do {                                                  \
  79     union _FP_UNION_E _flo; _flo.flt = (val);           \
  80                                                         \
  81     X##_f[2] = 0; X##_f[3] = 0;                         \
  82     X##_f[0] = _flo.bits.frac0;                         \
  83     X##_f[1] = _flo.bits.frac1;                         \
  84     X##_e  = _flo.bits.exp;                             \
  85     X##_s  = _flo.bits.sign;                            \
  86     if (!X##_e && (X##_f[1] || X##_f[0])                \
  87         && !(X##_f[1] & _FP_IMPLBIT_E))                 \
  88       {                                                 \
  89         X##_e++;                                        \
  90         FP_SET_EXCEPTION(FP_EX_DENORM);                 \
  91       }                                                 \
  92   } while (0)
  93
  94 #define FP_UNPACK_RAW_EP(X, val)                        \
  95   do {                                                  \
  96     union _FP_UNION_E *_flo =                           \
  97     (union _FP_UNION_E *)(val);                         \
  98                                                         \
  99     X##_f[2] = 0; X##_f[3] = 0;                         \
 100     X##_f[0] = _flo->bits.frac0;                        \
 101     X##_f[1] = _flo->bits.frac1;                        \
 102     X##_e  = _flo->bits.exp;                            \
 103     X##_s  = _flo->bits.sign;                           \
 104     if (!X##_e && (X##_f[1] || X##_f[0])                \
 105         && !(X##_f[1] & _FP_IMPLBIT_E))                 \
 106       {                                                 \
 107         X##_e++;                                        \
 108         FP_SET_EXCEPTION(FP_EX_DENORM);                 \
 109       }                                                 \
 110   } while (0)
 111
 112 #define FP_PACK_RAW_E(val, X)                           \
 113   do {                                                  \
 114     union _FP_UNION_E _flo;                             \
 115                                                         \
 116     if (X##_e) X##_f[1] |= _FP_IMPLBIT_E;               \
 117     else X##_f[1] &= ~(_FP_IMPLBIT_E);                  \
 118     _flo.bits.frac0 = X##_f[0];                         \
 119     _flo.bits.frac1 = X##_f[1];                         \
 120     _flo.bits.exp   = X##_e;                            \
 121     _flo.bits.sign  = X##_s;                            \
 122                                                         \
 123     (val) = _flo.flt;                                   \
 124   } while (0)
 125
 126 #define FP_PACK_RAW_EP(val, X)                          \
 127   do {                                                  \
 128     if (!FP_INHIBIT_RESULTS)                            \
 129       {                                                 \
 130         union _FP_UNION_E *_flo =                       \
 131           (union _FP_UNION_E *)(val);                   \
 132                                                         \
 133         if (X##_e) X##_f[1] |= _FP_IMPLBIT_E;           \
 134         else X##_f[1] &= ~(_FP_IMPLBIT_E);              \
 135         _flo->bits.frac0 = X##_f[0];                    \
 136         _flo->bits.frac1 = X##_f[1];                    \
 137         _flo->bits.exp   = X##_e;                       \
 138         _flo->bits.sign  = X##_s;                       \
 139       }                                                 \
 140   } while (0)
 141
 142 #define FP_UNPACK_E(X,val)              \
 143   do {                                  \
 144     FP_UNPACK_RAW_E(X,val);             \
 145     _FP_UNPACK_CANONICAL(E,4,X);        \
 146   } while (0)
 147
 148 #define FP_UNPACK_EP(X,val)             \
 149   do {                                  \
 150     FP_UNPACK_RAW_2_P(X,val);           \
 151     _FP_UNPACK_CANONICAL(E,4,X);        \
 152   } while (0)
 153
 154 #define FP_PACK_E(val,X)                \
 155   do {                                  \
 156     _FP_PACK_CANONICAL(E,4,X);          \
 157     FP_PACK_RAW_E(val,X);               \
 158   } while (0)
 159
 160 #define FP_PACK_EP(val,X)               \
 161   do {                                  \
 162     _FP_PACK_CANONICAL(E,4,X);          \
 163     FP_PACK_RAW_EP(val,X);              \
 164   } while (0)
 165
 166 #define FP_ISSIGNAN_E(X)        _FP_ISSIGNAN(E,4,X)
 167 #define FP_NEG_E(R,X)           _FP_NEG(E,4,R,X)
 168 #define FP_ADD_E(R,X,Y)         _FP_ADD(E,4,R,X,Y)
 169 #define FP_SUB_E(R,X,Y)         _FP_SUB(E,4,R,X,Y)
 170 #define FP_MUL_E(R,X,Y)         _FP_MUL(E,4,R,X,Y)
 171 #define FP_DIV_E(R,X,Y)         _FP_DIV(E,4,R,X,Y)
 172 #define FP_SQRT_E(R,X)          _FP_SQRT(E,4,R,X)
 173
 174 /*
 175  * Square root algorithms:
 176  * We have just one right now, maybe Newton approximation
 177  * should be added for those machines where division is fast.
 178  * This has special _E version because standard _4 square
 179  * root would not work (it has to start normally with the
 180  * second word and not the first), but as we have to do it
 181  * anyway, we optimize it by doing most of the calculations
 182  * in two UWtype registers instead of four.
 183  */
 184
 185 #define _FP_SQRT_MEAT_E(R, S, T, X, q)                  \
 186   do {                                                  \
 187     q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1);         \
 188     _FP_FRAC_SRL_4(X, (_FP_WORKBITS));                  \
 189     while (q)                                           \
 190       {                                                 \
 191         T##_f[1] = S##_f[1] + q;                        \
 192         if (T##_f[1] <= X##_f[1])                       \
 193           {                                             \
 194             S##_f[1] = T##_f[1] + q;                    \
 195             X##_f[1] -= T##_f[1];                       \
 196             R##_f[1] += q;                              \
 197           }                                             \
 198         _FP_FRAC_SLL_2(X, 1);                           \
 199         q >>= 1;                                        \
 200       }                                                 \
 201     q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1);         \
 202     while (q)                                           \
 203       {                                                 \
 204         T##_f[0] = S##_f[0] + q;                        \
 205         T##_f[1] = S##_f[1];                            \
 206         if (T##_f[1] < X##_f[1] ||                      \
 207             (T##_f[1] == X##_f[1] &&                    \
 208              T##_f[0] <= X##_f[0]))                     \
 209           {                                             \
 210             S##_f[0] = T##_f[0] + q;                    \
 211             S##_f[1] += (T##_f[0] > S##_f[0]);          \
 212             _FP_FRAC_DEC_2(X, T);                       \
 213             R##_f[0] += q;                              \
 214           }                                             \
 215         _FP_FRAC_SLL_2(X, 1);                           \
 216         q >>= 1;                                        \
 217       }                                                 \
 218     _FP_FRAC_SLL_4(R, (_FP_WORKBITS));                  \
 219     if (X##_f[0] | X##_f[1])                            \
 220       {                                                 \
 221         if (S##_f[1] < X##_f[1] ||                      \
 222             (S##_f[1] == X##_f[1] &&                    \
 223              S##_f[0] < X##_f[0]))                      \
 224           R##_f[0] |= _FP_WORK_ROUND;                   \
 225         R##_f[0] |= _FP_WORK_STICKY;                    \
 226       }                                                 \
 227   } while (0)
 228
 229 #define FP_CMP_E(r,X,Y,un)      _FP_CMP(E,4,r,X,Y,un)
 230 #define FP_CMP_EQ_E(r,X,Y)      _FP_CMP_EQ(E,4,r,X,Y)
 231
 232 #define FP_TO_INT_E(r,X,rsz,rsg)        _FP_TO_INT(E,4,r,X,rsz,rsg)
 233 #define FP_TO_INT_ROUND_E(r,X,rsz,rsg)  _FP_TO_INT_ROUND(E,4,r,X,rsz,rsg)
 234 #define FP_FROM_INT_E(X,r,rs,rt)        _FP_FROM_INT(E,4,X,r,rs,rt)
 235
 236 #define _FP_FRAC_HIGH_E(X)      (X##_f[2])
 237 #define _FP_FRAC_HIGH_RAW_E(X)  (X##_f[1])
 238
 239 #else   /* not _FP_W_TYPE_SIZE < 64 */
 240 union _FP_UNION_E
 241 {
 242   long double flt /* __attribute__((mode(TF))) */ ;
 243   struct {
 244 #if __BYTE_ORDER == __BIG_ENDIAN
 245     unsigned long pad : (_FP_W_TYPE_SIZE - 1 - _FP_EXPBITS_E);
 246     unsigned sign  : 1;
 247     unsigned exp   : _FP_EXPBITS_E;
 248     unsigned long frac : _FP_W_TYPE_SIZE;
 249 #else
 250     unsigned long frac : _FP_W_TYPE_SIZE;
 251     unsigned exp   : _FP_EXPBITS_E;
 252     unsigned sign  : 1;
 253 #endif
 254   } bits;
 255 };
 256
 257 #define FP_DECL_E(X)            _FP_DECL(2,X)
 258
 259 #define FP_UNPACK_RAW_E(X, val)                                 \
 260   do {                                                          \
 261     union _FP_UNION_E _flo; _flo.flt = (val);                   \
 262                                                                 \
 263     X##_f0 = _flo.bits.frac;                                    \
 264     X##_f1 = 0;                                                 \
 265     X##_e = _flo.bits.exp;                                      \
 266     X##_s = _flo.bits.sign;                                     \
 267     if (!X##_e && X##_f0 && !(X##_f0 & _FP_IMPLBIT_E))          \
 268       {                                                         \
 269         X##_e++;                                                \
 270         FP_SET_EXCEPTION(FP_EX_DENORM);                         \
 271       }                                                         \
 272   } while (0)
 273
 274 #define FP_UNPACK_RAW_EP(X, val)                                \
 275   do {                                                          \
 276     union _FP_UNION_E *_flo =                                   \
 277       (union _FP_UNION_E *)(val);                               \
 278                                                                 \
 279     X##_f0 = _flo->bits.frac;                                   \
 280     X##_f1 = 0;                                                 \
 281     X##_e = _flo->bits.exp;                                     \
 282     X##_s = _flo->bits.sign;                                    \
 283     if (!X##_e && X##_f0 && !(X##_f0 & _FP_IMPLBIT_E))          \
 284       {                                                         \
 285         X##_e++;                                                \
 286         FP_SET_EXCEPTION(FP_EX_DENORM);                         \
 287       }                                                         \
 288   } while (0)
 289
 290 #define FP_PACK_RAW_E(val, X)                                   \
 291   do {                                                          \
 292     union _FP_UNION_E _flo;                                     \
 293                                                                 \
 294     if (X##_e) X##_f0 |= _FP_IMPLBIT_E;                         \
 295     else X##_f0 &= ~(_FP_IMPLBIT_E);                            \
 296     _flo.bits.frac = X##_f0;                                    \
 297     _flo.bits.exp  = X##_e;                                     \
 298     _flo.bits.sign = X##_s;                                     \
 299                                                                 \
 300     (val) = _flo.flt;                                           \
 301   } while (0)
 302
 303 #define FP_PACK_RAW_EP(fs, val, X)                              \
 304   do {                                                          \
 305     if (!FP_INHIBIT_RESULTS)                                    \
 306       {                                                         \
 307         union _FP_UNION_E *_flo =                               \
 308           (union _FP_UNION_E *)(val);                           \
 309                                                                 \
 310         if (X##_e) X##_f0 |= _FP_IMPLBIT_E;                     \
 311         else X##_f0 &= ~(_FP_IMPLBIT_E);                        \
 312         _flo->bits.frac = X##_f0;                               \
 313         _flo->bits.exp  = X##_e;                                \
 314         _flo->bits.sign = X##_s;                                \
 315       }                                                         \
 316   } while (0)
 317
 318
 319 #define FP_UNPACK_E(X,val)              \
 320   do {                                  \
 321     FP_UNPACK_RAW_E(X,val);             \
 322     _FP_UNPACK_CANONICAL(E,2,X);        \
 323   } while (0)
 324
 325 #define FP_UNPACK_EP(X,val)             \
 326   do {                                  \
 327     FP_UNPACK_RAW_EP(X,val);            \
 328     _FP_UNPACK_CANONICAL(E,2,X);        \
 329   } while (0)
 330
 331 #define FP_PACK_E(val,X)                \
 332   do {                                  \
 333     _FP_PACK_CANONICAL(E,2,X);          \
 334     FP_PACK_RAW_E(val,X);               \
 335   } while (0)
 336
 337 #define FP_PACK_EP(val,X)               \
 338   do {                                  \
 339     _FP_PACK_CANONICAL(E,2,X);          \
 340     FP_PACK_RAW_EP(val,X);              \
 341   } while (0)
 342
 343 #define FP_ISSIGNAN_E(X)        _FP_ISSIGNAN(E,2,X)
 344 #define FP_NEG_E(R,X)           _FP_NEG(E,2,R,X)
 345 #define FP_ADD_E(R,X,Y)         _FP_ADD(E,2,R,X,Y)
 346 #define FP_SUB_E(R,X,Y)         _FP_SUB(E,2,R,X,Y)
 347 #define FP_MUL_E(R,X,Y)         _FP_MUL(E,2,R,X,Y)
 348 #define FP_DIV_E(R,X,Y)         _FP_DIV(E,2,R,X,Y)
 349 #define FP_SQRT_E(R,X)          _FP_SQRT(E,2,R,X)
 350
 351 /*
 352  * Square root algorithms:
 353  * We have just one right now, maybe Newton approximation
 354  * should be added for those machines where division is fast.
 355  * We optimize it by doing most of the calculations
 356  * in one UWtype registers instead of two, although we don't
 357  * have to.
 358  */
 359 #define _FP_SQRT_MEAT_E(R, S, T, X, q)                  \
 360   do {                                                  \
 361     q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1);         \
 362     _FP_FRAC_SRL_2(X, (_FP_WORKBITS));                  \
 363     while (q)                                           \
 364       {                                                 \
 365         T##_f0 = S##_f0 + q;                            \
 366         if (T##_f0 <= X##_f0)                           \
 367           {                                             \
 368             S##_f0 = T##_f0 + q;                        \
 369             X##_f0 -= T##_f0;                           \
 370             R##_f0 += q;                                \
 371           }                                             \
 372         _FP_FRAC_SLL_1(X, 1);                           \
 373         q >>= 1;                                        \
 374       }                                                 \
 375     _FP_FRAC_SLL_2(R, (_FP_WORKBITS));                  \
 376     if (X##_f0)                                         \
 377       {                                                 \
 378         if (S##_f0 < X##_f0)                            \
 379           R##_f0 |= _FP_WORK_ROUND;                     \
 380         R##_f0 |= _FP_WORK_STICKY;                      \
 381       }                                                 \
 382   } while (0)
 383
 384 #define FP_CMP_E(r,X,Y,un)      _FP_CMP(E,2,r,X,Y,un)
 385 #define FP_CMP_EQ_E(r,X,Y)      _FP_CMP_EQ(E,2,r,X,Y)
 386
 387 #define FP_TO_INT_E(r,X,rsz,rsg)        _FP_TO_INT(E,2,r,X,rsz,rsg)
 388 #define FP_TO_INT_ROUND_E(r,X,rsz,rsg)  _FP_TO_INT_ROUND(E,2,r,X,rsz,rsg)
 389 #define FP_FROM_INT_E(X,r,rs,rt)        _FP_FROM_INT(E,2,X,r,rs,rt)
 390
 391 #define _FP_FRAC_HIGH_E(X)      (X##_f1)
 392 #define _FP_FRAC_HIGH_RAW_E(X)  (X##_f0)
 393
 394 #endif /* not _FP_W_TYPE_SIZE < 64 */
 395
 396 #endif /* __MATH_EMU_EXTENDED_H__ */