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 #define FP_CMP_UNORD_E(r,X,Y)   _FP_CMP_UNORD(E,4,r,X,Y)
 228
 229 #define FP_TO_INT_E(r,X,rsz,rsg)        _FP_TO_INT(E,4,r,X,rsz,rsg)
 230 #define FP_FROM_INT_E(X,r,rs,rt)        _FP_FROM_INT(E,4,X,r,rs,rt)
 231
 232 #define _FP_FRAC_HIGH_E(X)      (X##_f[2])
 233 #define _FP_FRAC_HIGH_RAW_E(X)  (X##_f[1])
 234
 235 #else   /* not _FP_W_TYPE_SIZE < 64 */
 236 union _FP_UNION_E
 237 {
 238   long double flt /* __attribute__((mode(TF))) */ ;
 239   struct {
 240 #if __BYTE_ORDER == __BIG_ENDIAN
 241     unsigned long pad : (_FP_W_TYPE_SIZE - 1 - _FP_EXPBITS_E);
 242     unsigned sign  : 1;
 243     unsigned exp   : _FP_EXPBITS_E;
 244     unsigned long frac : _FP_W_TYPE_SIZE;
 245 #else
 246     unsigned long frac : _FP_W_TYPE_SIZE;
 247     unsigned exp   : _FP_EXPBITS_E;
 248     unsigned sign  : 1;
 249 #endif
 250   } bits;
 251 };
 252
 253 #define FP_DECL_E(X)            _FP_DECL(2,X)
 254
 255 #define FP_UNPACK_RAW_E(X, val)                                 \
 256   do {                                                          \
 257     union _FP_UNION_E _flo; _flo.flt = (val);                   \
 258                                                                 \
 259     X##_f0 = _flo.bits.frac;                                    \
 260     X##_f1 = 0;                                                 \
 261     X##_e = _flo.bits.exp;                                      \
 262     X##_s = _flo.bits.sign;                                     \
 263     if (!X##_e && X##_f0 && !(X##_f0 & _FP_IMPLBIT_E))          \
 264       {                                                         \
 265         X##_e++;                                                \
 266         FP_SET_EXCEPTION(FP_EX_DENORM);                         \
 267       }                                                         \
 268   } while (0)
 269
 270 #define FP_UNPACK_RAW_EP(X, val)                                \
 271   do {                                                          \
 272     union _FP_UNION_E *_flo =                                   \
 273       (union _FP_UNION_E *)(val);                               \
 274                                                                 \
 275     X##_f0 = _flo->bits.frac;                                   \
 276     X##_f1 = 0;                                                 \
 277     X##_e = _flo->bits.exp;                                     \
 278     X##_s = _flo->bits.sign;                                    \
 279     if (!X##_e && X##_f0 && !(X##_f0 & _FP_IMPLBIT_E))          \
 280       {                                                         \
 281         X##_e++;                                                \
 282         FP_SET_EXCEPTION(FP_EX_DENORM);                         \
 283       }                                                         \
 284   } while (0)
 285
 286 #define FP_PACK_RAW_E(val, X)                                   \
 287   do {                                                          \
 288     union _FP_UNION_E _flo;                                     \
 289                                                                 \
 290     if (X##_e) X##_f0 |= _FP_IMPLBIT_E;                         \
 291     else X##_f0 &= ~(_FP_IMPLBIT_E);                            \
 292     _flo.bits.frac = X##_f0;                                    \
 293     _flo.bits.exp  = X##_e;                                     \
 294     _flo.bits.sign = X##_s;                                     \
 295                                                                 \
 296     (val) = _flo.flt;                                           \
 297   } while (0)
 298
 299 #define FP_PACK_RAW_EP(fs, val, X)                              \
 300   do {                                                          \
 301     if (!FP_INHIBIT_RESULTS)                                    \
 302       {                                                         \
 303         union _FP_UNION_E *_flo =                               \
 304           (union _FP_UNION_E *)(val);                           \
 305                                                                 \
 306         if (X##_e) X##_f0 |= _FP_IMPLBIT_E;                     \
 307         else X##_f0 &= ~(_FP_IMPLBIT_E);                        \
 308         _flo->bits.frac = X##_f0;                               \
 309         _flo->bits.exp  = X##_e;                                \
 310         _flo->bits.sign = X##_s;                                \
 311       }                                                         \
 312   } while (0)
 313
 314
 315 #define FP_UNPACK_E(X,val)              \
 316   do {                                  \
 317     FP_UNPACK_RAW_E(X,val);             \
 318     _FP_UNPACK_CANONICAL(E,2,X);        \
 319   } while (0)
 320
 321 #define FP_UNPACK_EP(X,val)             \
 322   do {                                  \
 323     FP_UNPACK_RAW_EP(X,val);            \
 324     _FP_UNPACK_CANONICAL(E,2,X);        \
 325   } while (0)
 326
 327 #define FP_PACK_E(val,X)                \
 328   do {                                  \
 329     _FP_PACK_CANONICAL(E,2,X);          \
 330     FP_PACK_RAW_E(val,X);               \
 331   } while (0)
 332
 333 #define FP_PACK_EP(val,X)               \
 334   do {                                  \
 335     _FP_PACK_CANONICAL(E,2,X);          \
 336     FP_PACK_RAW_EP(val,X);              \
 337   } while (0)
 338
 339 #define FP_ISSIGNAN_E(X)        _FP_ISSIGNAN(E,2,X)
 340 #define FP_NEG_E(R,X)           _FP_NEG(E,2,R,X)
 341 #define FP_ADD_E(R,X,Y)         _FP_ADD(E,2,R,X,Y)
 342 #define FP_SUB_E(R,X,Y)         _FP_SUB(E,2,R,X,Y)
 343 #define FP_MUL_E(R,X,Y)         _FP_MUL(E,2,R,X,Y)
 344 #define FP_DIV_E(R,X,Y)         _FP_DIV(E,2,R,X,Y)
 345 #define FP_SQRT_E(R,X)          _FP_SQRT(E,2,R,X)
 346
 347 /*
 348  * Square root algorithms:
 349  * We have just one right now, maybe Newton approximation
 350  * should be added for those machines where division is fast.
 351  * We optimize it by doing most of the calculations
 352  * in one UWtype registers instead of two, although we don't
 353  * have to.
 354  */
 355 #define _FP_SQRT_MEAT_E(R, S, T, X, q)                  \
 356   do {                                                  \
 357     q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1);         \
 358     _FP_FRAC_SRL_2(X, (_FP_WORKBITS));                  \
 359     while (q)                                           \
 360       {                                                 \
 361         T##_f0 = S##_f0 + q;                            \
 362         if (T##_f0 <= X##_f0)                           \
 363           {                                             \
 364             S##_f0 = T##_f0 + q;                        \
 365             X##_f0 -= T##_f0;                           \
 366             R##_f0 += q;                                \
 367           }                                             \
 368         _FP_FRAC_SLL_1(X, 1);                           \
 369         q >>= 1;                                        \
 370       }                                                 \
 371     _FP_FRAC_SLL_2(R, (_FP_WORKBITS));                  \
 372     if (X##_f0)                                         \
 373       {                                                 \
 374         if (S##_f0 < X##_f0)                            \
 375           R##_f0 |= _FP_WORK_ROUND;                     \
 376         R##_f0 |= _FP_WORK_STICKY;                      \
 377       }                                                 \
 378   } while (0)
 379
 380 #define FP_CMP_E(r,X,Y,un)      _FP_CMP(E,2,r,X,Y,un)
 381 #define FP_CMP_EQ_E(r,X,Y)      _FP_CMP_EQ(E,2,r,X,Y)
 382
 383 #define FP_TO_INT_E(r,X,rsz,rsg)        _FP_TO_INT(E,2,r,X,rsz,rsg)
 384 #define FP_FROM_INT_E(X,r,rs,rt)        _FP_FROM_INT(E,2,X,r,rs,rt)
 385
 386 #define _FP_FRAC_HIGH_E(X)      (X##_f1)
 387 #define _FP_FRAC_HIGH_RAW_E(X)  (X##_f0)
 388
 389 #endif /* not _FP_W_TYPE_SIZE < 64 */