soft-fp/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 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_IMPLBIT_E           \
  43         ((_FP_W_TYPE)1 << (_FP_FRACBITS_E-1) % _FP_W_TYPE_SIZE)
  44 #define _FP_OVERFLOW_E          \
  45         ((_FP_W_TYPE)1 << (_FP_WFRACBITS_E % _FP_W_TYPE_SIZE))
  46
  47 #if _FP_W_TYPE_SIZE < 64
  48
  49 union _FP_UNION_E
  50 {
  51    long double flt;
  52    struct
  53    {
  54 #if __BYTE_ORDER == __BIG_ENDIAN
  55       unsigned long pad1 : _FP_W_TYPE_SIZE;
  56       unsigned long pad2 : (_FP_W_TYPE_SIZE - 1 - _FP_EXPBITS_E);
  57       unsigned long sign : 1;
  58       unsigned long exp : _FP_EXPBITS_E;
  59       unsigned long frac1 : _FP_W_TYPE_SIZE;
  60       unsigned long frac0 : _FP_W_TYPE_SIZE;
  61 #else
  62       unsigned long frac0 : _FP_W_TYPE_SIZE;
  63       unsigned long frac1 : _FP_W_TYPE_SIZE;
  64       unsigned exp : _FP_EXPBITS_E;
  65       unsigned sign : 1;
  66 #endif /* not bigendian */
  67    } bits __attribute__((packed));
  68 };
  69
  70
  71 #define FP_DECL_E(X)            _FP_DECL(4,X)
  72
  73 #define FP_UNPACK_RAW_E(X, val)                         \
  74   do {                                                  \
  75     union _FP_UNION_E _flo; _flo.flt = (val);           \
  76                                                         \
  77     X##_f[2] = 0; X##_f[3] = 0;                         \
  78     X##_f[0] = _flo.bits.frac0;                         \
  79     X##_f[1] = _flo.bits.frac1;                         \
  80     X##_e  = _flo.bits.exp;                             \
  81     X##_s  = _flo.bits.sign;                            \
  82     if (!X##_e && (X##_f[1] || X##_f[0])                \
  83         && !(X##_f[1] & _FP_IMPLBIT_E))                 \
  84       {                                                 \
  85         X##_e++;                                        \
  86         FP_SET_EXCEPTION(FP_EX_DENORM);                 \
  87       }                                                 \
  88   } while (0)
  89
  90 #define FP_UNPACK_RAW_EP(X, val)                        \
  91   do {                                                  \
  92     union _FP_UNION_E *_flo =                           \
  93     (union _FP_UNION_E *)(val);                         \
  94                                                         \
  95     X##_f[2] = 0; X##_f[3] = 0;                         \
  96     X##_f[0] = _flo->bits.frac0;                        \
  97     X##_f[1] = _flo->bits.frac1;                        \
  98     X##_e  = _flo->bits.exp;                            \
  99     X##_s  = _flo->bits.sign;                           \
 100     if (!X##_e && (X##_f[1] || X##_f[0])                \
 101         && !(X##_f[1] & _FP_IMPLBIT_E))                 \
 102       {                                                 \
 103         X##_e++;                                        \
 104         FP_SET_EXCEPTION(FP_EX_DENORM);                 \
 105       }                                                 \
 106   } while (0)
 107
 108 #define FP_PACK_RAW_E(val, X)                           \
 109   do {                                                  \
 110     union _FP_UNION_E _flo;                             \
 111                                                         \
 112     if (X##_e) X##_f[1] |= _FP_IMPLBIT_E;               \
 113     else X##_f[1] &= ~(_FP_IMPLBIT_E);                  \
 114     _flo.bits.frac0 = X##_f[0];                         \
 115     _flo.bits.frac1 = X##_f[1];                         \
 116     _flo.bits.exp   = X##_e;                            \
 117     _flo.bits.sign  = X##_s;                            \
 118                                                         \
 119     (val) = _flo.flt;                                   \
 120   } while (0)
 121
 122 #define FP_PACK_RAW_EP(val, X)                          \
 123   do {                                                  \
 124     if (!FP_INHIBIT_RESULTS)                            \
 125       {                                                 \
 126         union _FP_UNION_E *_flo =                       \
 127           (union _FP_UNION_E *)(val);                   \
 128                                                         \
 129         if (X##_e) X##_f[1] |= _FP_IMPLBIT_E;           \
 130         else X##_f[1] &= ~(_FP_IMPLBIT_E);              \
 131         _flo->bits.frac0 = X##_f[0];                    \
 132         _flo->bits.frac1 = X##_f[1];                    \
 133         _flo->bits.exp   = X##_e;                       \
 134         _flo->bits.sign  = X##_s;                       \
 135       }                                                 \
 136   } while (0)
 137
 138 #define FP_UNPACK_E(X,val)              \
 139   do {                                  \
 140     FP_UNPACK_RAW_E(X,val);             \
 141     _FP_UNPACK_CANONICAL(E,4,X);        \
 142   } while (0)
 143
 144 #define FP_UNPACK_EP(X,val)             \
 145   do {                                  \
 146     FP_UNPACK_RAW_2_P(X,val);           \
 147     _FP_UNPACK_CANONICAL(E,4,X);        \
 148   } while (0)
 149
 150 #define FP_PACK_E(val,X)                \
 151   do {                                  \
 152     _FP_PACK_CANONICAL(E,4,X);          \
 153     FP_PACK_RAW_E(val,X);               \
 154   } while (0)
 155
 156 #define FP_PACK_EP(val,X)               \
 157   do {                                  \
 158     _FP_PACK_CANONICAL(E,4,X);          \
 159     FP_PACK_RAW_EP(val,X);              \
 160   } while (0)
 161
 162 #define FP_ISSIGNAN_E(X)        _FP_ISSIGNAN(E,4,X)
 163 #define FP_NEG_E(R,X)           _FP_NEG(E,4,R,X)
 164 #define FP_ADD_E(R,X,Y)         _FP_ADD(E,4,R,X,Y)
 165 #define FP_SUB_E(R,X,Y)         _FP_SUB(E,4,R,X,Y)
 166 #define FP_MUL_E(R,X,Y)         _FP_MUL(E,4,R,X,Y)
 167 #define FP_DIV_E(R,X,Y)         _FP_DIV(E,4,R,X,Y)
 168 #define FP_SQRT_E(R,X)          _FP_SQRT(E,4,R,X)
 169
 170 /*
 171  * Square root algorithms:
 172  * We have just one right now, maybe Newton approximation
 173  * should be added for those machines where division is fast.
 174  * This has special _E version because standard _4 square
 175  * root would not work (it has to start normally with the
 176  * second word and not the first), but as we have to do it
 177  * anyway, we optimize it by doing most of the calculations
 178  * in two UWtype registers instead of four.
 179  */
 180
 181 #define _FP_SQRT_MEAT_E(R, S, T, X, q)                  \
 182   do {                                                  \
 183     q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1);         \
 184     _FP_FRAC_SRL_4(X, (_FP_WORKBITS));                  \
 185     while (q)                                           \
 186       {                                                 \
 187         T##_f[1] = S##_f[1] + q;                        \
 188         if (T##_f[1] <= X##_f[1])                       \
 189           {                                             \
 190             S##_f[1] = T##_f[1] + q;                    \
 191             X##_f[1] -= T##_f[1];                       \
 192             R##_f[1] += q;                              \
 193           }                                             \
 194         _FP_FRAC_SLL_2(X, 1);                           \
 195         q >>= 1;                                        \
 196       }                                                 \
 197     q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1);         \
 198     while (q)                                           \
 199       {                                                 \
 200         T##_f[0] = S##_f[0] + q;                        \
 201         T##_f[1] = S##_f[1];                            \
 202         if (T##_f[1] < X##_f[1] ||                      \
 203             (T##_f[1] == X##_f[1] &&                    \
 204              T##_f[0] <= X##_f[0]))                     \
 205           {                                             \
 206             S##_f[0] = T##_f[0] + q;                    \
 207             S##_f[1] += (T##_f[0] > S##_f[0]);          \
 208             _FP_FRAC_DEC_2(X, T);                       \
 209             R##_f[0] += q;                              \
 210           }                                             \
 211         _FP_FRAC_SLL_2(X, 1);                           \
 212         q >>= 1;                                        \
 213       }                                                 \
 214     _FP_FRAC_SLL_4(R, (_FP_WORKBITS));                  \
 215     if (X##_f[0] | X##_f[1])                            \
 216       {                                                 \
 217         if (S##_f[1] < X##_f[1] ||                      \
 218             (S##_f[1] == X##_f[1] &&                    \
 219              S##_f[0] < X##_f[0]))                      \
 220           R##_f[0] |= _FP_WORK_ROUND;                   \
 221         R##_f[0] |= _FP_WORK_STICKY;                    \
 222       }                                                 \
 223   } while (0)
 224
 225 #define FP_CMP_E(r,X,Y,un)      _FP_CMP(E,4,r,X,Y,un)
 226 #define FP_CMP_EQ_E(r,X,Y)      _FP_CMP_EQ(E,4,r,X,Y)
 227
 228 #define FP_TO_INT_E(r,X,rsz,rsg)        _FP_TO_INT(E,4,r,X,rsz,rsg)
 229 #define FP_FROM_INT_E(X,r,rs,rt)        _FP_FROM_INT(E,4,X,r,rs,rt)
 230
 231 #define _FP_FRAC_HIGH_E(X)      (X##_f[2])
 232 #define _FP_FRAC_HIGH_RAW_E(X)  (X##_f[1])
 233
 234 #else   /* not _FP_W_TYPE_SIZE < 64 */
 235 union _FP_UNION_E
 236 {
 237   long double flt /* __attribute__((mode(TF))) */ ;
 238   struct {
 239 #if __BYTE_ORDER == __BIG_ENDIAN
 240     unsigned long pad : (_FP_W_TYPE_SIZE - 1 - _FP_EXPBITS_E);
 241     unsigned sign  : 1;
 242     unsigned exp   : _FP_EXPBITS_E;
 243     unsigned long frac : _FP_W_TYPE_SIZE;
 244 #else
 245     unsigned long frac : _FP_W_TYPE_SIZE;
 246     unsigned exp   : _FP_EXPBITS_E;
 247     unsigned sign  : 1;
 248 #endif
 249   } bits;
 250 };
 251
 252 #define FP_DECL_E(X)            _FP_DECL(2,X)
 253
 254 #define FP_UNPACK_RAW_E(X, val)                                 \
 255   do {                                                          \
 256     union _FP_UNION_E _flo; _flo.flt = (val);                   \
 257                                                                 \
 258     X##_f0 = _flo.bits.frac;                                    \
 259     X##_f1 = 0;                                                 \
 260     X##_e = _flo.bits.exp;                                      \
 261     X##_s = _flo.bits.sign;                                     \
 262     if (!X##_e && X##_f0 && !(X##_f0 & _FP_IMPLBIT_E))          \
 263       {                                                         \
 264         X##_e++;                                                \
 265         FP_SET_EXCEPTION(FP_EX_DENORM);                         \
 266       }                                                         \
 267   } while (0)
 268
 269 #define FP_UNPACK_RAW_EP(X, val)                                \
 270   do {                                                          \
 271     union _FP_UNION_E *_flo =                                   \
 272       (union _FP_UNION_E *)(val);                               \
 273                                                                 \
 274     X##_f0 = _flo->bits.frac;                                   \
 275     X##_f1 = 0;                                                 \
 276     X##_e = _flo->bits.exp;                                     \
 277     X##_s = _flo->bits.sign;                                    \
 278     if (!X##_e && X##_f0 && !(X##_f0 & _FP_IMPLBIT_E))          \
 279       {                                                         \
 280         X##_e++;                                                \
 281         FP_SET_EXCEPTION(FP_EX_DENORM);                         \
 282       }                                                         \
 283   } while (0)
 284
 285 #define FP_PACK_RAW_E(val, X)                                   \
 286   do {                                                          \
 287     union _FP_UNION_E _flo;                                     \
 288                                                                 \
 289     if (X##_e) X##_f0 |= _FP_IMPLBIT_E;                         \
 290     else X##_f0 &= ~(_FP_IMPLBIT_E);                            \
 291     _flo.bits.frac = X##_f0;                                    \
 292     _flo.bits.exp  = X##_e;                                     \
 293     _flo.bits.sign = X##_s;                                     \
 294                                                                 \
 295     (val) = _flo.flt;                                           \
 296   } while (0)
 297
 298 #define FP_PACK_RAW_EP(fs, val, X)                              \
 299   do {                                                          \
 300     if (!FP_INHIBIT_RESULTS)                                    \
 301       {                                                         \
 302         union _FP_UNION_E *_flo =                               \
 303           (union _FP_UNION_E *)(val);                           \
 304                                                                 \
 305         if (X##_e) X##_f0 |= _FP_IMPLBIT_E;                     \
 306         else X##_f0 &= ~(_FP_IMPLBIT_E);                        \
 307         _flo->bits.frac = X##_f0;                               \
 308         _flo->bits.exp  = X##_e;                                \
 309         _flo->bits.sign = X##_s;                                \
 310       }                                                         \
 311   } while (0)
 312
 313
 314 #define FP_UNPACK_E(X,val)              \
 315   do {                                  \
 316     FP_UNPACK_RAW_E(X,val);             \
 317     _FP_UNPACK_CANONICAL(E,2,X);        \
 318   } while (0)
 319
 320 #define FP_UNPACK_EP(X,val)             \
 321   do {                                  \
 322     FP_UNPACK_RAW_EP(X,val);            \
 323     _FP_UNPACK_CANONICAL(E,2,X);        \
 324   } while (0)
 325
 326 #define FP_PACK_E(val,X)                \
 327   do {                                  \
 328     _FP_PACK_CANONICAL(E,2,X);          \
 329     FP_PACK_RAW_E(val,X);               \
 330   } while (0)
 331
 332 #define FP_PACK_EP(val,X)               \
 333   do {                                  \
 334     _FP_PACK_CANONICAL(E,2,X);          \
 335     FP_PACK_RAW_EP(val,X);              \
 336   } while (0)
 337
 338 #define FP_ISSIGNAN_E(X)        _FP_ISSIGNAN(E,2,X)
 339 #define FP_NEG_E(R,X)           _FP_NEG(E,2,R,X)
 340 #define FP_ADD_E(R,X,Y)         _FP_ADD(E,2,R,X,Y)
 341 #define FP_SUB_E(R,X,Y)         _FP_SUB(E,2,R,X,Y)
 342 #define FP_MUL_E(R,X,Y)         _FP_MUL(E,2,R,X,Y)
 343 #define FP_DIV_E(R,X,Y)         _FP_DIV(E,2,R,X,Y)
 344 #define FP_SQRT_E(R,X)          _FP_SQRT(E,2,R,X)
 345
 346 /*
 347  * Square root algorithms:
 348  * We have just one right now, maybe Newton approximation
 349  * should be added for those machines where division is fast.
 350  * We optimize it by doing most of the calculations
 351  * in one UWtype registers instead of two, although we don't
 352  * have to.
 353  */
 354 #define _FP_SQRT_MEAT_E(R, S, T, X, q)                  \
 355   do {                                                  \
 356     q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1);         \
 357     _FP_FRAC_SRL_2(X, (_FP_WORKBITS));                  \
 358     while (q)                                           \
 359       {                                                 \
 360         T##_f0 = S##_f0 + q;                            \
 361         if (T##_f0 <= X##_f0)                           \
 362           {                                             \
 363             S##_f0 = T##_f0 + q;                        \
 364             X##_f0 -= T##_f0;                           \
 365             R##_f0 += q;                                \
 366           }                                             \
 367         _FP_FRAC_SLL_1(X, 1);                           \
 368         q >>= 1;                                        \
 369       }                                                 \
 370     _FP_FRAC_SLL_2(R, (_FP_WORKBITS));                  \
 371     if (X##_f0)                                         \
 372       {                                                 \
 373         if (S##_f0 < X##_f0)                            \
 374           R##_f0 |= _FP_WORK_ROUND;                     \
 375         R##_f0 |= _FP_WORK_STICKY;                      \
 376       }                                                 \
 377   } while (0)
 378
 379 #define FP_CMP_E(r,X,Y,un)      _FP_CMP(E,2,r,X,Y,un)
 380 #define FP_CMP_EQ_E(r,X,Y)      _FP_CMP_EQ(E,2,r,X,Y)
 381
 382 #define FP_TO_INT_E(r,X,rsz,rsg)        _FP_TO_INT(E,2,r,X,rsz,rsg)
 383 #define FP_FROM_INT_E(X,r,rs,rt)        _FP_FROM_INT(E,2,X,r,rs,rt)
 384
 385 #define _FP_FRAC_HIGH_E(X)      (X##_f1)
 386 #define _FP_FRAC_HIGH_RAW_E(X)  (X##_f0)
 387
 388 #endif /* not _FP_W_TYPE_SIZE < 64 */