soft-fp/extended.h

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