soft-fp/extended.h

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